Improved guaranteed computable bounds on homogenized properties of periodic media by Fourier-Galerkin method with exact integration
IImproved guaranteed computable bounds on homogenizedproperties of periodic media by Fourier-Galerkin method withexact integration
Jaroslav Vondˇrejc
New Technologies for the Information Society, Faculty of Applied Sciences, Universityof West Bohemia, Univerzitn´ı 2732/8, 306 14 Plzeˇn, Czech Republic. E-mail: [email protected] Technische Universit¨at Braunschweig, Institute of Scientific Computing,M¨uhlenpfordstrasse 23, 38106 Braunschweig, Germany
Abstract
Moulinec and Suquet introduced FFT-based homogenization in 1994, and twenty yearslater, their approach is still effective for evaluating the homogenized properties arising fromthe periodic cell problem. This paper builds on the author’s (2013) variational reformulationapproximated by trigonometric polynomials establishing two numerical schemes: Galerkinapproximation (Ga) and a version with numerical integration (GaNi). The latter approach,fully equivalent to the original Moulinec-Suquet algorithm, was used to evaluate guaranteedupper-lower bounds on homogenized coefficients incorporating a closed-form double gridquadrature. Here, these concepts, based on the primal and the dual formulations, are em-ployed for the Ga scheme. For the same computational effort, the Ga outperforms the GaNiwith more accurate guaranteed bounds and more predictable numerical behaviors. Quadra-ture technique leading to block-sparse linear systems is extended here to materials definedvia high-resolution images in a way which allows for effective treatment using the FFT. Mem-ory demands are reduced by a reformulation of the double to the original grid scheme usingFFT shifts. Minimization of the bounds during iterations of conjugate gradients is effective,particularly when incorporating a solution from a coarser grid. The methodology presentedhere for the scalar linear elliptic problem could be extended to more complex frameworks.
Keywords:
Guaranteed bounds; Variational methods; Numerical homogenization; Galerkinapproximation; Trigonometric polynomials; Fourier Transform
This paper is devoted to FFT-based homogenization (Fourier-Galerkin method), a numericalmethod for evaluating homogenized (effective) material coefficients which are essential in multi-scale design. This method, which is an alternative to Finite Differences [1], Finite Elements [2, 3],Boundary Elements [4, 5], or Fast Multipole Methods [6, 7], Composite Finite Elements [8], X-FEM [9], or the Finite Cell Method [10], enables the direct treatment of material coefficientsdefined via high-resolution images.The method’s effectiveness relies on a Fast Fourier Transform (FFT) that is used for matrix-vector multiplication when solving linear systems; this can be provided in O ( N log N ) operationswhich outperforms most of the existing methods.The method’s reliability is provided by computable guaranteed upper-lower bounds on ho-mogenized properties [11], which are based on primal and dual variational formulations [12]along with a conforming approximation [13]. Surprisingly, this can be elegantly and efficiently1 a r X i v : . [ phy s i c s . c o m p - ph ] N ov reated also with the dual formulation since divergence-free fields can be constrained in theFourier domain, noticed e.g. by [14, 15]. Generally, the conforming approximation in the dualformulation is more difficult to provide because it requires special basis functions and techniquesespecially for vector-valued problems such as linearized elasticity.Up to now, the drawbacks of FFT-based methods compared to the above mentioned meth-ods included low adaptability originating from the use of regular discretization grids or indirectapplication for materials with holes. Difficulties also arose when treating various complex (non-linear, coupled) physical problems since the FFT-based method was originally formulated forthe Lippmann-Schwinger integral equation incorporating the Green function derived for a refer-ence medium as a parameter of the method. However, the method has already been applied toe.g. large deformations [16], viscoelasticity [17], thermo-elasticity [18], or fracture and damagemechanics [19]. For simplicity and clarity, the methodology is presented here only for a scalar linear ellipticproblem describing stationary heat transfer, electric conductivity, or diffusion.The model problem consists of an evaluation of homogenized properties A H ∈ R d × d thatcomply with the minimization problem A H E · E = inf u ∈ H , ( Y ) (cid:90) Y A ( x )[ E + ∇ u ( x )] · [ E + ∇ u ( x )] d x (1)for an arbitrary vector E ∈ R d . The region Y = (cid:81) dα =1 (cid:0) − , (cid:1) ⊂ R d accounts for a d -dimensionalcell where the material coefficients are defined through the bounded, symmetric, and uniformlyelliptic Y -periodic matrix function A : R d → R d × d . The trial space H , ( Y ) consists of the Y -periodic scalar functions u : Y → R with a square integrable gradient and zero mean.The original FFT-based method was proposed in 1994 by Moulinec and Suquet [20] as a newnumerical algorithm for solving the Lippmann-Schwinger equation, derived from (1). In thispaper, FFT-based methods build upon a variational reformulation by the author and co-workers[21, 13, 22] together with an approximation carried out with a truncated Fourier series space T d N with Fourier basis functions ϕ k ( x ) = exp(2 π i k · x ) having bounded frequencies k ∈ Z d N .Thus, the reference medium parameter is naturally avoided and two discretization schemes tothe cell problem (1) are revealed: Galerkin approximation (Ga) A H , N E · E = inf u N ∈ T N (cid:90) Y A ( x )[ E + ∇ u N ( x )] · [ E + ∇ u N ( x )] d x , (2a)and Galerkin approximation with numerical integration (GaNi) (cid:101) A GaNiH , N E · E = inf u N ∈ T N (cid:88) k ∈ Z d N A ( x kN )[ E + ∇ u N ( x kN )] · [ E + ∇ u N ( x kN )] , (2b)which is provided by a trapezoidal (rectangular) rule with integration points x kN for k ∈ Z d N located on a regular grid, introduced in section 3.1, particularly in (12) along with the index set Z d N in (11).The discretization of the cell problem (1) with trigonometric polynomials T d N is a standardnumerical approach, which has been used several times as the spectral Fourier-Galerkin or collo-cation method, e.g. [23, 14, 20, 24, 25, 26, 27], mostly within the Lippmann-Schwinger equation[28, 15] or in a standard variational setting [24]. A special attention is attributed to the uti-lization of the FFT algorithm in [23] or in [20] leading to the Moulinec and Suquet algorithmthat was interpreted as a trigonometric collocation method for the Lippmann-Schwinger equa-tion [29] or as a Fourier-Galerkin method with numerical integration (2b) in [13, section 5.2].2onnet in [15] incorporated not only FFT but also exact inclusion geometries into the Lippmann-Schwinger equation; it was interpreted in [13, section 4.2] as a Fourier-Galerkin method withexact integration (2a) in the standard variational framework.The theories regarding FFT-based methods — including discretization, convergence of ap-proximate solutions, and solution of the corresponding linear systems — have been providedonly recently by the author and co-workers [13] in the standard variational setting or later in[30] for rough material coefficients within the Lippmann-Schwinger equation.Another theoretically supported approach by Brisard and Dormieux in [31, 32] describesthe method using Galerkin approximation with piece-wise constant basis functions using aLippmann-Schwinger integral equation. However, the reference medium, a parameter of themethod, influences both the quality of the approximate solutions and the convergence of linearsolvers.Apart from the Fourier-Galerkin formulations (2) studied in [13, sections 4.2 and 4.3] and theBrisard and Dormieux approach [31, 32] — both of which lead to FFT-based schemes — thereare various other modifications and improvements. In [33], the authors provided a polarization-based scheme which can handle arbitrary phase contrast (voids and stiff inclusions); this can bealso managed by a numerical method based on augmented Lagrangians [34]. Recently, Willotet al. in [35, 36] adjusted the integral kernel in the Lippmann-Schwinger equation which led toimproved accuracy in approximate solutions, illustrated with a comparison using an analyticalsolution [37]. A possible improvement can also be achieved by smoothing of material coefficients[38],[11, section 8.3]. Significant attention was granted to improving the linear solvers leading tothe accelerated schemes [39, 18] or to the Krylov subspace methods, such as conjugate gradients[29, 31]; a comparison can be found in [40, 41]. The theory of guaranteed bounds on homogenized coefficients has been the subject of many stud-ies in analytical homogenization theories. These techniques employ the primal-dual formulationsof the cell problem (1) with limited — and often uncertain — information about the material co-efficients A . Specific examples include the Voigt [42], Reuss [43], and Hashin-Shtrikman bounds[44]; see the monographs [14, 45, 46, 47, 48] for a more complete overview. Because the boundsrely on limited data, their performance rapidly deteriorates for highly-contrasted media.Relatively less attention has been given to the computable upper-lower bounds arising froma conforming approximation to the cell problem (1). These bounds can be made arbitrarilyaccurate if the approximate solutions converge to the solution of (1). Moreover, the bounds areguaranteed if they allow for closed-form evaluation. Dvoˇr´ak and Haslinger, to our knowledge,specified the relevant ideas in their work [49, 50], which applied the approach to the p -versionof FEM. The application to the h -version of FEM occurs independently in [51] for an elastic-ity and to Fourier discretization within the Lippmann-Schwinger equation (Hashin-Shtrikmanfunctional) [28, 15] or the standard variational setting [24].The effective evaluation of upper-lower bounds using FFT has been provided recently in[52] for linear elasticity and later in [53] for permeability. Both frameworks rely on the Hashin-Shtrikman functional or the Lippmann-Schwinger equation, respectively, discretized with theBrisard–Dormieux method [32] assuming piecewise-constant approximations. However, the eval-uation of bounds is based on the summation of an infinite series leading to the loss of guarantee.The Fourier-Galerkin approaches used here in (2) or in [28, 24] provide guaranteed boundson homogenized properties, which is based on closed-form evaluation of corresponding bilinearforms in (2a) for material properties that have an analytical expression of Fourier coefficients.This approach has been already employed in [11] to GaNi (2b), while here the focus is onits generalization to and comparison with the Ga (2a). Moreover, these bounds can be made3rbitrarily accurate thanks to the convergence analysis of approximate solutions provided byVondˇrejc et al. [21, 13] and improved by Schneider [30] to account for rough coefficients.The guaranteed bounds incorporating trigonometric polynomials were mostly calculated withthe Hashin-Shtrikman functional [28]. The work [24] equivalent to (2a) is improved here byincorporating the FFT algorithm and by using a double grid quadrature, which leads to asparse structure according to [21, 11]. In [54], these ideas have been applied to linear elasticityillustrating that the classical variational formulation used here is always better than the Hashin-Shtrikman formulation. This paper is organized as follows. Notation and preliminaries to the periodic functions, Fouriertransform, and Helmholtz decomposition presented in section 2 are followed with the contin-uous homogenization problem in section 2.4. Then section 3 follows with Fourier-Galerkindiscretization for both primal and dual formulations leading to guaranteed bounds on homoge-nized coefficients; the structure between Ga and GaNi is established here. The methodology forevaluating guaranteed bounds is developed in section 4. Particularly in section 4.3, the doublegrid quadrature from [11] is generalized for materials defined via high-resolution images. Insection 4.4, the numerical scheme on the double grid is reduced to the original grid using shiftsof DFT (26). In section 4.5, the material properties without analytical expression of Fouriercoefficients are approximated in a way to still obtain the guaranteed bounds on the homogenizedproperties. Section 5 is dedicated to linear systems of Fourier-Galerkin schemes and the relatedcomputational aspects. Numerical examples in section 6 confirm the theoretical results andprovide a numerical comparison between the Ga (2a) and GaNi schemes (2b).
In the sections 2.1 and 2.2, the author introduces notation and recalls some useful facts relatedto matrix analysis and to spaces of periodic functions and the Fourier transform. Section 2.3 isdedicated to the Helmholtz decomposition of vector-valued periodic functions and its descriptionwith orthogonal projections, essential for the duality arguments in both discrete and continuoussettings.
In the subsequent section, d is reserved for the dimension of the model problem, assuming d = 2 ,
3. To keep the notation compact, X abbreviates the space of scalars, vectors, or matrices,i.e. R , R d , or R d × d , and ˆ X is used for their complex counterparts, i.e. C , C d , or C d × d . Vectorsand matrices are denoted by boldface letters, e.g. u , v ∈ R d or M ∈ R d × d , with Greek lettersused when referring to their entries; e.g. M = ( M αβ ) α,β =1 ,...,d . Matrix I = (cid:0) δ αβ (cid:1) αβ denotes theidentity matrix whereas the symbol δ αβ is reserved for the Kronecker delta, defined as δ αβ = 1for α = β and δ αβ = 0 otherwise.As usual, the matrix-vector product M v , the matrix-matrix product
M L , dot product u · v ,and the outer product u ⊗ v refer to( M u ) α = (cid:88) β M αβ u β , ( M L ) αβ = (cid:88) γ M αγ L γβ , u · v = (cid:88) α u α v α , ( u ⊗ v ) αβ = u α v β , where the author assumes that α and β range from 1 to d for the sake of brevity. Moreover, the4paces are endowed with the following inner product and norms, e.g. (cid:107) u (cid:107) C d = (cid:88) α u α u α , (cid:107) M (cid:107) C d × d = max u (cid:54) = (cid:107) M u (cid:107) C d (cid:107) u (cid:107) C d . For a unit cell Y = (cid:81) α (cid:0) − , (cid:1) , a function f : R d → X is Y -periodic if f ( x + k ) = f ( x ) for all x ∈ Y and all k ∈ Z d . According to [55, 56, 57], L p ( Y ; X ) = (cid:110) f : Y → X : f is Y -periodic, measurable, and (cid:107) f (cid:107) L p ( Y ; X ) < ∞ (cid:111) for p ∈ { , ∞} denotes the Lebesgue space equipped with the norm (cid:13)(cid:13) f (cid:13)(cid:13) L p ( Y ; X ) = ess sup x ∈Y (cid:13)(cid:13) f ( x ) (cid:13)(cid:13) X for p = ∞ , (cid:18)(cid:90) Y (cid:13)(cid:13) f ( x ) (cid:13)(cid:13) X d x (cid:19) / for p = 2 . The space L ( Y ; R d ) is also a Hilbert space with an inner product (cid:0) u , v (cid:1) L ( Y ; R d ) = (cid:90) Y u ( x ) · v ( x ) d x . For the sake of brevity, the author writes L p ( Y ) instead of L p ( Y ; R ), and often shortens L ( Y ; R d ) to L when referring to the norms and the inner product.Every function f ∈ L ( Y ; X ) can be expressed using Fourier series f ( x ) = (cid:88) k ∈ Z d (cid:98) f ( k ) ϕ k ( x )where the Fourier basis functions ϕ k and Fourier coefficients (cid:98) f ( k ) for k ∈ Z d are defined by ϕ k ( x ) = exp (cid:16) π i k · x (cid:17) (cid:98) f ( k ) = (cid:90) Y f ( x ) ϕ − k ( x ) d x ∈ (cid:98) X for x ∈ Y and k ∈ Z d , cf. [55, pp. 89–91]. The mean value of function f ∈ L ( Y ; X ) over periodic cell Y is denoted as (cid:104) f (cid:105) = (cid:90) Y f ( x ) d x = (cid:98) f ( ) ∈ X and corresponds to the zero-frequency Fourier coefficient. Operator ⊕ denotes the direct sum of mutually orthogonal subspaces, e.g. R d = U (1) ⊕ U (2) ⊕ . . . ⊕ U ( d ) for vectors U ( α ) = ( δ αβ ) β . According to the Helmholtz decomposition [56, pages 6–7], L ( Y ; R d ) admits an orthogonal decomposition L ( Y ; R d ) = U ⊕ E ⊕ J (3)into the subspaces of constant, zero-mean curl-free, and zero-mean divergence free fields U = { v ∈ L ( Y ; R d ) : v ( x ) = (cid:104) v (cid:105) for all x ∈ Y} , (4a) E = { v ∈ L ( Y ; R d ) : curl v = , (cid:104) v (cid:105) = } , (4b) J = { v ∈ L ( Y ; R d ) : div v = 0 , (cid:104) v (cid:105) = } . (4c)5ere, the differential operators curl and div are understood in the Fourier sense, so that(curl u ) αβ = (cid:88) k ∈ Z d π i (cid:0) k β (cid:98) u α ( k ) − k α (cid:98) u β ( k ) (cid:1) ϕ k , div u = (cid:88) k ∈ Z d π i k · (cid:98) u ( k ) ϕ k , cf. [56, pp. 2–3] and [57]. Furthermore, the constant functions from U are naturally identifiedwith vectors from R d .Alternatively, the subspaces arising in the Helmholtz decomposition (4) can be characterizedby the orthogonal projections introduced below. Definition 1.
Let G U , G E , and G J denote operators L ( Y ; R d ) → L ( Y ; R d ) defined via G • [ v ]( x ) = (cid:88) k ∈ Z d ˆΓ • ( k ) (cid:98) v ( k ) ϕ k ( x ) for • ∈ { U , E , J } , where the matrices of Fourier coefficients ˆΓ • ( k ) ∈ R d × d read ˆΓ U ( k ) = (cid:40) I ⊗ E ( k ) = (cid:40) ⊗ k ⊗ kk · k ˆΓ J ( k ) = (cid:40) ⊗ , for k = I − k ⊗ kk · k for k ∈ Z d \{ } . Lemma 2.
Operators G U , G E , and G J are mutually orthogonal projections with respect to theinner product on L ( Y ; R d ) , on U , E , and J . Here and in the next sections, the matrix field A : Y → R d × d spd is reserved for material coefficientswhich are required to be essentially bounded, symmetric, and uniformly elliptic: A ∈ L ∞ ( Y ; R d × d spd ) , c A (cid:13)(cid:13) v (cid:13)(cid:13) R d ≤ A ( x ) v · v ≤ C A (cid:13)(cid:13) v (cid:13)(cid:13) R d , (5)almost everywhere in Y and for all v ∈ R d with 0 < c A ≤ C A < + ∞ . By ρ A = C A /c A , theauthor denotes the contrast in coefficients A .Employing material coefficients, bilinear forms a, a − : L ( Y ; R d ) × L ( Y ; R d ) → R aredefined as a (cid:0) u , v (cid:1) := (cid:0) Au , v (cid:1) L ( Y ; R d ) a − (cid:0) u , v (cid:1) := (cid:0) A − u , v (cid:1) L ( Y ; R d ) . (6)Next, the author defines homogenization problem (1) in both primal and dual formulations. Definition 3 (Homogenized matrices) . Let material coefficients satisfy (5) . Then the primaland the dual homogenized matrices A H , B H ∈ R d × d satisfy (cid:0) A H E , E (cid:1) R d = min e ∈ E a (cid:0) E + e , E + e (cid:1) = a (cid:0) E + e ( E ) , E + e ( E ) (cid:1) , (7a) (cid:0) B H J , J (cid:1) R d = min ∈ J a − (cid:0) J + , J + (cid:1) = a − (cid:0) J + ( J ) , J + ( J ) (cid:1) (7b) for arbitrary quantities E , J ∈ R d and minimizers e ( E ) and ( J ) . Remark 4.
Notice that the primal formulation (7a) coincides with the problem (1) in theintroduction, because the subspace E from (3) admits an equivalent characterization E = {∇ f : f ∈ H ( Y , } . Remark 5 (Observations) . Thanks to the boundedness and ellipticity of material coefficients (5) , the homogenization problems in Definition 3 are well-posed with unique minimizers. Ho-mogenized matrices are symmetric, positive definite, and thus invertible in accordance with theperiodic homogenization theory, e.g. [58, 56, 59]. Moreover, they are mutually inverse A H = B − (8) in compliance with standard duality arguments [60]. Fourier-Galerkin discretization
This section describes the discretization of the homogenization problem (7) utilizing the Galerkinmethod in section 3.2 with an approximation space consisting of trigonometric polynomials insection 3.1. Then, the homogenized matrices defined with discrete problems are studied insection 3.3 to provide a structure of guaranteed bounds on homogenized properties.
Definition 6 (Trigonometric polynomials) . The number N ∈ N d such that N α is odd for all α (9) is an order of the approximation space of R d -valued trigonometric polynomials defined as T d N = (cid:110) (cid:88) k ∈ Z d N (cid:98) v k ϕ k : (cid:98) v k = (cid:98) v − k ∈ C d (cid:111) for ϕ k ( x ) = exp (2 π i k · x ) with x ∈ Y , (10) where the frequencies k are confined within the index set Z d N = (cid:110) k ∈ Z d : − N α ≤ k α < N α (cid:111) . (11)The approximation space (10) consists of real-valued functions only due to the prescribedHermitian symmetry of the Fourier coefficients, (cid:98) v k = (cid:98) v − k . Because of this and the necessity ofhaving a conforming approximation space in order to get guaranteed bounds on homogenizedproperties, the highest (Nyquist) frequencies k α = − N α / N according to (9).Next, I present the crucial property of trigonometric polynomials from space T d N : they canbe uniquely expressed on grid points x mM = (cid:88) α m α M α U ( α ) for M ∈ N d such that M α ≥ N α and m ∈ Z d M , (12)depicted in Figure 1(a). This can be shown with Discrete Fourier Transform (DFT) coefficientspresented here along with their orthogonality property ω knM = exp (cid:18) π i (cid:88) α k α n α M α (cid:19) , (cid:88) m ∈ Z d M ω − kmM ω mnM = | M | δ kn for k , n ∈ Z d M . (13)Therefore, a Fourier series of function u N ∈ T d N , noting that (cid:98) u N ( k ) = 0 for k ∈ Z d M \ Z d N and M ∈ N d such that M α ≥ N α , can be recast by substituting of the DFT coefficients (13) into u N ( x ) = (cid:88) k ∈ Z d M (cid:98) u N ( k ) ϕ k ( x ) = (cid:88) m ∈ Z d M (cid:88) n ∈ Z d M ω mnM (cid:98) u N ( n ) (cid:124) (cid:123)(cid:122) (cid:125) u N ( x mM ) (cid:88) k ∈ Z d M | M | ω − kmM ϕ k ( x ) (cid:124) (cid:123)(cid:122) (cid:125) ϕ M , m ( x ) , (14)which results in a fundamental trigonometric polynomial, see Figure 1(b), satisfying the Diracdelta property on grid points ϕ M , m ( x ) = 1 | M | (cid:88) k ∈ Z d M ω − kmM ϕ k ( x ) , ϕ M , m ( x nM ) = δ mn . (15)7his implies that the inverse DFT of Fourier coefficients (cid:98) u N ( n ) equals the function values ongrid points, i.e. (cid:80) n ∈ Z d M ω mnM (cid:98) u N ( n ) = u N ( x mM ). Thus, the trigonometric polynomial can beuniquely defined with the Fourier coefficients or with the values at grid points.In the following text, the letter N will be systematically used for odd approximation order(9) of trigonometric polynomials, while M will represent the number of grid points (12) withouta restriction to be odd. − . − . . . . Coordinate x − . − . . . . C oo r d i na t e x grid for N = [5 , grid for N = [15 , (a) Grid points − . − . . . . Coordinate x − . − . . . . F un c t i on v a l ue s real[ ϕ ]imag[ ϕ ] ϕ , grid points (b) Trigonometric polynomials
Figure 1: Examples of (a) two-dimensional grid points (12) and (b) one-dimensional complex-valued Fourier ϕ and real-valued fundamental ϕ , trigonometric polynomials (15) Generally, the difficulties in a conforming discretization of the continuous homogenization prob-lem (7) arise when approximating Helmholtz decomposition subspaces (3), especially in dualformulations for the space of divergence free fields J ; the approximation of curl-free space E in the primal formulation can be resolved with a conforming approximation of the space ofpotentials H , ( Y ).Luckily, in the case of trigonometric polynomials space (10), a Helmholtz decomposition T d N = U ⊕ E N ⊕ J N , E N = E (cid:92) T d N , J N = J (cid:92) T d N (16)can be easily performed with the orthogonal projections G U , G E , G J introduced in Definition 1for the continuous problem; the procedure is based on Lemma 2 and the property G • [ T d N ] ⊂ T d N observed directly from Definition 1. Moreover, these projections enable not only a properdiscretization but also an effective numerical treatment of discrete problems, cf. sections 5.1 and5.2.Then, the discrete homogenized problems are represented with two schemes, the Galerkinapproximation (Ga) and by its version with numerical integration (GaNi) corresponding tothe original Moulinec and Suquet scheme [20]. These Fourier-Galerkin schemes occurring in[23, 28, 24, 25, 26, 15, 27] as spectral methods were studied in the variational setting in [13,sections 4.2 and 4.3] for an odd approximation order (9) and in [11] for a general one. Definition 7 (Galerkin approximation — Ga) . The primal and the dual homogenization ma-trices A H , N , B H , N ∈ R d × d of Ga satisfy A H , N E · E = inf e N ∈ E N a (cid:0) E + e N , E + e N (cid:1) = a (cid:0) E + e ( E ) N , E + e ( E ) N (cid:1) (17a) B H , N J · J = inf N ∈ J N a − (cid:0) J + N , J + N (cid:1) = a − (cid:0) J + ( J ) N , J + ( J ) N (cid:1) (17b) for arbitrary quantities E , J ∈ R d . emark 8. Similarly to Remark 4 for a continuous setting, the primal formulation (17a) co-incides with the problem (2a) , because the subspace E N from (16) admits an equivalent charac-terization E N = {∇ f : f ∈ T N } . Evaluation of the integrals in (7), described in section 4, is generally unfeasible in closedform. This weakness can be compensated for with the trapezoidal (rectangular) integration ruleleading to the approximated bilinear forms (cid:101) a N , (cid:101) a − N : T d N × T d N → R expressed as a (cid:0) u N , v N (cid:1) ≈ (cid:101) a N (cid:0) u N , v N (cid:1) := (cid:88) k ∈ Z d N A ( x kN ) u N ( x kN ) v N ( x kN ) , (18a) a − (cid:0) u N , v N (cid:1) ≈ (cid:101) a − N (cid:0) u N , v N (cid:1) := (cid:88) k ∈ Z d N A − ( x kN ) u N ( x kN ) v N ( x kN ) . (18b)Noting that the objects related to this numerical integration are consistently denoted with thetilde symbol, e.g. (cid:101) a N , (cid:101) e ( E ) N , (cid:101) A H , N , or (cid:101) A boundH , N . Definition 9 (Galerkin approximation with numerical integration — GaNi) . Let material co-efficients (5) be additionally Riemann integrable or continuous A ∈ C ( Y ; R d × d spd ) . Then, theprimal and the dual homogenized coefficients (cid:101) A H , N , (cid:101) B H , N ∈ R d × d satisfy (cid:101) A H , N E · E = inf e N ∈ E N (cid:101) a N (cid:0) E + e N , E + e N (cid:1) = (cid:101) a N (cid:0) E + (cid:101) e ( E ) N , E + (cid:101) e ( E ) N (cid:1) , (19a) (cid:101) B H , N J · J = inf N ∈ J N (cid:101) a − N (cid:0) J + N , J + N (cid:1) = (cid:101) a − N (cid:0) J + (cid:101) ( J ) N , J + (cid:101) ( J ) N (cid:1) , (19b) for arbitrary quantities E , J ∈ R d . Remark 10 (Duality in GaNi) . Surprisingly, according to [11, Proposition 34], the duality ofhomogenized matrices from the continuous formulation (8) is inherited in the GaNi for the oddgrids (9) , particularly (cid:101) B − , N = (cid:101) A H , N . In [11, section 6], the minimizers of the GaNi scheme (cid:101) e ( E ) N , (cid:101) ( J ) N were used to a posteriori definethe approximate homogenized coefficients (cid:101) A boundH , N , (cid:101) B boundH , N ∈ R d × d satisfying for all E , J ∈ R d (cid:101) A boundH , N E · E = a (cid:0) E + (cid:101) e ( E ) N , E + (cid:101) e ( E ) N (cid:1) , (20a) (cid:101) B boundH , N J · J = a − (cid:0) J + (cid:101) ( J ) N , J + (cid:101) ( J ) N (cid:1) . (20b) In [11, section 6], we have shown that the matrices in (20) are guaranteed bounds on homogenizedcoefficients, i.e. (cid:0) (cid:101) B boundH , N (cid:1) − (cid:22) B − = A H (cid:22) (cid:101) A boundH , N . where the L¨owner partial order (cid:22) is used on a space of symmetric positive definite matrices R d × d spd according to [61, section 7.7] asfor L , M ∈ R d × d spd : L (cid:22) M if Lu · u ≤ M u · u for all u ∈ R d . Here, this structure is extended for the homogenized coefficients defined with the Ga scheme(17). 9 roposition 11 (Structure of guaranteed bounds) . The homogenized coefficients occurring inequations (7) , (17) , and (20) have the following structure (cid:0) (cid:101) B boundH , N (cid:1) − (cid:22) B − , N (cid:22) B − = A H (cid:22) A H , N (cid:22) (cid:101) A boundH , N . (21) Proof.
First of all, the proof of inequality A H (cid:22) A H , N is based on the conformity of approximatespace E N ⊂ E , recall (16). Indeed, as the minimization space of the homogenization problem(7a) is reduced, the minimum has to remain or increase, i.e. A H E · E = inf e ∈ E a (cid:0) E + e , E + e (cid:1) ≤ inf e N ∈ E N a (cid:0) E + e N , E + e N (cid:1) = A H , N E · E . Now, the following inequality A H , N (cid:22) (cid:101) A H , N is proven by substituting the minimizer e ( E ) N inthe GA scheme (17a) with the approximate minimizers (cid:101) e ( E ) N of the GaNi scheme (19a), i.e. A H , N E · E = inf e N ∈ E N a (cid:0) E + e N , E + e N (cid:1) ≤ a (cid:0) E + (cid:101) e ( E ) N , E + (cid:101) e ( E ) N (cid:1) = (cid:101) A boundH , N E · E . Since the primal inequalities A H (cid:22) A H , N (cid:22) (cid:101) A H , N are in hand, the dual formulation reveals theupper bounds on the dual matrix B H (cid:22) B H , N (cid:22) (cid:101) B boundH , N according to the same arguments. Theproof then arises from the inverse inequality, according to [61, Corollary 7.7.4.(a)], i.e.for L , M ∈ R d × d spd L (cid:22) M ⇐⇒ M − (cid:22) L − , (22)and from the duality of the continuous homogenization problem (8), i.e. B − = A H , whichfollows by the standard duality arguments in [60, 12] or [11, Proposition 7 and Corollary 9]. The computation of the guaranteed bounds on the homogenized coefficients (21) consists of theevaluation of bilinear forms (6) occurring in the formulation of the Ga scheme (17) or in the aposteriori estimate with the GaNi minimizers (20), i.e. integrals of the type (cid:0) Au N , v N (cid:1) L ( Y ; R d ) for A ∈ L ∞ ( Y ; R d × d ) and u N , v N ∈ T d N . (23)This integral evaluation based on double grid quadrature has already been analyzed in [11,section 6], which is summarized in section 4.2 for a matrix-inclusion composite (29). Then, insection 4.3, the methodology is generalized to an effective evaluation of integral (23) for grid-based composites, which are defined via high-resolution images assuming piece-wise constant orpiece-wise bilinear material coefficients.In section 4.4, I show that the double grid quadrature can be reduced to the original grid,resulting in a reduction of memory requirements. In section 4.5, I show that the guaranteedbounds on homogenized coefficients are not confined by a closed-form evaluation of (23) requiringthe knowledge of Fourier coefficient (cid:98) A ( k ); however, the approximation of A can be used, leadingto upper-upper and lower-lower bounds. A multi-index notation is systematically employed in which X M represents X M ×···× M d for M ∈ N d . Then the sets R d × M and (cid:2) R d × M (cid:3) or their complex counterparts represent the spaceof vectors and matrices denoted by bold serif font, e.g. u = (cid:0) u k α (cid:1) k ∈ Z d M α ∈ R d × M and U =( U km αβ ) k , m ∈ Z d M α,β ∈ (cid:2) R d × M (cid:3) with the index set Z d M defined in (11).10ub-vectors and sub-matrices are designated by superscripts, e.g. u k = ( u k α ) α ∈ R d or U km = ( U km αβ ) α,β ∈ R d × d . The scalar products on R d × M and C d × M are defined as (cid:0) u , v (cid:1) R d × M = 1 | M | (cid:88) k ∈ Z d M u k · v k , (cid:0) u , v (cid:1) C d × M = (cid:88) k ∈ Z d M u k · v k , where | M | = (cid:81) α M α stands for the number of discretization points. Moreover, the matrix-vectoror matrix-matrix multiplications follow from( Uv ) k = (cid:88) m ∈ Z d N U km v m ∈ R d or ( UV ) km = (cid:88) n ∈ Z d N U kn V nm ∈ R d × d , for k , m ∈ Z d N and V ∈ (cid:2) R d × M (cid:3) .The discretization operator that stores the values of a function at grid points (12) is definedby I M : C ( Y ; R d ) → R d × M , I M [ u ] = (cid:0) u α ( x mM ) (cid:1) m ∈ Z d M α =1 ,...,d . (24)When the polynomial of order N is expressed on a grid with the higher number of points M , u N = I M [ u N ] ∈ R d × M for u N ∈ T d N and M , N ∈ R d such that M α > N α , the discrete representation of the polynomial u N is kept with subscript N to emphasize thepolynomial order rather than the vector size. The actual dimension of u N is understood im-plicitly from the context, so that terms such as (cid:0) A M u N , u N (cid:1) R d × M with A M ∈ (cid:2) R d × M (cid:3) remainwell-defined.Using the discretization operator (24), the relation (14) between trigonometric polynomialsexpressed in Fourier and space domain can be represented in compact form as u N = (cid:88) k ∈ Z d N (cid:98) u kN ϕ k = (cid:88) m ∈ Z d M u mN ϕ M , m with u N = I M [ u N ] and (cid:98) u N = F N u N ∈ C d × N , (25)where matrices of vector-valued DFT and inverse DFT are defined as F N = 1 | N | (cid:0) δ αβ ω − mkN (cid:1) m , k ∈ Z d N α,β ∈ (cid:2) C d × N (cid:3) , F − N = (cid:0) δ αβ ω mkN (cid:1) m , k ∈ Z d N α,β ∈ (cid:2) C d × N (cid:3) . (26) Lemma 12.
For odd N according to (9) and M ∈ N d such that M α ≥ N α , discretizationoperator I M is an isometric isomorphism (one-to-one map that preserves distances) betweentrigonometric polynomials T d N defined in (10) and vector space I M [ T d N ] ⊆ R d × M , with equality I N [ T d N ] = R d × N for N = M . Here, the basic concepts regarding numerical integration of (23) are summarized according to[11]. The evaluation of bilinear forms in the conventional FEM leads to sparse matrices, whichalso is the case for the GaNi (19). However, more complicated scenarios arise with the Ga (17).A direct integration on the original grid leads to a fully populated matrix, see Lemma 14, whilea double-grid quadrature produces a sparse matrix with the same block diagonal structure asthe GaNi; see Lemma 15 and compare with Remark 13.
Remark 13 (Rectangular integration rule) . Using the notation about discrete spaces from sec-tion 4.1, the numerical integration in GaNi, recall (18) , leads to the following expression (cid:101) a N (cid:0) u N , v N (cid:1) = (cid:0)(cid:101) A N u N , v N (cid:1) R d × N for u N , v N ∈ T d N and u N = I N [ u N ] , v N = I N [ v N ] with (cid:101) A klN = δ kl A ( x kN ) for k , l ∈ Z d N . emma 14 (Fully populated expression) . According to [11, Lemma 37], the integral (23) equals (cid:0) Au N , v N (cid:1) L ( Y ; R d ) = (cid:0)(cid:98) A full N (cid:98) u N , (cid:98) v N (cid:1) C d × N = (cid:0) A full N u N , v N (cid:1) R d × N where vectors u N , v N ∈ R d × N and (cid:98) u N , (cid:98) v N ∈ C d × N are defined via u N = I N [ u N ] , v N = I N [ v N ] , (cid:98) u N = F N u N , (cid:98) v N = F N v N , and matrices A full N ∈ (cid:2) R d × N (cid:3) and (cid:98) A full N ∈ (cid:2) C d × N (cid:3) as (cid:0)(cid:98) A full N (cid:1) lk = (cid:90) Y A ( x ) ϕ k ( x ) ϕ − l ( x ) d x for l , k ∈ Z d N , A full N = F N (cid:98) A full N F − N . Lemma 15 (Sparse expression on a double grid, Lemma 39 in [11]) . For an odd order N according to (9) and M such that M α ≥ N α − , the integral (23) equals (cid:0) Au N , v N (cid:1) L ( Y ; R d ) = (cid:0) A M u N , v N (cid:1) R d × M where u N = I M [ u N ] , v N = I M [ v N ] ∈ R d × M , and A M ∈ (cid:2) R d × M (cid:3) with components ( A M ) kl αβ = δ kl (cid:88) n ∈ Z d N − ω knM (cid:98) A αβ ( n ) for k , l ∈ Z d M and α, β = 1 , . . . , d. (27)Because of the requirements in section 4.4, this lemma generalizes [11, Lemma 39] thatconfines M = 2 N − ; thus, a proof is presented here. Proof.
Because the product of two trigonometric polynomials u N v N = (cid:0) u N ,α v N ,α (cid:1) α ∈ T d N − has an order bounded by 2 N − , it can be expressed, in accordance with (14), on any grid M such that M α ≥ N α −
1, i.e. u N ,β v N ,α = (cid:88) m ∈ Z d M (cid:92) u N ,β v N ,α ( m ) ϕ m = (cid:88) m ∈ Z d M u N ,β ( x mM ) v N ,α ( x mM ) ϕ M , m where (cid:92) u N ,β v N ,α ( m ) = 0 for m ∈ Z d M \ Z d N − .Substitution into (23) and direct computation reveals (cid:0) Au N , v N (cid:1) L = (cid:88) α,β (cid:88) m ∈ Z d N − (cid:98) A αβ ( − m ) (cid:92) u N ,β v N ,α ( m )= (cid:88) α,β (cid:88) m ∈ Z d N − (cid:98) A αβ ( − m ) (cid:88) k ∈ Z d M ω − mkM | M | u N ,β ( x kM ) v N ,α ( x kM ) . The statement of the lemma follows by substituting n with − m . Remark 16 (Material coefficients leading to guaranteed bounds) . The closed-form evaluationof integral (23) leading to the fully discrete matrix (27) rests upon recognition of the Fourierseries expansion for material coefficients A . The space of functions having analytical expressionof Fourier coefficients constitutes a linear space thanks to the linearity of the integration. Threesuitable examples, which are also used for the numerical examples in section 6, are introducedhere. et h ∈ R d and r ∈ R be parameters such that < h α ≤ and r ≤ , then the periodicfunctions rect h , tri h , circ r ∈ L ∞ ( Y ; R ) are defined on Y along with their Fourier coefficients as rect h ( x ) = (cid:40) if | x α | < h α for all α otherwise , (cid:100) rect h ( k ) = (cid:89) α h α sinc ( h α k α ) , (28a)tri h ( x ) = (cid:89) α max { − | x α h α | , } , (cid:99) tri h ( k ) = (cid:89) α h α sinc ( h α k α ) , (28b)circ r ( x ) = (cid:40) for (cid:107) x (cid:107) < r otherwise , (cid:100) circ r ( k ) = (cid:40) πr for k = r B (2 πr (cid:107) k (cid:107) ) r (cid:107) k (cid:107) otherwise , (28c) where sinc( x ) = (cid:40) for x = 0 sin( πx ) πx for x (cid:54) = 0 and B is the Bessel function of the first kind. In [11, Lemma 38], the incorporation of characteristic functions (28) has been elaborated indetail for the inclusion-matrix composites, characterized by coefficients in the form A ( x ) = A (0) + n (cid:88) i =1 f ( i ) ( x − x ( i ) ) A ( i ) (29)where matrices A ( i ) ∈ R d × d for i = 0 , . . . , n represent the coefficients of the matrix phase andinclusions, functions f ( i ) ∈ L ∞ ( Y ) quantify the distribution of coefficients within inclusions,centered at x ( i ) , along with their topology. Their discrete coefficients (27) read as A klM = δ kl A (0) + n (cid:88) i =1 A ( i ) (cid:88) m ∈ Z d N − ω kmM ϕ − m ( x ( i ) ) (cid:98) f ( i ) ( m ) for k , l ∈ Z d M . (30) As an alternative to the inclusion-matrix composite (29), the grid-based composite is expressedas a linear combination of characteristic functions ψ concentrated on grid points x mP for m ∈ Z d P and P ∈ N d , i.e. A ( x ) = (cid:88) p ∈ Z d P ψ ( x − x pP ) C pP for P ∈ N d , x ∈ Y , and C P ∈ R d × d × P . (31)The material resolution is systematically expressed by the symbol P , which corresponds to thenumber of grid points and is independent of the order N ∈ N d of trigonometric polynomials T d N . Remark 17.
The material coefficients (31) are represented with a pixel- or voxel-based imagewith a resolution P ∈ N d . If the basis function ψ is taken as (28a) or (28b) , the coefficients arepiece-wise constant or bilinear, respectively. Lemma 18 (Grid-based composites) . The matrix (27) for material coefficients (31) is given by (cid:0) A M (cid:1) kl αβ = δ kl | P | (cid:88) m ∈ Z d N − ω kmM (cid:98) ψ ( m ) (cid:88) p ∈ Z d P ω − mpP | P | C pP ,αβ (32) for k , l ∈ Z d M and α, β = 1 , . . . , d , where (cid:98) ψ ( m ) for m ∈ Z d N − are the Fourier coefficients of ψ . roof. Using an affine substitution, the following formula for m , p ∈ Z d is deduced (cid:90) Y ψ ( x − x pP ) ϕ − m ( x ) d x = (cid:90) Y ψ ( x ) ϕ − m ( x + x pP ) d x = ω − mpP (cid:98) ψ ( m ) . This, (27), and (31) enable computation (cid:0) A M (cid:1) kl αβ = δ kl (cid:88) m ∈ Z d N − ω kmM (cid:90) Y A αβ ( x ) ϕ − m ( x ) d x = δ kl (cid:88) m ∈ Z d N − ω kmM (cid:88) p ∈ Z d P C pP ,αβ (cid:90) Y ψ ( x − x pP ) ϕ − m ( x ) d x = δ kl (cid:88) m ∈ Z d N − ω kmM (cid:88) p ∈ Z d P C pP ,αβ ω − mpP (cid:98) ψ ( m ) . In Lemma 15, the integral (23) occurring in the Ga scheme (17) has been evaluated on a doublegrid of size M such that M α ≥ N α −
1. Here in Lemma 19, it is reformulated using DFTshifts to the original size N occurring also in the GaNi scheme, cf. Remark 13. This procedurereduces memory requirements which is discussed in Remark 30. Lemma 19 (Reduction to the original grid) . Let S d = (cid:81) α { , } and A N be a matrix fromLemma 15. Then the fully populated matrices in Lemma 14 can be reformulated to A full N = 2 − d (cid:88) s ∈ S d F − N S ∗ N ( s ) F N A N ( s ) F − N S N ( s ) F N , (33a) (cid:98) A full N = 2 − d (cid:88) s ∈ S d S ∗ N ( s ) F N A N ( s ) F − N S N ( s ) , (33b) where the matrices S N ( s ) , A N ( s ) ∈ (cid:2) C d × N (cid:3) are defined for s ∈ S d as S N ( s ) = ( δ αβ δ mn ω − sm N ) m , n ∈ Z d N α,β , A N ( s ) = ( δ kr A (2 r − s )(2 r − s )2 N ,αβ ) k , r ∈ Z d N α,β . Proof.
According to Lemma 14 and 15, integral (23) can be expressed in three ways (cid:0) Au N , v N (cid:1) L = (cid:0) A full N u N , v N (cid:1) R d × N = (cid:0)(cid:98) A full N (cid:98) u N , (cid:98) v N (cid:1) C d × N = (cid:0) A N ˘ u N , ˘ v N (cid:1) R d × N (34)where u N = I − N [ u N ] = I − N [ F − N (cid:98) u N ] = I − N [˘ u N ] , v N = I − N [ v N ] = I − N [ F − N (cid:98) v N ] = I − N [˘ v N ] . Now, the matrix A N will be decomposed to meet (33). In order to reduce the double gridsizing 2 N to original grid sizing N , index k ∈ Z d N is uniquely split into k = 2 r − s with r ∈ Z d N and s ∈ S . (35)Using the connection of the two representations of trigonometric polynomials via DFT statedin (25), it holds for k ∈ Z d N according to (35) that˘ u kN = u N ( x k N ) = (cid:88) m ∈ Z d N ω km N (cid:98) u N ( m ) = (cid:88) m ∈ Z d N ω rmN ω − sm N (cid:98) u N ( m ) = (cid:2) F − N S N ( s ) (cid:98) u N (cid:3) r , Z d N instead of Z d N since the trigonometric polynomial u N belongs to T d N . Using the notation A k N := A kk N ∈ R d × d for k ∈ Z d N , the bilinear form (34)can be reformulated into (cid:0) A N ˘ u N , ˘ v N (cid:1) R d × N = 1 | N | (cid:88) k ∈ Z d N (cid:0) A k N ˘ u kN , ˘ v kN (cid:1) R d = 1 | N | (cid:88) α,β (cid:88) s ∈ S d (cid:88) r ∈ Z d N A r − s N ,αβ ˘ u r − sN ,β ˘ v r − sN ,α = 12 d | N | (cid:88) α,β (cid:88) s ∈ S d (cid:88) r , m , n ∈ Z d N A r − s N ,αβ (cid:0) ω rmN ω − sm N (cid:98) u mN ,β (cid:1) (cid:0) ω rnN ω − sn N (cid:98) v nN ,α (cid:1) . In the latter brackets, the sign of index n can be changed thanks to the symmetry of Z d N forodd grids (9), i.e. n ∈ Z d N ⇒ − n ∈ Z d N , which allow for expressing summations as a scalarproduct on C d × N , so (cid:0) A N ˘ u N , ˘ v N (cid:1) R d × N = 12 d | N | (cid:88) α,β (cid:88) s ∈ S d (cid:88) r , m , n ∈ Z d N A r − s N ,αβ (cid:0) ω rmN ω − sm N (cid:98) u mN ,β (cid:1) (cid:16) ω − rnN ω sn N (cid:98) v − nN ,α (cid:17) = 12 d (cid:88) α,β (cid:88) s ∈ S d (cid:88) r , m , n ∈ Z d N ω − sn N ω − rnN | N | A r − s N ,αβ (cid:0) ω rmN ω − sm N (cid:98) u mN ,β (cid:1) (cid:16)(cid:98) v nN ,α (cid:17) = 12 d (cid:88) s ∈ S d (cid:0) S N ( s ) F N A N ( s ) F − N S N ( s ) (cid:98) u N , (cid:98) v N (cid:1) C d × N where the conjugate symmetry of (cid:98) v nN and ω sn N has been used. Using substitution (cid:98) u N = F N u N and (cid:98) v N = F N v N , it can be reformulated as a scalar product on R d × N (cid:0) A N ˘ u N , ˘ v N (cid:1) R d × N = 12 d (cid:88) s ∈ S d (cid:0) F − N S N ( s ) F N A N ( s ) F − N S N ( s ) F N u N , v N (cid:1) R d × N . Throughout section 4, the author presents the methodology for evaluating guaranteed bounds onhomogenized properties relying on determination of the matrix (27). For closed-form evaluation,an analytical expression of Fourier coefficients (cid:98) A αβ ( m ) for m ∈ Z d N − is required, cf. Lemma 15.For general material coefficients, an approximate evaluation of the integral (23) can violatethe structure of guaranteed bounds (21); this can be resolved by appropriate adjustment ofthe material coefficients presented in the following Lemma 20 with a particular example inRemark 21. This approach is inspired by [49, 50], which incorporated outer approximation ofinclusion topology in the FEM framework. Lemma 20 (Upper-upper and lower-lower guaranteed bounds) . Let A , A ∈ L ∞ ( Y ; R d × d ) beupper and lower approximations of material coefficients (5) satisfying A ( x ) (cid:22) A ( x ) (cid:22) A ( x ) for almost all x ∈ Y , (36) and let a, a − : L ( Y ; R d ) × L ( Y ; R d ) → R be corresponding bilinear forms a (cid:0) u , v (cid:1) := (cid:0) Au , v (cid:1) L ( Y ; R d ) , a − (cid:0) u , v (cid:1) := (cid:0) A − u , v (cid:1) L ( Y ; R d ) . hen matrices A H , A H , A H , N , A H , N ∈ R d × d defined for arbitrary quantities E , J ∈ R d (cid:0) A H E , E (cid:1) R d = inf e ∈ E a (cid:0) E + e , E + e (cid:1) , (cid:0) B H J , J (cid:1) R d = inf ∈ J a − (cid:0) J + , J + (cid:1) , (cid:0) A H , N E , E (cid:1) R d = inf e ∈ E N a (cid:0) E + e , E + e (cid:1) , (cid:0) B H , N J , J (cid:1) R d = inf ∈ J N a − (cid:0) J + , J + (cid:1) , (cid:0) (cid:101) A H , N E , E (cid:1) R d = a (cid:0) E + (cid:101) e ( E ) N , E + (cid:101) e ( E ) N (cid:1) , (cid:0) (cid:101) B H , N J , J (cid:1) R d = a − (cid:0) J + (cid:101) ( J ) N , J + (cid:101) ( J ) N (cid:1) with minimizers (cid:101) e ( E ) N , (cid:101) ( J ) N of the GaNi scheme (19) comply with the following structure of guar-anteed bounds, i.e. (cid:101) B − , N (cid:22) B − , N (cid:22) B − (cid:22) B − = A H (cid:22) A H (cid:22) A H , N (cid:22) (cid:101) A H , N B − , N (cid:22) B − , N (cid:22) B − = A H (cid:22) A H , N (cid:22) A H , N . Proof.
It is possible to prove only the inequalities coming from the primal formulations sincethe dual part follows from the inverse inequality (22). The inequalities A H (cid:22) A H , N and A H (cid:22) A H , N (cid:22) (cid:101) A H , N have already been proven in Proposition 11 for material coefficients A and A ,respectively.In order to prove the rest, the following inequality is deduced for arbitrary v ∈ X ⊆ E inf e ∈X ⊆ E a (cid:0) E + e , E + e (cid:1) ≤ a (cid:0) E + v , E + v (cid:1) ≤ a (cid:0) E + v , E + v (cid:1) , where (36) and the monotonicity of the Lebesgue integration are used for the latter inequality.Since the first term is independent of v , it is possible to add an infimum, i.e. inf e ∈X ⊆ E a (cid:0) E + e , E + e (cid:1) ≤ inf e ∈X ⊆ E a (cid:0) E + e , E + e (cid:1) . The proof of A H (cid:22) A H and A H , N (cid:22) A H , N now followsfor the choices X = E and X = E N , respectively. Remark 21 (Choice of A and A ) . To comply with requirement (36) in the previous lemma, apossible choice of material coefficients consists of local approximations with piece-wise constantfunctions in a grid-based composite (31) . This material is then characterized with a pixel- orvoxel-based image defined via the following formula A = (cid:88) p ∈ Z d P rect h ( x − x pP ) C pP with C pP = (cid:13)(cid:13) A ( · − x pP ) (cid:13)(cid:13) L ∞ (Ω h ; R d × d ) · I , where vector P ∈ N d denotes an image resolution and where region Ω h = Π α (cid:0) − h α , h α (cid:1) for h α = P α represents a pixel or voxel placed at the origin with characteristic function rect h defined in (28a) . Factor (cid:13)(cid:13) A ( ·− x kP ) (cid:13)(cid:13) L ∞ (Ω h ; R d × d ) then indicates the largest eigenvalue of materialcoefficients A ( x ) ∈ R d × d over a pixel or voxel ( x pP + Ω h ) located at the corresponding grid point x pP . Remark 22.
The previous approximation of material coefficients, according to Lemma 20, leadsto guaranteed bounds. The following approximation with piece-wise bilinear functions A ( x ) ≈ (cid:88) p ∈ Z d P tri h ( x − x pP ) A ( x pP ) for h α = 1 P α enables the closed-form computation of bilinear forms; however, the guaranteed bounds are onlyapproximated. Linear systems and computational aspects
The focus of this section is on the resolution of minimizers defined by the Galerkin approxima-tions in section 3.2. Using the results about numerical integration in section 4, the linear systemsare described in section 5.2 with the help of discretization spaces of trigonometric polynomials,introduced in section 5.1. Then, computational aspects are discussed in section 5.3.
Definition 23 (Discrete spaces) . Using discretization operator I M and DFT matrix F definedin (24) and (26) respectively, the discrete spaces are introduced as U M = I M [ U ] , E N , M = I M [ E N ] , J N , M = I M [ J N ] , (37a) (cid:98) U N = F N [ U N ] , (cid:98) E N = F N [ E N , N ] , (cid:98) J N = F N [ J N , N ] . (37b)Thanks to the operator I M being an isometric isomorphism, see Lemma 12, Helmholtzdecomposition (16) for trigonometric polynomials is transformed to discrete spaces I N [ T d N ] = U N ⊕ E N , N ⊕ J N , N = R d × N , I M [ T d N ] = U M ⊕ E N , M ⊕ J N , M ⊂ R d × M . Now, the discrete projections on subspaces (37) are defined; for better orientation among oper-ators and subspaces, see following diagram (cid:98) E N F N ←− E N , N I N ←− E N I M −→ E N , M (cid:98) G E N , N − → G E N , N − → G E − → G E N , M − → F N [ R d × N ] F N ←− R d × N I N ←− T d N I M −→ I M [ T d N ] ⊆ R d × M (cid:98) G J N , N ← − G J N , N ← − G J ← − G J N , M ← − (cid:98) J N F N ←− J N , N I N ←− J N I M −→ J N , M . (38) Definition 24 (Discrete projections) . Let N ∈ R d satisfy the odd grid assumption (9) , M ∈ R d be a vector such that M α ≥ N α for all α . Then for • ∈ { U , E , J } , matrices (cid:98) G • N , M , G • N , M ∈ (cid:2) R d × M (cid:3) are defined as (cid:16)(cid:98) G • N , M (cid:17) kl = (cid:40) ˆΓ • ( k ) δ kl , for k , l ∈ Z d N for k , l ∈ Z d M \ Z d N , G • N , M = F − M (cid:98) G • N , M F M . where the matrices ˆΓ • ( k ) ∈ R d × d are Fourier coefficients of continuous projections introduced inDefinition 1, and where F M is the DFT matrix from (26) . Lemma 25 (Discrete projections) . For • ∈ { U , E , J } ,(i) operators (cid:98) G • N , N : C d × N → C d × N are orthogonal projections on (cid:98) U N , (cid:98) E N , and (cid:98) J N ,(ii) operators G • N , M : R d × M → R d × M are orthogonal projections on U M , E N , M , and J N , M .Proof. The fact that operators (cid:98) G • N , N and G • N , M are mutually orthogonal projections followsfrom direct calculation; the properties are inherited from continuous projections in Definition 1.The images of individual projections follow from the properties of operators I M , F N , G • alongwith the definition of subspaces (37) and (16). Indeed, the discrete projections can be expressedas (cid:98) G • N , N (cid:98) v N = F N ◦ I N ◦ G • ◦ I − N ◦ F − N (cid:98) v N , G • N , M v N = I M ◦ G • ◦ I − M v N for (cid:98) v N ∈ F N [ R d × N ] and v N ∈ I M [ T d N ]; see the structure of operators and subspaces (38).17 .2 Linear systems This section deals with resolutions of discrete minimizers from linear systems. This topic hasalready been studied in [13, section 5] and [11, section 7] for the GaNi scheme (19). Here, theconcept is summarized and extended to the Ga scheme (17).
Proposition 26 (From minimization to linear system) . Let H be a Hilbert space with a non-trivial orthogonal decomposition H = ˚ U ⊕ ˚ E ⊕ ˚ J , where ˚ U is isometrically isomorphic with R d .Next, let bilinear form ˚ a : H × H → R be defined as ˚ a (cid:0) u , v (cid:1) = (cid:0) ˚ Au , v (cid:1) H , for the symmetric,coercive, and bounded linear operator ˚ A : H → H , i.e. there exist c ˚ A > and C ˚ A > suchthat c ˚ A (cid:107) u (cid:107) H ≤ (cid:0) ˚ Au , u (cid:1) H ≤ C ˚ A (cid:107) u (cid:107) H for all u ∈ H . Then a problem for E ∈ ˚ U to find aminimizer ˚ e ( E ) ∈ ˚ E of ˚ e ( E ) = arg min ˚ e ∈ ˚ E ˚ a (cid:0) E + ˚ e , E + ˚ e (cid:1) (39a) is equivalent to finding the solution ˚ e ( E ) ∈ ˚ E of the following equation in H ˚ G ˚ A ˚ e ( E ) = − ˚ G ˚ AE , (39b) where ˚ G is an orthogonal projection on ˚ E .Proof. The proof starts with an optimality condition of (39a), namely ˚ a (cid:0) ˚ e ( E ) N , v (cid:1) = − ˚ a (cid:0) E , v (cid:1) for all v ∈ ˚ E . Then, the projection is incorporated in order to enlarge the space of test functions˚ a (cid:0) ˚ e ( E ) , ˚ Gv (cid:1) = − ˚ a (cid:0) E , ˚ Gv (cid:1) ∀ v ∈ H , (cid:0) ˚ G ˚ A ˚ e ( E ) N , v (cid:1) H = − (cid:0) ˚ G ˚ AE , v (cid:1) H ∀ v ∈ H , where the orthogonality (symmetry) of ˚ G has also been used. Now, it is possible to remove thescalar product and deduce the required (39b). Remark 27 (Linear systems for the GaNi) . According to [13, Proposition 12], with the notationfrom Remark 13, the minimizers I N [ (cid:101) e ( E ) N ] = (cid:101) e ( E ) N ∈ E N , N and I N [ (cid:101) ( J ) N ] = (cid:101) j ( J ) N ∈ J N , N for E , J ∈ R d in the GaNi scheme (19) satisfy the following equations G E N , N (cid:101) A N (cid:101) e ( E ) N = − G E N , N (cid:101) A N E , G J N , N (cid:101) A − N (cid:101) j ( J ) N = − G J N , N (cid:101) A − N J , (40) where G E N , N and G J N , N are projection matrices from Def. 24 and (cid:101) A kmN = δ km (cid:98) A ( k ) for k , m ∈ Z d N . Corollary 28 (Linear systems for the Ga) . Let G • N , M , (cid:98) G • N , M for • ∈ { E , J } be projectionmatrices from Definition 24. Then, for the minimizers e ( E ) N ∈ E N and ( J ) N ∈ J N of the Gascheme (17) for E , J ∈ R d , the following hold:(i) For M = 2 N − , the minimizers e ( E ) N := I M [ e ( E ) N ] ∈ E N , M and j ( J ) N := I M [ ( J ) N ] ∈ J N , M satisfy G E N , M A M e ( E ) N = − G E N , M A M E , G J N , M B M j ( J ) N = − G J N , M B M J , (41a) where A M , B M ∈ (cid:2) R d × M (cid:3) are defined in (27) , particularly in (30) and (32) for aninclusion-based (29) and grid-based (31) composite.(ii) The minimizers (cid:98) e ( E ) N := F N I N [ e ( E ) N ] ∈ (cid:98) E N and (cid:98) j ( J ) N := F N I N [ ( J ) N ] ∈ (cid:98) J N satisfy (cid:98) G E N , N (cid:98) A full N (cid:98) e ( E ) N = − (cid:98) G E N , N (cid:98) A full N (cid:98) E , (cid:98) G J N , N (cid:98) B full N (cid:98) j ( J ) N = − (cid:98) G J N , N (cid:98) B full N (cid:98) J , (41b) where matrices (cid:98) A full N , (cid:98) B full N ∈ (cid:2) C d × N (cid:3) are defined according to sparse decomposition (33b) . .3 Computational and implementation issues Here, practical aspects regarding the resolution of minimizers from linear systems are discussed.
Remark 29 (Solution by conjugate gradients) . The discrete problems, see the Ga (17) andGaNi (19) schemes, can be effectively solved with Krylov subspace methods [62, 63], particu-larly conjugate gradients [63, Algorithm 6.18]. It was pointed out in [29, 31] and explained byvariational reformulation in [22, 13] for the GaNi scheme (19) .Using the general notation from Proposition 26, the minimization problems of both discreteschemes (17) and (19) rely on the quadratic functional (39a) with a symmetric and positivedefinite matrix. Thus, conjugate gradients (CG) can be employed as the minimization oversubspace ˚ E is carried out with projection operator ˚ G .The minimization process also corresponds to the solution of the linear system (39b) with aninitial approximation ˚ e ( E )(0) from the minimization space ˚ E , which ensures that a residual vector r ( k ) = − ˚ G ˚ A ˚ e ( E )( k ) − ˚ G ˚ AE is from the subspace ˚ E for arbitrary k -th iteration. Then, the CG algorithm is interpreted as aminimization ˚ e ( E )( k ) = arg min ˚ e ∈ ˚ e ( E )(0) + K ( i ) ˚ a (cid:0) E + ˚ e , E + ˚ e (cid:1) (42) over Krylov subspaces defined for i = 1 , , . . . as K ( i ) = span (cid:110) r (0) , ˚ G ˚ A r (0) , . . . , (˚ G ˚ A ) i − r (0) (cid:111) satisfying K ( i ) ⊆ K ( i +1) ⊆ ˚ E = ˚ G [ H ] . The application of the CG algorithm only requires the implementation of the matrix-vector multi-plication of the linear system. For the GaNi (40) and the Ga (41) , it is outlined in Algorithms 1,2, and 3.
Algorithm 1
Matrix-vector multiplication for the primal formulation in the GaNi (40)
Require: A ∈ R d × d × N ← nonzero elements of (cid:101) A N ∈ (cid:2) R d × N (cid:3) (cid:46) A k α,β = A αβ ( x kN ) for α, β = 1 , . . . , d and k ∈ Z d N ; optionally (cid:98) G ∈ R d × d × N ← nonzeroelements of (cid:98) G E N , N ∈ (cid:2) R d × N (cid:3) procedure multiplication ( A , x ) (cid:46) calculates y = G E N , N (cid:101) A N x ∈ R d × N for x ∈ R d × N y ← (cid:101) A N x (cid:46) y k α = (cid:80) β A k αβ x k β for all α and k ∈ Z d N y ← F N y (cid:46) y α = FFT N ( y α ) for all α with an FFT of size N y ← (cid:98) G E N , N y (cid:46) y k α = (cid:80) β k α k β (cid:107) k (cid:107) R d y k β for all α and k ∈ Z d N y ← F − N y (cid:46) y α = iFFT N ( y α ) for all α with an inverse FFT of size N end procedure: return y Remark 30 (Memory and computational requirements) . For the approximation order N oftrigonometric polynomials, the linear systems in (40) and (41) have different sizes leading todifferent memory requirements, see Table 1. Noting that memory demands can be further re-duced by incorporating the symmetry of coef. matrix A and by calculating projection matrix (cid:98) G instead of storing, when needed. Despite the different sizes of linear systems, the number ofindependent unknowns remains the same and is equal to the dimension of approximation spaces E N , J N or their discrete relatives E N , N , J N , N . The subspace for primal formulations has a lgorithm 2 Matrix-vector multiplication for the primal formulation in the double grid Ga(41a)
Require: A ∈ R d × d × M ← nonzero elements of A M ∈ (cid:2) R d × M (cid:3) for M = 2 N − (cid:46) for anevaluation, see Remark 31; optionally (cid:98) G ∈ R d × d × N ← nonzero elements of (cid:98) G E N , M ∈ (cid:2) R d × M (cid:3) procedure multiplication ( A , x ) (cid:46) calculates y = G E N , M A M x ∈ R d × M for x ∈ R d × M y ← A M x (cid:46) y k α = (cid:80) β A k αβ x k β for α = 1 , . . . , d and k ∈ Z d M y ← F M y (cid:46) y α = FFT M ( y α ) for all α with an FFT of size M y ← (cid:98) G E N , M y (cid:46) y k α = (cid:80) β k α k β (cid:107) k (cid:107) R d y k β for k ∈ Z d N and y k α = 0 for k ∈ Z d M \ Z d N and all α y ← F − M y (cid:46) y α = iFFT M ( y α ) for all α with an inverse FFT of size M end procedure: return y Algorithm 3
Matrix-vector multiplication for the primal formulation in the reduced Ga (41b)
Require: A ∈ R d × d × N ← nonzero elements of A N ∈ (cid:2) R d × N (cid:3) (cid:46) for an evaluation, seeRemark 31; optionally (cid:98) G ∈ R d × d × N ← nonzero elements of (cid:98) G E N , N ∈ (cid:2) R d × N (cid:3) procedure multiplication ( A , x ) (cid:46) calculates y = (cid:98) G E N , N (cid:98) A full N x ∈ C d × N for x ∈ C d × N y ← ∈ C d × N for s ∈ S d = { , } d do z ← S N ( s ) x (cid:46) z k α = ω − sk N x k α for all α and k ∈ Z d N z ← F − N z (cid:46) z α = iFFT( z α ) for all α with an inverse FFT of size N z ← A N ( s ) z (cid:46) z k α = (cid:80) β A k − s αβ z k β for all α and k ∈ Z d N z ← F N z (cid:46) z α = FFT( z α ) for all α with an FFT of size N y ← y + 2 − d S ∗ N ( s ) z (cid:46) y k α = y k α + 2 − d ω sk N z k α for all α and k ∈ Z d N end for y ← (cid:98) G E N , N y (cid:46) y k α = (cid:80) β k α k β (cid:107) k (cid:107) R d y k β for all α and k ∈ Z d N end procedure: return y d | N | d | N − | ≈ d d | N | d | N | matrix of material coef. A d | N | d | N − | ≈ d d | N | d | N | = d d | N | projection matrix (cid:98) G d | N | d | N | d | N | indepen. unknowns in primal form. | N | − | N | − | N | − d − | N | −
1) ( d − | N | −
1) ( d − | N | − dimension dim E N = | N |− since this can be expressed using potential with zero-mean. Becausethe dimensions of constant fields and the whole space are dim U N = d and dim R d × N = d | N | ,the dimension of the dual space is equal to dim J N = ( d − | N | − .The linear systems for GaNi (40) and for Ga (41a) possess exactly the same mathematical struc-ture with a block-diagonal matrix of material coefficients (for isotropic material, only diagonal);see Algorithms 1 and 2; compare Remark 13 with Lemma 15. However, the Ga has a double sizeof vectors and matrices in the linear system. The corresponding higher memory and computa-tional requirements are outperformed with higher accuracy for the Ga scheme, see section 6.3for a comparison.In accordance with Lemma 19, the reduced Ga scheme (41b) benefits from the size reduction of anunknown vector, which is amplified when more vectors are stored (conjugate gradients, nonlinearproblems and solvers, etc.). Furthermore, the computational requirements remain approximatelythe same, which is illustrated in Figure 2. No. of dofs | N | − − − T i m e i n s e c ond s Ga, d = 2 Ga reduced, d = 2 Ga, d = 3 Ga reduced, d = 3 O ( | N | log | N | ) No. of dofs | N | . . . . . . . . T i m e r a t i o : G a r edu c ed vs . G a dimension d = 2 dimension d = 3 Figure 2: Comparison of matrix-vector multiplication for Ga (41a) with its reduced version(41b)
Remark 31 (Evaluation of material coefficients matrices with FFT) . The matrix (32) de-rived for the grid-based composite (31) can be evaluated efficiently using the FFT algorithm;for the inclusion-matrix composite (29) , the effective evaluation of (30) was discussed in [13,Remark 48].In (32) , the sum over Z d P and Z d N − is provided by P -sized FFT and M -sized inverseFFT algorithm resp., whereas the factor (cid:98) ψ ( m ) for m ∈ Z d N − occurs as an element-wisemultiplication. However, for P (cid:54) = 2 N − , the additional treatment has to be provided. For P α > N α − , the vector (cid:16) (cid:88) p ∈ Z d P ω − mpP | P | C pP ,αβ (cid:17) m ∈ Z d P ∈ C P s truncated to C N − , while for P α < N α − , it is periodically enlarged to C N − thanks tothe periodicity of ω · pP . This section is dedicated to numerical examples that confirm the properties of guaranteed bounds(21) with an emphasis on the comparison of Ga (17) with GaNi (19) and (20).
Problem 32.
A two-dimensional problem with material coefficients defined on a periodic cell Y = ( − , × ( − , ⊂ R is considered and defined via A ( x ) = I [1 + ρf • ( x )] for x ∈ Y , where I ∈ R × is the identity matrix, f : Y → R is a scalar nonnegative function which controlsthe shape of inclusions (recall Remark 16 for specific examples), and ρ > is a parametercorresponding to the phase contrast. Two types of inclusions, square and circle, are considered,namely f square ( x ) = (cid:40) for (cid:107) x (cid:107) ∞ < s otherwise , f circle ( x ) = (cid:40) for (cid:107) x (cid:107) < s otherwise , (43) where parameter s corresponds to an inclusion size, the side of the square and the radius,respectively. The problem is discretized with odd grids (9) with an example shown in Figure 3along with inclusion interfaces for both geometries (43) . − . − . . . . Coordinate x − . − . . . . C oo r d i na t e x inclusion interfacegrid for N = [5 , (a) Square inclusion − . − . . . . Coordinate x − . − . . . . C oo r d i na t e x inclusion interfacegrid for N = [5 , (b) Circle inclusion
Figure 3: Cells with grid and inclusion interfaces for size s = Remark 33.
All the computations have been provided using Python software
FFTHomPy avail-able at: https: // github. com/ vondrejc/ FFTHomPy. git . The linear systems presented insection 5.2 have been solved by conjugate gradients; a convergence criterion on the norm ofresiduum has been chosen with a relatively small tolerance − in order to suppress algebraicerror. The numerical examples are separated into the following parts: section 6.1 explores sensitiv-ity of homogenized properties in regard to inclusion size, section 6.2 describes an evolution ofupper-lower bounds for an increase in grid points, section 6.3 treats the behavior with differentphase contrasts, and section 6.4 shows the progress of guaranteed bounds during iterations ofconjugate gradients. 22 .1 Numerical sensitivity for the inclusion size
Here, homogenized properties are investigated with regard to an inclusion size s . Figure 4 depictsthe results for a relatively small number of discretization points N = (5 ,
5) which highlight thedifference between the Ga (17) and GaNi (19) schemes. . . . . . . Inclusion size (parameter s ) H o m ogen i z ed m a t r i x ( c o m ponen t ) A H , N B − , N e A boundH , N (cid:0) e B boundH , N (cid:1) − e A H , N (a) Square inclusion . . . . . . Inclusion size (parameter s ) H o m ogen i z ed m a t r i x ( c o m ponen t ) A H , N B − , N e A boundH , N (cid:0) e B boundH , N (cid:1) − e A H , N (b) Circle inclusion
Figure 4: Sensitivity of homogenized properties for ρ = 10 and N = (5 , A H , N , B − , N using the Ga scheme (17) when compared to the guaranteed bounds (cid:101) A boundH , N , (cid:0) (cid:101) B boundH , N (cid:1) − of the GaNi (20). The GaNi matrix (cid:101) A H , N in (19) together with its guaranteed bounds (20)has already been studied in [11], where the authors pointed out that the homogenized matrix (cid:101) A H , N , in some cases, underestimates or overestimates its own guaranteed bounds (cid:101) A boundH , N and (cid:0) (cid:101) B boundH , N (cid:1) − , respectively. The GaNi scheme (19) is influenced by inaccurate numerical integra-tion which disregards exact inclusion shapes because the scheme is defined only on grid points.As a result of exact integration, the homogenized matrices A H , N , B H , N change smoothly inrelation to the inclusion size s . This section is dedicated to the behavior of homogenized properties for an increase in the ap-proximation order of trigonometric polynomials N , see Definition 6. It is depicted in Figures 5and 6 for homogenized properties and also for their guaranteed errors defined as η N := tr (cid:18) A H , N − B − , N (cid:19) , (cid:101) η bound N := tr (cid:18) (cid:101) A boundH , N − (cid:0) (cid:101) B boundH , N (cid:1) − (cid:19) . (44)All the homogenized coefficients from both the Ga and GaNi schemes support the structure ofguaranteed bounds (21) and converge to homogenized matrix A H for an increasing number of gridpoints, which has been proven theoretically in [13, section 4.2] for the Ga scheme; the convergencefor GaNi is provided in [13, section 4.3] along with a regularization for discontinuous materialcoefficients according to [21, Section 3, pp. 115–117] or later in [30] for Riemann integrablecoefficients. Moreover, thanks to the hierarchy of approximation spaces E N ⊆ E M ⊂ E and J N ⊆ J M ⊂ J for N α ≤ M α , (45)the homogenized matrices of the Ga scheme A H , N , B − , N evolve monotonically as opposed to thehomogenized matrices of the GaNi scheme, which suffer, as already noticed in previous section,from inexact numerical integration causing the so-called ”variational crime” [64].23
20 40 60 80 100 120 140 160
No. of discretization points in each direction . . . . . . . H o m ogen i z ed m a t r i x ( c o m ponen t ) A H , N B − , N e A boundH , N (cid:0) e B boundH , N (cid:1) − e A H , N (a) Square inclusion
No. of discretization points in each direction . . . . . . . H o m ogen i z ed m a t r i x ( c o m ponen t ) A H , N B − , N e A boundH , N (cid:0) e B boundH , N (cid:1) − e A H , N (b) Circle inclusion
Figure 5: Bounds on homogenized matrix for inclusion size s = 0 . ρ = 10 No. of discretization points in each direction − − − E rr o r s i nho m ogen i z edp r ope r t i e s η N e η bound N rate = 1 (a) Square inclusion No. of discretization points in each direction − − − E rr o r s i nho m ogen i z edp r ope r t i e s η N e η bound N rate = 1 (b) Circle inclusion
Figure 6: Errors in homogenized properties (44) for inclusion size s = 0 . ρ = 10 This section investigates the homogenized properties in terms of normalized errors (44) foran increase in phase contrast ρ (see Figure 7). Moreover, it enables fair comparion of Ga withGaNi in terms of computational and memory requirements along with the accuracy of individualmethods. Indeed, the Ga (41a) and the GaNi (40) linear systems possess the same structure withblock-diagonal matrices of material coefficients; however, the Ga is evaluated on a double grid,resulting in higher computational and memory requirements for the same approximation order N ; see Remark 30 for a detailed discussion. Because of this, the GaNi is calculated with a doubleorder N than the Ga scheme; for this choice, the computational demands are approximatelythe same, while the memory requirements are even slightly lower for Ga, especially when thereduced version (41b) is used.Independently of inclusion shapes, the Ga (7) progresses with sharply better rates than theGaNi (19). Moreover, for the same computational demands, the Ga scheme produces tighterguaranteed bounds on homogenized properties, which is amplified for higher phase contrasts. Phase ratio ρ − − − N o r m a li z ede rr o r s η N ; N α = 203 e η bound N ; N α = 405 η N ; N α = 1823 e η bound N ; N α = 3645 rate = 1 rate = (a) Square inclusion Phase ratio ρ − − − − N o r m a li z ede rr o r s η N ; N α = 203 e η bound N ; N α = 405 η N ; N α = 1823 e η bound N ; N α = 3645 rate = 1 rate = (b) Circle inclusion
Figure 7: Normalized errors (44) for an increase in phase ratio ρ , s = 0 . Here, the author investigates the evolution of bounds during iterations of conjugate gradients(CG). In each iteration, a guaranteed bound is evaluated using the corresponding quadraticform as in (42). The results are shown in Figure 8 for primal formulation (upper bound), bothtopologies, and a relatively high phase contrast ρ = 10 to highlight the behavior.According to the standard results summarized in Remark 29, CG minimize the quadraticfunctional corresponding to the upper bound; the monotonic evolution of homogenized propertiesis confirmed in Figure 8. For all grid sizes, since the initial approximation for CG is taken as azero vector, the bounds begin from a Voigt bound (cid:104) A (cid:105) , the mean of material coefficients.25 Iterations of Conjugate gradients U ppe r bound s N α = 15 ; e N , (0) = N α = 45 ; e N , (0) = N α = 135 ; e N , (0) = N α = 45 ; e N , (0) = I N ◦ I − N / [ e N / ] N α = 135 ; e N , (0) = I N ◦ I − N / [ e N / ] (a) Square inclusion
Iterations of Conjugate gradients U ppe r bound s N α = 15 ; e N , (0) = N α = 45 ; e N , (0) = N α = 135 ; e N , (0) = N α = 45 ; e N , (0) = I N ◦ I − N / [ e N / ] N α = 135 ; e N , (0) = I N ◦ I − N / [ e N / ] (b) Circle inclusion
Figure 8: Upper bounds during iterations of conjugate gradients; ρ = 10 ; s = 0 . I N ◦ I − M : R d × M → R d × N is defined ongeneral grids N , M ∈ N d , N α > M α using discretization operator (24). The initial approxima-tion e ( E ) N , (0) = I N ◦ I − M [ e ( E ) M ] on a fine grid is then calculated from the solution of a linear systemon a coarse grid e ( E ) M ∈ R d × M with an FFT of size M and an inverse FFT of size N ; in thecase of Figure 8, coarse grid M is chosen to be N /
3. Note that no approximation is made inthis step because the corresponding trigonometric polynomial on the coarse grid equals the oneon the fine grid, i.e. I − M [ e ( E ) M ] = I − N [ e ( E ) N , (0) ] ∈ E M ⊂ E N . Here, the author shows how these methods can be applied to a complex material consisting ofalkali-activated fly ash foam. The coefficients, according to [65], A ( x ) = (cid:2) . · f ( x ) + 0 . · (cid:0) − f ( x ) (cid:1)(cid:3) · I for x ∈ Y , are defined via a fly ash phase characteristic function f : Y → R depicted in Figure 9 as avoxel-based image with resolution N = [99 , ,
99] corresponding to 970299 points. . . . Figure 9: A frontal view on a three-dimensional cellThe models were calculated on a conventional PC (Intel c (cid:13)
Core TM i7-4790 CPU @ 3.60GHzand 32 GB of RAM) within less than half an hour for both the GaNi and the Ga schemes. The26esults are represented for eigenvalues of homogenized coefficients because they also satisfy thestructure of upper-lower bounds (21), i.e. for the Ga scheme (17)eig A H , N = (cid:2) . . . (cid:3) , (46a)eig B − , N = (cid:2) . . . (cid:3) , (46b)for the GaNi scheme (19)eig (cid:101) A H , N = (cid:2) . . . (cid:3) , (47a)eig( (cid:101) B − , N ) = (cid:2) . . . (cid:3) , (47b)and for their corresponding guaranteed bounds (20)eig (cid:101) A boundH , N = (cid:2) . . . (cid:3) , (48a)eig (cid:0) (cid:101) B boundH , N (cid:1) − = (cid:2) . . . (cid:3) . (48b)The eigenvalues of the GaNi formulation (47) differ only because of an algebraic error and thisconfirms the duality of the GaNi scheme stated in [11, Propositin 34] (see Remark 10 for anoverview). Moreover, they are located between the guaranteed bounds obtained by both theGa (46) and the GaNi (48) schemes, and thus the GaNi provides an applicable prediction ofhomogenized properties. Because the guaranteed bounds comply with the energetic norms ofminimizers, the Ga (46) signifies a better approximation of local fields than the GaNi (48). Thisgap is accentuated in highly-contrasted media. This paper focuses on the numerical solution to the variational form of the unit cell problem (7),describing the homogenized properties of periodic heterogeneous materials. For discretization,two Fourier-Galerkin schemes were used and studied: Galerkin approximation (Ga) in (17) andits version with numerical integration (GaNi) in (19). In [11], the computable guaranteed boundson homogenized properties were introduced for the latter scheme. The approach, consisting inan exact evaluation of the primal-dual variational formulation for materials with an analyticalexpression of Fourier coefficients, is generalized here and applied to the Ga scheme, also resultingin a comparison with the GaNi. Theoretical results are confirmed with numerical examples. Tosummarize the most important findings: • The structure of the guaranteed bounds on homogenized properties, originating from Gaand GaNi, was established in Proposition 11, section 3.3. • In Lemma 18, section 4.3, the methodology for efficient double grid quadrature from theauthor’s previous work [11, section 6] is generalized for a grid-based composite (31). Thesematerial coefficients, defined via high-resolution images assuming e.g. piece-wise constantor bilinear approximation, can be effectively treated using FFT, see Remark 31. • Both the Ga (17) and GaNi (19) schemes lead to discrete formulations with a very similarblock-sparse structure; compare Remark 13 with Lemma 15 and linear system (40) with(41). However, the Ga is primarily evaluated on a double grid which can be recast tothe original grid using shifts of DFT, see Lemma 19. The memory and computationalrequirements discussed in Remark 30 are higher for the linear systems of Ga (41) thanthe GaNi (40). Nevertheless, the recast Ga (41b) leads to reduced memory requirementscompared to the original Ga (41a) without impacting computational costs involved insolving linear systems. 27
The Ga scheme (17) outperforms the GaNi (19). Under matching computational costs forboth schemes, the Ga provides more accurate results; guaranteed bounds on homogenizedproperties are more tight. The gap between the two schemes is accentuated in highly-contrasted media, section 6.3. • Evaluation of guaranteed bounds depends on knowing the Fourier coefficients of materialproperties A , see [11]. So, the approximation of A is proposed as a way to produceupper-upper and lower-lower guaranteed bounds, section 4.5. • Contrary to GaNi, the Ga scheme exhibits monotonous behavior, without oscillations inhomogenized properties, for an increase in grid points and for a change in inclusion size,sections 6.1 and 6.2. • Both schemes have the same rate of convergence of both minimizers and homogenizedproperties, which confirms the theoretical results in [13]. From the rate of convergence, itis possible to predict the grid size for a required level of accuracy, section 6.2. • The Ga scheme can be effectively solved using conjugate gradients providing monotonousimprovements of guaranteed bounds during iterations. Moreover, an approximate solu-tion on a coarse grid can be easily transferred to a fine grid to significantly improve theconvergence of the solution to the linear system, sections 5.2 and 6.4.To conclude, I recommend using the Ga scheme because it leads to more accurate approximationsfor the same computational effort. Moreover, the numerical behavior of the Ga is more smoothand predictable than the GaNi.The methodology used here is also valid for linearized elasticity. When using engineeringnotation (e.g. Mandel’s notation) in topological dimension 3, elasticity corresponds to a scalarproblem treated here for dimension d = 6 along with a different projection operator (cid:98) G . Never-theless, additional investigation is required for more complex problems. Acknowledgement
This work has been supported by project EXLIZ – CZ.1.07/2.3.00/30.0013 which is co-financedby the European Social Fund and the national budget of the Czech Republic and by the CzechScience Foundation through project No. P105/12/0331.
References [1] Flaherty JE, Keller JB. Elastic behavior of composite media.
Communications on Pureand Applied Mathematics (4):565–580, doi:10.1002/cpa.3160260409. URL http://doi.wiley.com/10.1002/cpa.3160260409 .[2] Guedes JM, Kikuchi N. Preprocessing and postprocessing for materials based on the ho-mogenization method with adaptive finite element methods. Computer Methods in AppliedMechanics and Engineering (2):143–198, doi:10.1016/0045-7825(90)90148-F.[3] Geers M, Kouznetsova V, Brekelmans W. Multi-scale computational homogenization:Trends and challenges. Journal of Computational and Applied Mathematics (7):2175–2182, doi:10.1016/j.cam.2009.08.077.[4] Eischen J, Torquato S. Determining elastic behavior of composites by the boundary elementmethod.
Journal of Applied Physics (1):159–170.285] Proch´azka P, ˇSejnoha J. A BEM formulation for homogenization of composites with ran-domly distributed fibers. Engineering analysis with boundary elements (2):137–144.[6] Greengard L, Lee J. Electrostatics and heat conduction in high contrast composite materi-als. Journal of Computational Physics (1):64–76.[7] Helsing J. The effective conductivity of arrays of squares: large random unit cells andextreme contrast ratios.
Journal of Computational Physics (20):7533–7547.[8] Hackbusch W, Sauter SA. Composite finite elements for the approximation of PDEs ondomains with complicated micro-structures.
Numerische Mathematik (4):447–472,doi:10.1007/s002110050248.[9] Legrain G, Cartraud P. An X-FEM and level set computational approach for image-1basedmodelling: Application to homogenization. International Journal for Numerical Methodsin Engineering (7):915–934.[10] D¨uster A, Sehlhorst HG, Rank E. Numerical homogenization of heterogeneous and cellularmaterials utilizing the finite cell method. Computational Mechanics jan 2012; (4):413–431, doi:10.1007/s00466-012-0681-2.[11] Vondˇrejc J, Zeman J, Marek I. Guaranteed upper-lower bounds on homogenized propertiesby FFT-based Galerkin method. Computer Methods in Applied Mechanics and Engineering :258–291, doi:10.1016/j.cma.2015.09.003.[12] Suquet P. Une m´ethode duale en homog´en´eisation: application aux milieux ´elastiques.
Journal de M´ecanique th´eorique et Appliqu´ee (Special issue)
Computers & Mathematics with Applications (3):156–173, doi:10.1016/j.camwa.2014.05.014.[14] Nemat-Nasser S, Hori M. Micromechanics: overall properties of heterogeneous materials .North-Holland: Amsterdam, 1993.[15] Bonnet G. Effective properties of elastic periodic composite media with fibers.
Journal ofthe Mechanics and Physics of Solids (5):881–899, doi:10.1016/j.jmps.2006.11.007.[16] Kabel M, B¨ohlke T, Schneider M. Efficient fixed point and Newton–Krylov solvers for FFT-based homogenization of elasticity at large deformations. Computational Mechanics
Cement and Concrete Research (2):197–207.[18] Vinogradov V, Milton GW. An accelerated FFT algorithm for thermoelastic and non-linear composites. International Journal for Numerical Methods in Engineering (11):1678–1695, doi:10.1002/nme.[19] Li J, Tian XX, Abdelmoula R. A damage model for crack prediction in brittle and quasi-brittle materials solved by the FFT method. International Journal of Fracture (2):135–146, doi:10.1007/s10704-011-9671-1.[20] Moulinec H, Suquet P. A fast numerical method for computing the linear and nonlinearmechanical properties of composites.
Comptes rendus de l’Acad´emie des sciences. S´erie II,M´ecanique, physique, chimie, astronomie (11):1417–1423.2921] Vondˇrejc J. FFT-based method for homogenization of periodic media: Theory and appli-cations. PhD Thesis, Czech Technical University in Prague, 2013. URL http://mech.fsv.cvut.cz/wiki/images/4/49/PhD_dissertation_Vondrejc_2013.pdf .[22] Vondˇrejc J, Zeman J, Marek I. Analysis of a Fast Fourier Transform Based Method forModeling of Heterogeneous Materials.
Large-Scale Scientific Computing , Lecture Notes inComputer Science , vol. 7116, Lirkov I, Margenov S, Wa´sniewski J (eds.). Springer: Berlin,Heidelberg, 2012; 512–522, doi:10.1007/978-3-642-29843-1 58.[23] Dykaar BB, Kitanidis PK. Determination of the effective hydraulic conductivity for hetero-geneous porous media using a numerical spectral approach: 1. Method.
Water ResourcesResearch (4):1155–1166, doi:10.1029/91WR03084.[24] Luciano R, Sacco E. Variational methods for the homogenization of periodic hetero-geneous media. European Journal of Mechanics - A/Solids (4):599–617, doi:10.1016/S0997-7538(99)80024-2.[25] Vainikko G. Fast solvers of the Lippmann-Schwinger equation. Direct and Inverse Problemsof Mathematical Physics :423–440.[26] Næss OF, Eckhoff KS. A Modified Fourier–Galerkin Method for the Poisson andHelmholtz Equations. Journal of Scientific Computing (1-4):529–539, doi:10.1023/A:1015162328151.[27] Cai H, Xu Y. A Fast Fourier–Galerkin Method for Solving Singular Boundary Integral Equa-tions. SIAM Journal on Numerical Analysis (4):1965–1984, doi:10.1137/070703478.[28] Nemat-Nasser S, Yu N, Hori M. Bounds and estimates of overall moduli of compositeswith periodic microstructure. Mechanics of Materials (3):163–181, doi:10.1016/0167-6636(93)90016-K.[29] Zeman J, Vondˇrejc J, Nov´ak J, Marek I. Accelerating a FFT-based solver for numericalhomogenization of periodic media by conjugate gradients. Journal of Computational Physics (21):8065–8071, doi:10.1016/j.jcp.2010.07.010.[30] Schneider M. Convergence of FFT-based homogenization for strongly heterogeneous media.
Mathematical Methods in the Applied Sciences (13):2761–2778, doi:10.1002/mma.3259.[31] Brisard S, Dormieux L. FFT-based methods for the mechanics of composites: A generalvariational framework. Computational Materials Science (3):663–671.[32] Brisard S, Dormieux L. Combining Galerkin approximation techniques with the principleof Hashin and Shtrikman to derive a new FFT-based numerical method for the homog-enization of composites. Computer Methods in Applied Mechanics and Engineering :197–212, doi:10.1016/j.cma.2012.01.003.[33] Monchiet V, Bonnet G. A polarization-based FFT iterative scheme for computing the ef-fective properties of elastic composites with arbitrary contrast.
International Journal forNumerical Methods in Engineering (11):1419–1436, doi:10.1002/nme.3295.[34] Michel JC, Moulinec H, Suquet P. A computational method based on augmented La-grangians and fast Fourier transforms for composites with high contrast. CMES: ComputerModeling in Engineering & Sciences (2):79–88.3035] Willot F, Abdallah B, Pellegrini YP. Fourier-based schemes with modified Green operatorfor computing the electrical response of heterogeneous media with accurate local fields. International Journal for Numerical Methods in Engineering (7):518–533, doi:10.1002/nme.4641.[36] Willot F. Fourier-based schemes for computing the mechanical response of composites withaccurate local fields. Comptes Rendus M´ecanique :232–245, doi:10.1016/j.crme.2014.12.005.[37] Craster RV, Obnosov YV. Four-phase checkerboard composites.
SIAM Journal on AppliedMathematics (6):1839–1856.[38] Merkert D, Andr¨a H, Kabel M, Schneider M, Simeon B. Voxel-based fast solution ofthe Lippmann-Schwinger equation with smooth material interfaces. Proceedings in AppliedMathematics and Mechanics , vol. 14, 2014; 579–580, doi:10.1002/pamm.201410277.[39] Eyre DJ, Milton GW. A fast numerical scheme for computing the response of compositesusing grid refinement.
The European Physical Journal Applied Physics (1):41–47.[40] Moulinec H, Silva F. Comparison of three accelerated FFT-based schemes for omputing themechanical response of composite materials. International Journal for Numerical Methodsin Engineering (13):960–985, doi:10.1002/nme.4614.[41] Mishra N, Vondˇrejc J, Zeman J. A comparative study on low-memory iterative solvers forFFT-based homogenization of periodic media. arXiv:1508.02045 , 2015.[42] Voigt W. Lehrbuch der kristallphysik , vol. 34. BG Teubner, 1910.[43] Reuss A. Berechnung der Fließgrenze von Mischkristallen auf Grund der Plas-tizit¨atsbedingung f¨ur Einkristalle.
ZAMM - Zeitschrift f¨ur Angewandte Mathematik undMechanik (1):49–58.[44] Hashin Z, Shtrikman S. A variational approach to the theory of the elastic behaviour ofmultiphase materials. Journal of the Mechanics and Physics of Solids (2):127–140.[45] Cherkaev A. Variational methods for structural optimization . Springer-Verlag: New York,2000.[46] Milton GW.
The Theory of Composites . Cambridge University Press: Cambridge, UK,2002.[47] Torquato S.
Random heterogeneous materials: microstructure and macroscopic properties ,Springer-Verlag: New York, 2002.[48] Dvorak GJ.
Micromechanics of Composite Materials . Springer: Netherlands, 2012.[49] Dvoˇr´ak J. Optimization of Composite Materials. Master’s Thesis, Charles University inPrague, 1993.[50] Haslinger J, Dvoˇr´ak J. Optimum composite material design.
RAIRO-Mathematical Mod-elling and Numerical Analysis-Modelisation Mathematique et Analyse Numerique (6):657–686.[51] Wi¸eckowski Z. Dual Finite Element Methods in Mechanics of Composite Materials. Journalof Theoretical and Applied Mechanics (33):233–252.3152] Kabel M, Andr¨a H. Fast numerical computation of precise bounds of effective elastic moduli. In: Berichte des Fraunhofer ITWM (224):1–16. URL http://math2market.de/Publications/2013ReportFraunhoferITWM_Nr224.pdf .[53] Bignonnet F, Dormieux L. FFT-based bounds on the permeability of complex microstruc-tures.
International Journal for Numerical and Analytical Methods in Geomechanics (16):1707–1723, doi:10.1002/nag.2278.[54] Monchiet V. Combining FFT methods and standard variational principles to computebounds and estimates for the properties of elastic composites. Computer Methods in AppliedMechanics and Engineering :454–473, doi:10.1016/j.cma.2014.10.005.[55] Rudin W.
Real and complex analysis . third edn., McGraw-Hill: New York, 1986.[56] Jikov VV, Kozlov SM, Oleinik OA.
Homogenization of Differential Operators and IntegralFunctionals . Springer-Verlag: Berlin, Heidelberg, 1994.[57] Saranen J, Vainikko G.
Periodic Integral and Pseudodifferential Equations with NumericalApproximation . Springer Monographs Mathematics: Berlin, Heidelberg, 2002.[58] Bensoussan A, Lions JL, Papanicolaou G.
Asymptotic Analysis for Periodic Structures .North Holland: Amsterdam, 1978.[59] Cioranescu D, Donato P.
An Introduction to Homogenization . Oxford Lecture Series inMathematics and Its Applications, Oxford University Press, 1999.[60] Ekeland I, T´emam R.
Convex Analysis and Variational Problems . North-Holland: Amster-dam, 1976.[61] Horn RA, Johnson CR.
Matrix analysis . Second edn., Cambridge University Press: NewYork, NY, USA, 2013.[62] Trefethen LN, Bau D.
Numerical linear algebra . SIAM: Philadelphia, PA, USA, 1997.[63] Saad Y.
Iterative Methods for Sparse Linear Systems . Second edn., SIAM: Philadelphia,PA, USA, 2003.[64] Strang G. Variational crimes in the finite element method.
The mathematical foundations ofthe finite element method with applications to partial differential equations
Journal of the EuropeanCeramic Society35