Development and comparison of spectral algorithms for numerical modeling of the quasi-static mechanical behavior of inhomogeneous materials
DDevelopment and comparison of spectral algorithms fornumerical modeling of the quasi-static mechanicalbehavior of inhomogeneous materials
M. Khorrami , J. R. Mianroodi , , P. Shanthraj , , B. Svendsen , ∗ Material Mechanics, RWTH Aachen University, Aachen, Germany Microstructure Physics and Alloy Design, Max-Planck Institute for Iron Research, Düsseldorf,Germany The School of Materials, University of Manchester, Manchester, U.K.
Abstract
In the current work, a number of algorithms are developed and compared for thenumerical solution of periodic (quasi-static) linear elastic mechanical boundary-value problems (BVPs) based on two di ff erent discretizations of Fourier series.The first is standard and based on the trapezoidal approximation of the Fouriermode integral, resulting in trapezoidal discretization (TD) of the truncated Fourierseries. Less standard is the second discretization based on piecewise-constant ap-proximation of the Fourier mode integrand, yielding a piecewise-constant dis-cretization (PCD) of this series. Employing these, fixed-point algorithms areformulated via Green-function preconditioning (GFP) and finite-di ff erence dis-cretization (of di ff erential operators; FDD). In particular, in the context of PCD,this includes an algorithm based on the so-called "discrete Green operator" (DGO)recently introduced by Eloh et al. (2019), which employs GFP, but not FDD. Forcomputational comparisons, the (classic) benchmark case of a cubic inclusion em-bedded in a matrix (e.g., Suquet, 1997; Willot, 2015) is employed. Both discon-tinuous and smooth transitions in elastic sti ff ness at the matrix-inclusion (MI)interface are considered. In the context of both TD and PCD, a number of GFP-and FDD-based algorithms are developed. Among these, one based on so-calledaveraged-forward-backward-di ff erencing (AFB) is shown to result in the greatestimprovement in convergence rate. As it turns out, AFB is equivalent to the "ro- ∗ Corresponding author
Preprint submitted to ... September 9, 2020 a r X i v : . [ c s . C E ] S e p ated scheme" (R) of Willot (2015) in the context of TD. In the context of PCD,comparison of the performance and convergence behavior of AFB / R- and DGO-based algorithms shows that the former is more e ffi cient than the latter for largerphase contrasts. Keywords: periodic matrix-inclusion boundary value problem, Fourier seriesdiscretization, Green-function preconditioning, finite-di ff erence discretization
1. Introduction
Numerical modeling of the mechanical behaviour of inhomogeneous materialsis often based on the assumption that their microstructure is periodic. In thiscase, a well-known alternative to finite di ff erence, finite element, or isogeometric,methods for the numerical solution of the corresponding boundary value problems(BVPs) is o ff ered by spectral methods (e.g., Gottlieb and Orszag, 1977; Canutoet al., 1988; Trefethen, 2000; Press et al., 2007; Kopriva, 2009). Besides the pos-sibility of exponential convergence, spectral methods generally o ff er much higherspatial resolution than finite-di ff erence- or finite-element-based ones. Indeed, asdiscussed for example by Press et al. (2007, §20.7), when applicable, these are themethods of choice for this purpose. Another di ff erence here is that, in contrast tofinite-di ff erence and finite-element methods, which approximate / discretize themodel field relations, spectral methods approximate / discretize their solution. Inparticular, in the Fourier case, this latter is based on discretization of truncatedFourier series, referred to as Fourier series discretization in what follows.For the current case of numerical solution of the (quasi-static) momentum balanceand corresponding mechanical BVPs, Fourier series discretization has tradition-ally been combined with di ff erential-operator preconditioning (e.g., Moulinec andSuquet, 1994; Suquet, 1997), formally analogous to di ff erential-operator inversionfor analytic solution of BVPs in mathematical physics based on Green functions(e.g., Lippmann-Schwinger). This is referred to as "Green-function precondition-ing" (GFP) in what follows. Together with continuous and discrete Fourier trans-formation, such methods have a long history of application in material scienceand in particular in continuum micromechanics (e.g., Khachaturyan, 1983; Mura,1987; Moulinec and Suquet, 1994; Suquet, 1997; Chen, 2002) to the modelingof material microstructure. The corresponding mechanical BVP is formulated atthe level of a unit cell or representative volume element. Traditionally, these arestrain-based, geometrically linear, and based on fixed-point iterative solution (e.g.,2uquet, 1997; Eyre and Milton, 1999; Michel et al., 2000, 2001; Lebensohn, 2001;Eisenlohr et al., 2013; Lebensohn and Rollett, 2020). More recently, these havebeen extended conjugate-gradient- and Newton-Krylov-based numerical solutionand geometric nonlinearity (e.g., Kabel et al., 2014; Shanthraj et al., 2015; Mishraet al., 2016).As well-known, in the case of material heterogeneity due to discontinuous mate-rial properties, Fourier series discretizationn su ff ers from oscillations due to theGibbs phenomena and to aliasing (e.g., Trefethen, 2000; Kopriva, 2009). In gen-eral, both have an adverse e ff ect on convergence behavior, resulting in a significantloss of algorithmic e ffi ciency and robustness. To address these issues, a number ofalgorithms going beyond the original "basic scheme" of Suquet (1997) have beendeveloped. One class of such algorithms combines Fourier series discretizationwith finite-di ff erence discretization (FDD) of di ff erential operators (e.g., Dreyerand Müller, 2000; Brown et al., 2002; Moulinec and Silva, 2014; Willot, 2015),resulting in a reduction of the e ff ect of oscillations on convergence behavior due tomodes at wavelengths shorter than the nodal or grid spacing. As shown in particu-lar by Willot (2015), among such "accelerated schemes" (e.g., Moulinec and Silva,2014), his "rotated scheme" is most e ff ective in improving algorithmic e ffi ciencyand robustness. This scheme has been employed for numerical solution of BVPsfor heterogeneous materials in a number of works (e.g., for field dislocation me-chanics in Djaka et al., 2017). More recently, Schneider et al. (2017) showed thatfinite-element discretization based on linear hexahedral elements with reduced in-tegration is equivalent to the rotated scheme of Willot (2015). Quite recently,Lucarini and Segurado (2019) developed a displacement-based approach via di-rect Fourier transformation and real-function-based Fourier transform symmetryreduction, yielding an algorithmic system based on a full-rank associated Her-mitian matrix amenable to preconditioning. For the latter purpose, they workedwith a Green function based on the average elastic sti ff ness. In contrast to Willot(2015) and the current work, FDD was not employed.All algorithms or "schemes" discussed up to this point employ Fourier seriesdiscretization in "standard" form, i.e., based on trapezoidal approximation / dis-cretization (TD) of the integral for Fourier modes with respect to the unit cell. TDis the basis of Fourier interpolation and the standard form of the discrete Fouriertransform (e.g., Trefethen, 2000; Press et al., 2007; Kopriva, 2009). Alternatively,Brisard and Dormieux (2010) and Anglin et al. (2014) discretize by assumingthe solution is piecewise-constant, resulting in piecewise-constant discretization(PCD). This approach has been exploited recently by Eloh et al. (2019), who3ombine PCD with GFP to obtain an algorithm depending on a so-called discreteGreen operator (DGO) for numerical solution of periodic mechanical BVPs formaterially heterogeneous elastostatics via fixed-point interation.One purpose of the current work is to develop additional algorithms in the contextof Fourier series discretization based on TD and PCD. All of these exploit GFP,and some of them FDD. A second purpose is to carry out detailed theoretical andcomputational comparisons of these with existing algorithms, in particular withthe "rotated" (R) scheme of Willot (2015), and with the DGO-based scheme ofEloh et al. (2019). The corresponding computational comparisons are based onthe classic benchmark case of a periodic matrix-inclusion microstructure with cu-bic inclusions having sharp corners. The material inhomogeneity takes the form of(both discontinuous and smooth) phase contrast in elastic sti ff ness (e.g., Suquet,1997; Nemat-Nasser and Hori, 1999; Willot, 2015) across the matrix-inclusion(MI) interface. In contrast, Eloh et al. (2019) focus mainly on material hetero-geneity due to eigenstrains (see also e.g., Anglin et al., 2014); in one example,however, spherical inclusions are considered. Likewise, Lucarini and Segurado(2019) focus (solely) on spherical inclusions.To establish relevant basic concepts and issues, the work begins in Section 2 withbasic considerations in one dimension (1D). These include analytic solution of theperiodic matrix-inclusion mechanical BVP and its truncated Fourier series repre-sentation. Attention is focused next on the development of algorithms in 1D basedon TD and PCD of truncated Fourier series in Section 3. As discussed above, allof these employ GFP, and some of them FDD. Computational comparison of theseis then carried out in Section 4 for the matrix-inclusion benchmark case. In partic-ular, material heterogeneity in the form of both (i) discontinuous and (ii) smoothcompliance / sti ff ness distributions are considered. Multidimensional forms of the1D algorithms are obtained via direct tensor product generalization in Section 5.Corresponding computational comparsions are presented in Section 6. The workends with a summary and discussion in Section 7.In this work, (three-dimensional) Euclidean vectors are represented by lower-casebold italic characters a , b , . . . , . In particular, let i , i , i represent the Cartesianbasis vectors. Second-order tensors are represented by upper-case bold italic char-acters A , B , . . . , with I the second-order identity. By definition, any A maps any b linearly into a vector Ab . Fourth-order Euclidean tensors A , B , . . . , are de-noted by upper-case slanted sans-serif characters. By definition, any A maps any B linearly into a second-order tensor A B . Let A · B = A i jk ... B i jk ... (summationconvention) represent the scalar product of two arbitrary-order tensors A and B .4sing this, A T b · c : = b · Ac defines the transpose A T of A . Let sym A : = ( A + A T )represent the symmetric part of A . The tensor products ( A (cid:3) B ) C : = ACB and( A (cid:52) B ) C : = AC T B will also be employed. Additional notation will be introducedas needed in what follows.
2. Basic considerations in 1D
Let the interval [0 , l ] of length l > f = ¯ f + ˜ f , ¯ f : = l (cid:90) l f ( x ) dx , ˜ f : = f − ¯ f , (1)of any (integrable) f into mean ¯ f and "fluctuating" ˜ f parts.Let u represent the displacement field on the unit cell and H = ∇ u the correspond-ing distortion. Integration of ∇ u = ¯ H + ˜ H yields the split u = u h + u p , u h ( x ) : = c + ¯ H x , u p ( x ) : = (cid:90) ˜ H ( x ) dx , (2)of u into "homogeneous" u h and "particular" u p parts, with c the constant of in-tegration. In what follows, c and ¯ H are assumed given or known. Then u h ( x ) isdetermined, and u p ( x ) represents the primary unknown field.Excluding cracking, kinematic compatibility requires u , and so u p , to be continu-ous everywhere, in particular at MI interfaces. In addition, quasi-static mechanicalequilibrium div T = T to be continuous everywhere and infact constant in 1D; then ˜ T = T = ¯ T . In this case, the linear elastic relation E = S T = C − T for the 1D strain E ≡ H in terms of the compliance S = C − andsti ff ness C results in the equilibrium relations ¯ E = ¯ S ¯ T and ˜ E = ˜ S ¯ T for E = ¯ E + ˜ E in 1D. In turn, ˜ E ( x )¯ E = ˜ S ( x )¯ S , u p ( x )¯ E = (cid:90) ˜ S ( x )¯ S dx , (3)then hold for ˜ E ( x ) and u p ( x ), the latter via (2) .In the classic composite case (e.g., Suquet, 1997), the interface between two bulkphases is idealized as "sharp" in the sense that material properties like S are as-sumed to vary discontinuously across the interface. On the other hand, phase field5odels idealize such interfaces as a mixture of the bulk phases in which propertieslike S vary smoothly from one bulk phase to the other. Assuming the classic caseof a double well potential and "flat" interface region of half-width (cid:15) centered at x = c , solution of the equilibrium Ginzburg-Landau relation at zero stress yields φ (cid:15) ( x ; c ) = + tanh(( x − c ) /(cid:15) ) (4)for the phase field of the MI system varying between φ (cid:15) = φ (cid:15) = φ ( x ; c ) : = lim (cid:15) → φ (cid:15) ( x ; c ) = θ ( x − c ) in terms of the modified Heaviside stepfunction θ ( x ) (i.e., θ ( x < = θ ( x =
0) : = / θ ( x > =
1; e.g., Bracewell,2000, Chapter 4). Assuming that the left (right) MI interface is at x = c l ( x = c r ), ν (cid:15) ( x ) : = φ ( x ; c l , (cid:15) ) − φ ( x ; c r , (cid:15) ) (5)represents the 1D "volume density" of the inclusion. In terms of this, S = (1 − ν (cid:15) ) S M + ν (cid:15) S I = S M + S IM ν (cid:15) , S IM : = S I − S M , (6)holds for the elastic compliance of the MI unit cell (note S M / S I = C I / C M ). Sub-stituting (6) into (3) yields˜ E ( x ) = s IM ˜ ν (cid:15) ( x ) ¯ E , u p ( x ) = s IM [ λ (cid:15) ( x ) − λ (cid:15) (0) − ¯ ν (cid:15) x / l ] ¯ El , s IM : = S IM / ¯ S , (7)with λ (cid:15) ( x ) : = [ (cid:96) (cid:15) ( x ; c l ) − (cid:96) (cid:15) ( x ; c r )] / l and (cid:96) (cid:15) ( x ; c ) : = (cid:82) φ (cid:15) ( x ; c ) dx . Note that u p (0) = without loss of physical generality. ˜ E ( x ) and u p ( x ) aredisplayed in Figure 1. (cid:1) (cid:1)(cid:2) (cid:1)(cid:3) - (cid:2)(cid:1)(cid:2) (cid:1) / (cid:2)(cid:3) (cid:3) (cid:1) (cid:1)(cid:2) (cid:1)(cid:3) (cid:1) (cid:1)(cid:2) (cid:1) / (cid:2)(cid:3) (cid:4) (cid:4) (cid:2) Figure 1: ˜ E / ¯ E (left) and u p / ¯ El (right) from (7) on the left half [0 , l /
2] of the unit cell [0 , l ] forsharp (red: (cid:15)/ l =
0) and smooth (green: (cid:15)/ l = / (cid:15)/ l = /
50) MI interfaces. Theseand all 1D results to follow are based on c l = l / c r = l /
4, and s IM = − /
100 (i.e., C I / C M = S M / S I =
6s expected from (3), ˜ S determines a continuous and smooth ˜ E for (cid:15) >
0, and adiscontinuous ˜ E for (cid:15) =
0. Likewise, u p is continuous and smooth for (cid:15) >
0, butonly continuous for (cid:15) = As a first step toward Fourier-series-based numerical solution of the above BVP,consider next the truncated Fourier series approximation of ˜ E and u p from (7).The truncated Fourier series of any f ∈ L (0 , l ) is given by F m f ( x ) : = (cid:88) m κ = − m e ı k κ x ˆ f κ , ˆ f κ : = l (cid:90) l e − ı k κ x f ( x ) dx , (8)(e.g., Trefethen, 2000; Kopriva, 2009; Liu, 2011) with ı : = √− k κ : = πκ/ l .As is well-known, F ∞ f ( x ) (cid:44) f ( x ) in general at one or more x ∈ [0 , l ]. In particular,for f discontinuous on [0 , l ], then, F m f deviates from f due to (i) F ∞ f ( x ) (cid:44) f ( x )at points of discontinuity x ∈ [0 , l ], and (ii) truncation of F ∞ f . Since (i) is relatedto the Gibbs phenomenon, it will be referred to as Gibbs error in what follows. Inthe sharp ( (cid:15) =
0) MI interface case, then, F m ˜ E is a ff ected by (i) and (ii), whereas F m u p is a ff ected only by (ii). These are displayed in Figure 2. (cid:1) (cid:1)(cid:2) (cid:1)(cid:3) - (cid:2)(cid:1)(cid:2) (cid:1) / (cid:2)(cid:3) (cid:1) (cid:4) (cid:4) (cid:1) (cid:1)(cid:2) (cid:1)(cid:3) (cid:1) (cid:1)(cid:2) (cid:1) / (cid:2)(cid:3) (cid:1) (cid:4) (cid:4) (cid:5) (cid:2) Figure 2: Normalized Fourier series approximation F m ˜ E ( x ) / ¯ E (left) and F m u p ( x ) / ¯ El (right) of ˜ E ( x )and u p ( x ), respectively, on the left half unit cell [0 , l /
2] for discontinuous ( (cid:15) =
0) compliance anddi ff erent degrees of truncation m < ∞ . In particular, m = m =
15 (green), and m = E ( x ) / ¯ E and u p ( x ) / ¯ El from (7) are shown(black) for comparison. As expected, lim m →∞ F m ˜ E ( x ) (cid:44) ˜ E ( x ) at x = l /
4; rather, lim m →∞ F m ˜ E ( l / = lim x ↑ l / ˜ E ( x ) + lim x ↓ l / ˜ E ( x ). Being due to approximation of discontinuous˜ E by F ∞ ˜ E , note that Gibbs error is una ff ected by truncation of F ∞ ˜ E to F m ˜ E , i.e.,independent of m . In addition, it is independent of the magnitude of S IM .7 . Numerical solution algorithms in 1D As discussed in the introduction, algorithms for the numerical solution of the cur-rent BVP are based on (i) discretization of (8) in 1D, and (ii) tensor-product gen-eralization of these to 3D. More specifically, two discretizations of (8) are con-sidered and compared in this work. The most common of these is trapezoidaldiscretization of (8) based on that [0 , l ] = (cid:83) ni = [ x i , x i + ] of [0 , l ] with nodes at x i = ( i − h and nodal spacing h = l / n . In this case, (8) and ˆ f κ discretizes toˆ f t κ : = n − / ˇ f t κ , ˇ f t κ : = n − / (cid:88) ni = e − ı k κ x i f t i , (9)("t" stands for "trapezoidal"). As detailed in Appendix A, (9) induces the trape-zoidal discretization (TD) f T ( x ) : = ϕ t µ ( x ) ˇ f t µ , ϕ t µ ( x ) : = n − / e q µ x , f t i = F t + µ i ˇ f t µ , ˇ f t µ = F t − µ i f t i , (10)of (8) via (A.6) for i , µ = , . . . , n . Here, ˇ f t µ ≡ ˇ f t κ = µ − − m , q µ : = ı k µ − − m , f t i ≡ f ( x i ), and F t ± µ i : = n − / e ± q µ x i from (A.7). If ¯ f is prescribed, note that ˇ f t µ = m + = ˇ f t κ = = n − / (cid:80) ni = f t i = n / ¯ f is determined for a given discretization. For notationalsimplicity, the "t" superscript for "trapezoidal" is dropped in what follows.As discussed in the introduction, following previous work (e.g., Dreyer and Müller,2000; Moulinec and Silva, 2014; Willot, 2015), finite di ff erence discretization(FDD) of di ff erential operators in the context of (10) is employed here for im-proved robustness in convergence. For example, u pT ( x ) = ϕ µ ( x ) ˇ u p µ , ∇ a u pT ( x ) = ϕ µ ( x ) ˇ u p µ ⊗ q a µ , (11)via FDD of ∇ in terms of the e ff ective wavenumber q a µ . In particular, ∇ FD u pT ( x ) : = h − [ u pT ( x + h ) − u pT ( x )] = ϕ µ ( x ) ˇ u p µ ⊗ q FD µ = : ˜ E FD ( x ) , ∇ BD u pT ( x ) : = h − [ u pT ( x ) − u pT ( x − h )] = ϕ µ ( x ) ˇ u p µ ⊗ q BD µ = : ˜ E BD ( x ) , ∇ CD u pT ( x ) : = h − [ u pT ( x + h ) − u pT ( x − h )] = ϕ µ ( x ) ˇ u p µ ⊗ q CD µ = : ˜ E CD ( x ) , ∇ hC u pT ( x ) : = h − [ u pT ( x + h ) − u pT ( x − h ))] = ϕ µ ( x ) ˇ u p µ ⊗ q hC µ = : ˜ E hC ( x ) , (12) In (10) and all relations to follows, repeated indices, e.g., i or µ , imply a sum over these from1 to n when this makes sense . a = FD), backward ( a = BD), central ( a = CD), and half central( a = hC), di ff erence discretization, respectively, with q FD µ : = h − ( e q µ h − , q BD µ : = h − (1 − e − q µ h ) , (13)and q CD µ : = h − sinh( q µ h ) , q hC µ : = h − sinh( q µ h ) . (14)Analogous to (11), TD of T ( x ) yields T T ( x ) = ϕ µ ( x ) ˇ T µ , div b T T ( x ) = ϕ µ ( x ) ˇ T µ q b µ , (15)via FDD of the divergence operator. Together, (11) and (15) induce the precon-ditioned form − div b C H [ ∇ a u pT ( x )] = ϕ µ ( x ) ˇ G − ( q a µ , q b µ ) ˇ u p µ (16)of the homogenous elliptic operator withˇ G − ( q a µ , q b µ ) ˇ u p µ : = − C H [ ˇ u p µ ⊗ q a µ ] q b µ , ˇ G H ( q a µ , q b µ ) = − / ( C H q a µ q b µ ) , (17)the corresponding Green function ( q a µ (cid:44) , q b µ (cid:44) q b µ is determined from a given choice of a in two di ff erent ways.The first is based on conjugacy, i.e., q b µ = q a ∗ µ (18)(e.g., Willot, 2015, Equations (19) and (24)). For example, q b µ = − q BD µ for a = FDfrom (13), or q b µ = − q a µ for a = F , CD , hC via (14). The second way is based on thechoice a = FD and the observation that E FD ( x ) = ∇ FD u ( x ) determines T ( x + h )via the mean-value theorem. The choice div BD T ( x + h ) and identity ∇ BD f ( x + h ) = ∇ hC f ( x ) (19)then result in the combination ( a , b ) = (FD , hC). For this choice, then, both thestress divergence and displacement are calculated at x . Note that this choice doesnot satisfy (18).The above results are incorporated in Algorithm TD1.9 lgorithm TD1
1. given: n , C ( x ), ¯ E , C H , ˇ G ab H µ : = ˇ G H ( q a µ , q b µ )2. initialization ι = • for i = , . . . , nC i = C ( x i ); T ( ι ) i = C i ¯ E ; • for µ = , . . . , n ˇ u p( ι ) µ =
0; ˇ T ( ι ) µ = F − µ i T ( ι ) i ; ∆ ˇ u p( ι ) µ = ˇ G ab H µ ˇ T ( ι ) µ q b µ ; • ι ++ ;3. while (cid:80) n µ = | ∆ ˇ u p( ι ) µ − ∆ ˇ u p( ι − µ | (cid:62) tol & ι (cid:54) maxit • for i = , . . . , n ˜ E ( ι ) i = F + µ i ˇ u p( ι ) µ ⊗ q a µ ; T ( ι ) i = T (0) i + C i ˜ E ( ι ) i ; • for µ = , . . . , n ˇ T ( ι ) µ = F − µ i T ( ι ) i ; ∆ ˇ u p( ι ) µ = ˇ G ab H µ ˇ T ( ι ) µ q b µ ; ˇ u p( ι ) µ + = ∆ ˇ u p( ι ) µ ; • ι ++ ;The direct Fourier case is included here by defining q F µ : = q µ ("F" for "Fourier").Rather than being strain-based like the "basic scheme" of Suquet (1997) andMichel et al. (2000, 2001), note that Algorithm TD1 is displacement-based. Asevident from (17), ˇ G − ( q a µ , q b µ ) is not invertible when either q a µ or q b µ vanish. For n even, (13) and (14) imply q a , b µ = µ = m + a , b = F , FD , BD , hC, and at µ = , m + a , b = CD. For this latter case, we follow Willot (2015, Equation(29)) and formally set ˇ G H ( q a µ , q b µ ) ≡ µ =
1. For n odd, q a , b µ = µ = m + Assume now that f ( x ) in (8) is constant on each subinterval [ x i , x i + ] of [0 , l ] withvalue f c i ≡ f ( x c i ) at the mid-point or "center" x c i = ( i − ) h of [ x i , x i + ]. In thiscase, ˆ f κ in (8) reduces toˆ f c κ : = n − / s κ ˇ f c κ , s κ : = sinc( πκ/ n ) , ˇ f c κ : = n − / (cid:88) ni = e − ı k κ x c i f c i , (20)with sinc( x ) : = sin( x ) / x the (unnormalized) sinc function. As shown by compar-ison of ˆ f c κ with ˆ f t κ from (9), di ff erences between these include (i) the dependenceof ˆ f c κ on the sinc "filter" s κ , and (ii) di ff erent discretizations of [0 , l ]. Leaving the10etails to Appendix B, (20) results in the piecewise-constant discretization (PCD) f PC ( x ) : = ϕ c µ ( x ) ˇ f c µ , ϕ c µ ( x ) : = s c µ ϕ t µ ( x ) , f c i = F c + µ i ˇ f c µ , ˇ f c µ = F c − µ i f c i , (21)of (8) via (B.1) for i , µ = , . . . , n , analogous to (10) in the TD case, with s c µ : = s µ − − m and F c ± µ i : = n − / e ± q µ x c i . In particular, s c m + = s = f c µ = m + = ˇ f c κ = = n − / (cid:80) ni = f c i = n / ¯ f . Analogous to (11) and (15), then, we have u pPC ( x ) = ϕ c µ ( x ) ˇ u pc µ , T PC ( x ) = ϕ c µ ( x ) ˇ T c µ , (22)and ∇ a u pPC ( x ) = ϕ c µ ( x ) ˇ u pc µ ⊗ q a µ , div b T cPC ( x ) = ϕ c µ ( x ) ˇ T c µ q b µ . (23)Further, − div b C H [ ∇ a u pPC ( x )] = ϕ c µ ( x ) ˇ G − ( q a µ , q b µ ) ˇ u pc µ , (24)analogous to (16). Completely analogous to the case of Algorithm TD1, then,these relations and the fact that s c µ (cid:44) Algorithm PCD1
1. given: n , C ( x ), ¯ E , C H , ˇ G ab H µ : = ˇ G H ( q a µ , q b µ )2. initialization ι = • for i = , . . . , nC c i = C ( x c i ); T c( ι ) i = C c i ¯ E ; • for µ = , . . . , n ˇ u pc( ι ) µ =
0; ˇ T c( ι ) µ = F − µ i T c( ι ) i ; ∆ ˇ u pc( ι ) µ = ˇ G ab H µ ˇ T c( ι ) µ q b µ ; • ι ++ ;3. while (cid:80) n µ = | ∆ ˇ u pc( ι ) µ − ∆ ˇ u pc( ι − µ | (cid:62) tol & ι (cid:54) maxit • for i = , . . . , n ˜ E c( ι ) i = F c + µ i ˇ u pc( ι ) µ ⊗ q a µ ; T c( ι ) i = T c(0) i + C c i ˜ E c( ι ) i ; • for µ = , . . . , n ˇ T c( ι ) µ = F c − µ i T c( ι ) i ; ∆ ˇ u pc( ι ) µ = ˇ G ab H µ ˇ T c( ι ) µ q b µ ; ˇ u pc( ι ) µ + = ∆ ˇ u pc( ι ) µ ; • ι ++ ;As evident, this algorithm is independent of s c µ = s µ − − m = sinc( π ( µ − − m ) / n )from (20) . This is in contrast to the PCD-based algorithm of Eloh et al. (2019),to which we now turn. 11 .3. Algorithm of Eloh et al. (2019) Following Brisard and Dormieux (2010), Eloh et al. (2019) recently developed analgorithm also based on PCD and (20) quite di ff erent than Algorithm PCD1. Tothis end, they formulate their algorithm based on F c ∞ f ( x ) : = (cid:80) ∞ κ = −∞ e ı k κ x ˆ f c κ from(8) for m = ∞ and (20), and truncate the result. Equivalently, as done here, onecan work with the truncated form F c p f ( x ) : = (cid:88) p − κ = − p e ı k κ x ˆ f c κ = n − / (cid:88) n − ω = e ı k ω x (cid:88) m − ν = − m e ı k ν n x s ν n + ω ˇ f c ν n + ω (25)of F c ∞ f ( x ) from the start (for p = nm even). Combining this with (the second of)ˇ f c ν n + ω = ( − ν ˇ f c ω , n − / (cid:88) n − ω = e ı k ω x c i ˇ f c ω = f ( x c i ) = f c i , (26)from (20) results in the PCD f El ( x ) : = ϕ El ω ( x ) ˇ f c ω , ϕ El ω ( x ) : = n − / (cid:88) m − ν = − m e ı k ν n + ω x ( − ν s ν n + ω , (27)(sum over repeated ω ) of (8) for ω = , . . . , n −
1, analogous to (21) ("El" standsfor "Eloh"). As in the basic scheme (e.g., Suquet, 1997), they work with (i) strain˜ E El ( x ) = ϕ El ω ( x ) ˇ˜ E c ω as the primary discretant, and (ii) the Fourier transform ˆ Γ H ( k ) ˆ T ( k ) = − ˆ G H ( k ) ˆ T ( k ) k ⊗ k = , k (cid:44) , (28)of mechanical equilbrium Γ H ∗ T : = ∇ G H ∗ div T = To be more precise, again like Suquet (1997), Eloh et al. (2019) work with the polarizationstress T − C H E instead of T directly as done in Algorithm DGO1. This is also true for the 3D caseand Algorithm DGO2 below. lgorithm DGO1
1. given: n , m , C ( x ), ¯ E , C H , ˇ Γ cH ω = =
0, ˇ Γ cH ω (cid:44) = (cid:80) m − ν = − m s ν n + ω ˆ Γ H ( k ν n + ω )2. initialization ι = • for i = , . . . , nC c i = C ( x c i ); T c( ι ) i = C c i ¯ E ; • for ω = , . . . , n − E c( ι ) ω =
0; ˇ T c( ι ) ω = F − ω i T c( ι ) i ; ∆ ˇ E c( ι ) ω = ˇ Γ cH ω ˇ T c( ι ) ω ; • ι ++ ;3. while (cid:80) n − ω = | ∆ ˇ E c( ι ) ω − ∆ ˇ E c( ι − ω | (cid:62) tol & ι (cid:54) maxit • for i = , . . . , n ˜ E c( ι ) i = ˇ F c + ω i ˇ E c( ι ) ω ; T c( ι ) i = T c(0) i + C c i ˜ E c( ι ) i ; • for ω = , . . . , n − T c( ι ) ω = ˇ F c − ω i T c( ι ) i ; ∆ ˇ E c( ι ) ω = ˇ Γ cH ω ˇ T c( ι ) ω ; ˇ E c( ι ) ω + = ∆ ˇ E c( ι ) ω ; • ι ++ ;Here, ˇ F c ± ω i : = n − / e ± ı k ω x c i and k ν n + ω = π ( ν n + ω ) / l = π ( ν + ω/ n ) / h . Again, a sumon repeated i , ω is employed. Note the dependence of the discretized Lippmann-Schwinger operator ˇ Γ cH ω on s − mn + ω , . . . , s ( m − n + ω . As already noted above, this isone di ff erence between the two PCD-based algorithms PCD1 and DGO1.
4. Computational comparisons in 1D
Employing the algorithms just discussed, the 1D BVP for the periodic MI case isnow solved numerically and compared with the analytical solution. To this end,the unit cell from the analytic case in Section 2 as based on C I / C M =
100 isemployed, and the homogeneous sti ff ness C H is determined by C H = ( C M + C I ).Since C ( x ) is not defined on the MI interface in the discontinuous case, note thatthis interface always lies between two nodes, one in the matrix, and the other inthe inclusion. This is in contrast to the smooth case, for which C ( x ) is definedeverywhere. The sti ff ness profile C ( x ) = / S ( x ) for the smooth interface is basedon ν ( x ) with (cid:15)/ l = / E =
1. Inthe context of (18), results are obtained and compared for a = F , CD , FD. Theseare obtained as well for ( a , b ) = (FD , hC). In all cases, the tolerance tol is machineprecision. 13o begin, consider the results for ˜ E ( x ) in Figure 3. Figure 3: Comparison of results for ˜ E a ( x ) = ϕ µ ( x ) ˇ u p µ ⊗ q a µ obtained from Algorithm TD1 with theanalytic result ˜ E ( x ) (black curve) from (7) across a sharp (left) and a smooth (right) MI interfaceat x / l = / n =
10 ( m =
5; blue curve), (ii) n =
50 ( m =
25; greencurve), (iii) n =
90 ( m =
45; red curve). As done in Figures 1 and 2 above, results are shown hereand in what follows for the left half of the unit cell.
As implied by these results, all choices considered for q a µ and q b µ in the contextof Algorithm TD1 yield the same solution for ˜ E a ( x ) at both sharp (Figure 3, left)and smooth (Figure 3, right) MI interfaces. Moreover, as expected from Fourierseries approximation of the analytic solution and series truncation (Figure 2, left),˜ E a ( x ) at the sharp MI interface (left) exhibits Gibbs error whose magnitude isindependent of numerical resolution. This is in contrast to ˜ E a ( x ) at the smooth MIinterface (Figure 3, right), which does converge in this fashion, again as expected.Consider next the comparison of results from Algorithms TD1 and PCD1 for thecase of a sharp MI interface shown in Figure 414 (cid:1)(cid:2) (cid:1)(cid:3) (cid:1) (cid:1)(cid:2) (cid:1) / (cid:2)(cid:3) (cid:4) (cid:4) (cid:2) (cid:1) (cid:1)(cid:2) (cid:1)(cid:3) - (cid:2)(cid:1)(cid:2) (cid:1) / (cid:2)(cid:4) (cid:4) Figure 4: Left: comparison of results for u pT ( x ) = ϕ µ ( x ) ˇ u p µ from Algorithm TD1 (red points, curve)and for u pPC ( x ) = ϕ c µ ( x ) ˇ u pc µ from Algorithm PCD1 (blue points, curve) with analytic solution (blackcurve) across a sharp MI interface at x / l = .
25. Right: comparison of results for ˜ E T ( x ) : = ∇ F u pT ( x ) = ϕ µ ( x ) ˇ u p µ ⊗ q F µ from Algorithm TD1 (red points, curve) and for ˜ E PC ( x ) : = ∇ F u pPC ( x ) = ϕ c µ ( x ) ˇ u pc µ ⊗ q F µ from Algorithm PCD1 (blue points, curve) with analytic solution (black curve).Results are based on discretizations of n =
30 ( m =
15; red points) and n =
32 ( m =
16; bluepoints).
As seen in particular in the strain results on the right, the dependence of ϕ c µ ( x ) in(21) on s c µ dampens Gibbs oscillations / error in ˜ E PC ( x ) in comparison to ˜ E T ( x ).As expected, this has no e ff ect on the numerical (i.e., collocation) solution at thenodes, in contrast to the analogous results from Algorithm DGO1 in Figure 5. Figure 5: Left: Results for nodal ˜ E c i from Algorithm DGO1 across a sharp MI interface at x / l = / n =
10 ( m =
5; blue squares), (ii) n =
50 ( m =
25; red circles), (iii) n =
90 ( m =
45; green triangles), compared with the analytic result ˜ E ( x ) (black curve) from (7) .Righ: Comparison of results for nodal ˜ E i from Algorithm TD1 (blue asterisks) and those for ˜ E c i from Algorithm DGO1 (red circles) for the discretization n =
50 ( m = As shown on the left, the nodal results from Algorithm DGO1 deviate slightlyfrom those of the anayltic solution in the (relatively soft) matrix, with maximumdeviation in the matrix next to the MI interface. Note that this deviation decreases15ith increasing numerical resolution. This is in contrast to the analogous solutionsfrom Algorithm TD1 (and so Algorithm PCD1) as shown by the comparison onthe right. One di ff erence between the PCD-based Algorithms PCD1 and DGO1that is likely playing a role here is the dependence of the discretized Lippmann-Schwinger operator ˇ Γ cH ω on s ν n + ω in Algorithm DGO1. These and other di ff erencesin 1D become even more evident in the context of tensor-product-based general-ization of the above 1D algorithms to 3D, to which we now turn.
5. Algorithms for strong mechanical equilibrium in 3D
Let u , H = ∇ u , and E = sym H represent the 3D displacement, distortion andstrain, respectively. In this case, let f = ¯ f + ˜ f now represent the 3D generalizationof (1) with respect to the unit cell U . Likewise, u = u h + u p , u h ( x ) = c + ¯ Hx , u p ( x ) = (cid:90) ˜ H ( x ) d x , ∇ u p = ˜ H , (29)is the generalization of (2) to 3D with ¯ u p = , and so ¯ u h = c + ¯ H ¯ x = ¯ u . Asin the 1D case, c , ¯ H , and so u h , are known, and attention is focused on u p . Onthis basis, quasi-static mechanical equilibrium and isotropic, linear elastic ma-terial behavior continue to apply. Let T = C E be the linear elastic stress, and C ( x ) = λ ( x ) I ⊗ I + µ ( x ) ( I (cid:3) I + I (cid:52) I ) the isotropic sti ff ness. For the case of dis-continuous C ( x ), the following hold at the MI interface (with unit normal n ): (i)no cracking u I = u M , (ii) kinematic compatibility H I = H M + h ⊗ n , and (iii)mechanical equilibrium T I n = T M n . Both λ ( x ) and µ ( x ) are determined by directtensor-product generalization of analogous 1D relations like (6) based on θ in thediscontinuous, and on ν in the continuous, case. Let r = , ,
3. As in the 1D case, consider a uniform discretization of [0 , l r ] basedon n r + x i r = i r h r with i r = , . . . , n r and spacing h r such that l r = n r h r ( n r subintervals). Direct tensor-product-based generalization of (10) to 3D yields f T ( x ) = ϕ µ ( x ) ϕ µ ( x ) ϕ µ ( x ) ˇ f µ µ µ = : ϕ µ ( x ) (cid:63) ˇ f µ (30)16sum on repeated indices) in terms of Rayleigh product notation with x : = ( x , x , x )and µ r = , . . . , n r , r = , ,
3. Here, f i : = f i i i = F + µ i F + µ i F + µ i ˇ f µ µ µ = : F + µ i (cid:63) ˇ f µ , ˇ f µ : = ˇ f µ µ µ = F − µ i F − µ i F − µ i f i i i = : F − µ i (cid:63) f i , (31)(sum on repeated indices) via (10) , . Given these, generalization of AlgorithmTD1 to 3D is based in particular on those ∇ a u pT ( x ) = ϕ µ ( x ) (cid:63) ˇ u p µ ⊗ q a µ , div b T T ( x ) = ϕ µ ( x ) (cid:63) ˇ T µ q b µ , (32)of (11) and (15) , respectively. Here, ˇ u p µ : = ( ˇ u p1 µ , ˇ u p2 µ , ˇ u p3 µ ), ˇ T µ is the 3 × q a µ : = ( q a µ , q a µ , q a µ ) with q µ r : = ı k µ r − − m r .As in the 1D case (17), 3D quasi-static momentum balance is formulated algo-rithmically in preconditioned form based on the Green function ˇ G H ( q a µ , q b µ ) of thecorresponding operator − div b C H ∇ ah u ( x ) = − ϕ µ ( x ) (cid:63) C H [ ˇ u µ ⊗ q a µ ] q b µ = ϕ µ ( x ) (cid:63) ˇ G − ( q a µ , q b µ ) ˇ u µ . (33)Given these relations, direct componentwise generalization of Algorithm TD1yields Algorithm TD2. Algorithm TD2
1. given: n , . . . , n , C ( x ), ¯ E , C H , ˇ G ab H µ = ˆ G H ( q a µ , q b µ )2. initialization ι = • for i = , . . . , n , . . . , i = , . . . , n C i = C ( x i ); T ( ι ) i = C i ¯ E ; • for µ = , . . . , n , . . . , µ = , . . . , n ˇ u p( ι ) µ = ; ˇ T ( ι ) µ = F − µ i (cid:63) T ( ι ) i ; ∆ ˇ u p( ι ) µ = ˇ G ab H µ ˇ T ( ι ) µ q b µ ; • ι ++ ;3. while (cid:80) n µ = · · · (cid:80) n µ = | ∆ ˇ u p( ι ) µ − ∆ ˇ u p( ι − µ | (cid:62) tol & ι (cid:54) maxit • for i = , . . . , n , . . . , i = , . . . , n ˜ E ( ι ) i = F + µ i (cid:63) sym( ˇ u p( ι ) µ ⊗ q a ( ι ) µ ); T ( ι ) i = T (0) i + C i ˜ E a ( ι ) i ; • for µ = , . . . , n , . . . , µ = , . . . , n ˇ T ( ι ) µ = F − µ i (cid:63) T ( ι ) i ; ∆ ˇ u p( ι ) µ = ˇ G ab H µ ˇ T ( ι ) µ q b µ ; ˇ u p( ι ) µ + = ∆ ˇ u p( ι ) µ ; • ι ++ ; 17or the current isotropic case, C H = λ H I ⊗ I + µ H ( I (cid:3) I + I (cid:52) I ), and soˇ G − ( q a µ , q b µ ) : = − µ H ( q a µ · q b µ ) I − µ H q a µ ⊗ q b µ − λ H q b µ ⊗ q a µ . (34)Since det ˇ G − ( q a µ , q b µ ) = − ( λ H + µ H +
3) ( q a µ · q b µ ), note that ˆ G − is invertible for q a µ · q b µ (cid:44)
0. In the context of conjugacy (e.g., (45) below in 3D), note that q b µ ⊗ q a µ = q a µ ⊗ q b µ hold when q a r ∗ µ r = − q a r µ r (e.g., a r = ACD , AhC), and so ˇ G − ( q a µ , q b µ ), issymmetric. For other cases (e.g., AFB / R; see below), this is not true in general.Restriction to λ H = µ H (e.g., Willot, 2015), however, does result in symmetricˇ G − ( q a µ , q b µ ) for all q a µ and q b µ . As a first step toward 3D generalization of the 1D FDDs, consider first the 2Dcase. To this end, it is useful to work with the di ff erence operators δ FD h i f ( . . . , x i , . . . ) : = h − i [ f ( . . . , x i + h i , . . . ) − f ( . . . , x i , . . . )] ,δ BD h i f ( . . . , x i , . . . ) : = h − i [ f ( . . . , x i , . . . ) − f ( . . . , x i − h i , . . . )] ,δ CD h i f ( . . . , x i , . . . ) : = h − i [ f ( . . . , x i + h i , . . . ) − f ( . . . , x i − h i , . . . )] ,δ hC h i f ( . . . , x i , . . . ) : = h − i [ f ( . . . , x i + h i , . . . ) − f ( . . . , x i − h i , . . . )] . (35)Given these, consider the 2D grid "cell" with corners at ( x , x ), ( x , x + h ),( x + h , x + h ), and ( x + h , x ). With respect to this cell, average forwarddi ff erencing (AFD) is defined as ∇ AFD i f ( x , x ) : = [ δ FD h f ( x , x ) + δ FD h f ( x , x + h )] , ∇ AFD i f ( x , x ) : = [ δ FD h f ( x , x ) + δ FD h f ( x + h , x )] . (36)Analogously, average backward di ff erencing (ABD) is defined as ∇ ABD i f ( x , x ) : = [ δ BD h f ( x , x ) + δ BD h f ( x , x − h )] , ∇ ABD i f ( x , x ) : = [ δ BD h f ( x , x ) + δ BD h f ( x − h , x )] , (37)with respect to the grid cell with corners at ( x , x ), ( x , x − h ), ( x − h , x − h ),and ( x − h , x ). In the same fashion, we have ∇ ACD i f ( x , x ) : = [ δ CD h f ( x , x − h ) + δ CD h f ( x , x + h )] , ∇ ACD i f ( x , x ) : = [ δ CD h f ( x − h , x ) + δ CD h f ( x + h , x )] , (38)18or average central di ff erencing (ACD), and ∇ AhC i f ( x , x ) : = [ δ hC h f ( x , x − h ) + δ hC h f ( x , x + h )] , ∇ AhC i f ( x , x ) : = [ δ hC h f ( x − h , x ) + δ hC h f ( x + h , x )] , (39)for average half-central di ff erencing (AhC). From (37) and (39) follow the directgeneralization ∇ ABD i , f ( x + h , x + h ) = ∇ AhC i , f ( x , x ) (40)of (19). To obtain the modal forms of these, the notation f ( x + c h , . . . , x d + c d h d ) = e c q µ h · · · e c d q µ d h d f ( x , . . . , x d ) (41)based on f ( x , . . . , x d ) = ˇ f µ ··· µ d e q µ x · · · e q µ d x d is useful. In terms of this,ˇ ∇ AFD i r : = h − r ( e q µ r h r −
1) ( e q µ s h s + = q FD µ r ( e q µ s h s + = : q AFD µ r , ˇ ∇ ABD i r : = h − r (1 − e − q µ r h r ) ( e − q µ s h s + = q BD µ r ( e − q µ s h s + = : q ABD µ r , ˇ ∇ ACD i r : = h − r ( e q µ r h r − e − q µ r h r ) ( e q µ s h s + e − q µ s h s ) = q CD µ r cosh( q µ s h s ) = : q ACD µ r , ˇ ∇ AhC i r : = h − r ( e q µ r h r − e − q µ r h r ) ( e q µ s h s + e − q µ s h s ) = q hC µ r cosh( q µ s h s ) = : q AhC µ r , (42)via (13)-(14) with s = , r = ,
2. Direct generalization of (42) to 3D yieldsˇ ∇ AFD i r : = q FD µ r ( e q µ s h s +
1) ( e q µ t h t + = : q AFD µ r , ˇ ∇ ABD i r : = q BD µ r ( e − q µ s h s +
1) ( e − q µ t h t + = : q ABD µ r , ˇ ∇ ACD i r : = q CD µ r cosh( q µ s h s ) cosh( q µ t h t ) = : q ACD µ r , ˇ ∇ AhC i r : = q hC µ r cosh( q µ s h s ) cosh( q µ t h t ) = : q AhC µ r . (43)with ( s , t ) = (2 , , (3 , , (1 ,
2) for r = , ,
3. Yet another e ff ective wavenumber q R µ r : = q W µ r ( e q µ s h s + e q µ t h t + , q W µ : = h − tanh (cid:32) q µ h (cid:33) ( e q µ h + , (44)(in the current notation) was obtained by Willot (2015, Equation (36)) in the con-text of his "rotated scheme" (R). As it turns out, q W µ = q FD µ for µ = , . . . , n (recall that q µ h = ı k µ − − m h = πı ( µ − − m ) / n ). On the other hand, for µ = q W µ (cid:44) q FD µ since q W µ = is indeterminate (note q µ = h = − πı for m = n / µ → q W µ = q FD µ = does hold. Consequently, AFB- and R-based algorithms will betreated as as equivalent in what follows.Analogous to the 1D case, these results can be employed to formulate two typesof FDDs. Again following Willot (2015), the first type is based on the directgeneralization q b r µ r = q a r ∗ µ r (45)of (18) to 3D. In particular, note that q b r µ r = − q ABD µ r for a r = AFD and q b r µ r = − q a r µ r for a r = ACD , AhC. The second type generalizes the choice ( a , b ) = (FD , hC)in 1D to ( a r , b r ) = (AFD , AhC) in 3D, and is referred to as the "average forward-backward / rotated" (AFB / R) FDD in what follows.
Direct tensor-product-based generalization of (21) to 3D yields f PC ( x ) = ϕ c µ ( x ) ϕ c µ ( x ) ϕ c µ ( x ) ˇ f c µ µ µ = : ϕ c µ ( x ) (cid:63) ˇ f c µ (46)again in terms of Rayleigh product notation with f c i : = f c i i i = F c + µ i F c + µ i F c + µ i ˇ f c µ µ µ = : F c + µ i (cid:63) ˇ f c µ , ˇ f c µ : = ˇ f c µ µ µ = F c − µ i F c − µ i F c − µ i f c i i i = : F c − µ i (cid:63) f c i , (47)(again sum on repeated indices) via (21) , and analogous to (31). On this basis, Algorithm PCD2
1. given: n , . . . , n , C ( x ), ¯ E , C H , ˇ G ab H µ = ˆ G H ( q a µ , q b µ )2. initialization ι = • for i = , . . . , n , . . . , i = , . . . , n C c i = C ( x c i ); T c( ι ) i = C c i ¯ E ; • for µ = , . . . , n , . . . , µ = , . . . , n ˇ u pc( ι ) µ = ; ˇ T c( ι ) µ = F c − µ i (cid:63) T c( ι ) i ; ∆ ˇ u pc( ι ) µ = ˇ G ab H µ ˇ T c( ι ) µ q b µ ; • ι ++ ;3. while (cid:80) n µ = · · · (cid:80) n µ = | ∆ ˇ u pc( ι ) µ − ∆ ˇ u pc( ι − µ | (cid:62) tol & ι (cid:54) maxit • for i = , . . . , n , . . . , i = , . . . , n ˜ E c( ι ) i = F c + µ i (cid:63) sym( ˇ u pc( ι ) µ ⊗ q a µ ); T c( ι ) i = T c(0) i + C c i ˜ E c( ι ) i ; • for µ = , . . . , n , . . . , µ = , . . . , n ˇ T c( ι ) µ = F c − µ i (cid:63) T c( ι ) i ; ∆ ˇ u pc( ι ) µ = ˇ G ab H µ ˇ T c( ι ) µ q b µ ; ˇ u pc( ι ) µ + = ∆ ˇ u pc( ι ) µ ; • ι ++ ; 20 lgorithm DGO2
1. given: n , . . . , n , C ( x ), ¯ E , C H , ˇ Γ cH = ,ˇ Γ cH ω (cid:44) = (cid:80) m − ν = − m s ν n + ω · · · (cid:80) m − ν = − m s ν n + ω ˆ Γ H ( k ν n + ω , k ν n + ω , k ν n + ω )2. initialization ι = • for i = , . . . , n , . . . , i = , . . . , n C c i = C ( x c i ); T c( ι ) i = C c i ¯ E ; • for ω = , . . . , n − . . . , ω = , . . . , n − E c( ι ) ω = ; ˇ T c( ι ) ω = F c − ω i (cid:63) T c( ι ) i ; ∆ ˇ E c( ι ) ω = ˇ Γ cH ω ˇ T c( ι ) ω ; • ι ++ ;3. while (cid:80) n − ω = · · · (cid:80) n − ω = | ∆ ˇ E c( ι ) ω − ∆ ˇ E c( ι − ω | (cid:62) tol & ι (cid:54) maxit • for i = , . . . , n , . . . , i = , . . . , n ˜ E c( ι ) i = F c + ω i (cid:63) ˇ E c( ι ) ω ; T c( ι ) i = T c(0) i + C c i ˜ E c( ι ) i ; • for ω = , . . . , n − . . . , ω = , . . . , n − T c( ι ) ω = F c − ω i (cid:63) T c( ι ) i ; ∆ ˇ E c( ι ) ω = ˇ Γ cH ω ˇ T c( ι ) ω ; ˇ E c( ι ) ω + = ∆ ˇ E c( ι ) ω ; • ι ++ ;are obtained in 3D via direct componentwise generalization of Algorithm PCD1and Algorithm DGO1, respectively. In particular, both algorithms are based onthe Cartesian components of the Fourier transformˆ G H ( k ) = µ H | k | (cid:34) I − + λ H /µ H + (1 + λ H /µ H ) k | k | ⊗ k | k | (cid:35) , k (cid:44) , (48)of the isotropic Green tensor. In addition, Algorithm DGO2 utilizes the Cartesiancomponents of the Fourier transform ˆ Γ H of the 3D Lippmann-Schwinger operator Γ H , defined by ˆ Γ H ( k ) A : = − sym[ ˆ G H ( k ) A ( k ⊗ k )] for all A .
6. Computational comparisons in 3D
In this section, the algorithms formulated in the last section for solution of thestrong BVP are compared in the 3D matrix-inclusion (MI) case. As in the 1D caseabove, both sharp and smooth MI interfaces are considered. More specifically,in the context of TD and Algorithm TD2, results are compared for the choices a r = F , CD , ACD in the context of (45), as well as for ( a r , b r ) = (AFD , AhC)(i.e., AFB / R). As well, in the context of PCD, Algorithm PCD2 for ( a r , b r ) = (AFD , AhC) is compared with Algorithm DGO2. For reference, the results from21hese in the context of the strong BVP are compared with analogous results fromthe numerical solution of the corresponding weak BVP via standard finite ele-ment (SFE) discretization. In this latter case, note that the sharp MI interface isdiscretized by element boundaries.For comparability with the literature, the algorithmic comparisons just discussedare carried out in the sequel employing the MI benchmark example from Willot(2015). In particular, this is based on the choices λ M = µ M = . λ H = λ I ,and µ H = µ I . The 3D computational domain Ω of Willot (2015) is discretizeduniformly based on L nodes and unit grid spacing h =
1. Since his L then cor-responds to n r here, n = n = n = n and h = h = h = h = T /µ M in whatfollows for di ff erent phase contrasts χ = λ I /λ M = µ I /µ M . (49)In terms of χ , note that λ ( x ) /µ M = ( λ M /µ M ) f ( x , χ ) and µ ( x ) /µ M = f ( x , χ ) hold,with f ( x , χ ) : = + ν (cid:15) ( x ) ν (cid:15) ( x ) ν (cid:15) ( x ) ( χ − E , with ¯ E xy =
1, and all other components zero.
Results for the stress field at a sharp MI interface with a phase contrast of χ = igure 6: Comparison of discrete results for T xy ( x i , x i , z I ) /µ M for a sharp MI interface with χ = ff erent choices of a r , b r in Algorithm TD2 and at di ff erent numerical resolutions n . Inparticular, n =
21 (left), n =
41 (middle left), n =
81 (middle right), and n =
161 (right). The samesolutions are obtained for the even choices n =
22 (left), n =
42 (middle left), n =
82 (middleright), and n =
162 (right). The discrete results shown are from the grid points in the ( x , x ) planeat x i = z I inside the inclusion and adjacent to the MI interface. The analogous results for T xy at a smooth MI interface are shown in Figure 7.23 igure 7: Same as Figure 6 for a smooth MI interface. Here, T xy ( x i , x i , z MI ) /µ M is shown fromthe grid points in the ( x , x ) plane on the MI interface at x i = z MI . As in the 1D case and Figure3, the smooth transition in sti ff ness between matrix and inclusion is based on (cid:15)/ l = /
100 in eachdirection.
As seen in Figure 6 (top row), for the discontinuous MI interface case in thecontext of (45), the choice a r = F yields qualitatively incorrect results at theresolutions considered. Indeed, as shown by Willot (2015, Figure 4), who workedwith much finer discretizations of L = , , ff erence choice a r = CD (Figure 6, second row from the top) yieldsa stress field which is qualitatively correct and nearly converges at the highestresolution investigated here. Among the choices based on conjugacy (45), a r = ACD (Figure 6, middle row) converges almost as quickly as the non-conjugate24FB / R choice ( a r , b r ) = (AFD , AhC) (Figure 6, second row from the bottom),which is closest to the SFE-based reference results (Figure 6, bottom row).Analogous to the 1D case (Figure 3), and as expected, a quite di ff erent pictureemerges for the convergence behavior and stress field at a smooth MI interfacein 3D. Indeed, as for the 1D results in Figure 3 (right), the convergence behaviordisplayed in Figure 7 in the 3D case for all choices of ( a r , b r ) is a ff ected predomi-nantly by numerical resolution. In particular, as in the discontinuous MI interfacecase, among the choices based on conjugacy (45), a r = ACD converges almost asquickly as the non-conjugate AFB / R choice ( a r , b r ) = (AFD , AhC), which againis closest to the SFE-based reference case (bottom row).As evident in the results from Figures 6 and 7, the stress field at the MI corners ismost sensitive to numerical resolution. To look at this in more detail, consider theresults in Figure 8.
Figure 8: Discrete AFB / R (left, triangles) and SFE (right, circles) results for T xy ( x i , y I , z I ) across aMI corner for the sharp (left) MI interface. Results shown are at grid points in the x direction forfixed x i = y I , x i = z I just inside the inclusion and for the numerical resolutions of n r =
42 (blue), n r =
82 (green), n r =
162 (red).
As expected, in the sharp interface case, the solution does not converge with in-creasing resolution, i.e., is mesh-dependent, in contrast to the smooth interfacecase. The convergence behavior of the ACD ( a r = ACD) algorithm based on con-jugacy (45), as well the AFB / R (( a r , b r ) = (AFD , AhC)), are compared for boththe sharp and smooth MI interfaces are compared in Figure 9.25 -3 -2 -1 Figure 9: Convergence behavior of selected algorithms as a function of phase contrast at sharp(solid curves) and smooth (dashed curves) MI interfaces for n r = , / R.Green: a r = ACD. No convergence results for a r = F are shown because this case does notconverge in less than 10 iterations at the highest contrasts. As shown, for phase contrasts up to a factor of 10 (i.e., χ = − , ), little di ff er-ence in convergence rate is apparent. For higher contrasts, however, the e ff ect ofoscillations due to Gibbs and aliasing on the convergence rate for the sharp inter-face cases becomes apparent. Indeed, as expected from the stress results in Figure6, the e ff ective low-pass filtering e ff ect of finite-di ff erence discretization of di ff er-ential operators for a r = ACD and AFB / R reduces the e ff ect of oscillations dueto Gibbs and aliasing on convergence rate. In particular, for high phase contrasts,the convergence rate for AFB / R at n r = ,
162 is about 10 times faster than for a r = F, and about 5 times faster than for a r = ACD, for the current benchmarkcase.
Lastly, results from Algorithm PCD2 for ( a r , b r ) = (AFD , AhC) (i.e., AFB / R) arecompared with corresponding ones from Algorithm DGO2 in Figure 10.26 igure 10: PCD-based results from AFB / R (above), the DGO (middle) for T c xy ( x c i , x c i , z cI ) /µ M andSFE (below), in the sharp MI interface case for di ff erent phase contrasts χ at a resolution of n r =
40. The discrete results shown are from the grid points in the inclusion adjacent ( x c i = z cI ) tothe MI interface. See text for discussion. Except for the smaller phase contrasts and consequently lower stress levels, thePCD-based results for AFB / R in Figure 10 (upper row) agree with the TD-basedresults in Figure 6 (middle left row). For the chosen nodal resolution of n r = m r = n r =
20 in AFB / R. For comparability, then, m r =
20 waschosen to obtain the DGO-based results. Except for the much higher stress at theMI corners obtained with the DGO, the results from both algorithms agree (evenquantitatively) well for χ =
10. With increasing χ , however, the solution based onthe DGO is much sti ff er than than based on AFB / R, indicating lack of convergenceat this resolution. Consequently, the DGO converges more slowly than AFB / R inthe context of PCD.As mentioned in the introduction, in contrast to the current case of cubic inclu-sions, Eloh et al. (2019, §§4.3-4.4) consider spherical inclusions for which ana-lytic solutions exist (e.g., Mura, 1987). In addition, following Anglin et al. (2014),they work with phase contrasts χ of 0 . κ : = ( λ + µ/ λ I = χλ M and µ I = χµ M from (49), note that κ I = χ κ M also27olds. Consequently, they worked with smaller phase contrasts than χ = κ H = ( κ M + κ I ), ν H = ν M = ν I , and λ M /µ M = .
2, the latter in contrast to λ M /µ M =
7. Summary and discussion
In the current work, a number of algorithms have been developed and comparedfor the numerical solution of periodic boundary-value problems (BVPs) for thequasi-static mechanical behavior of heterogeneous linear elastic materials basedon Fourier series discretization. In particular, two such discretizations are con-sidered and compared here. The first is based on trapezoidal approximation ofthe integral (8) over the unit cell determining the Fourier modes ˆ f κ in the trun-cated Fourier series (8) , resulting in the trapezoidal discretization (10) of thisseries. The second is based on approximation of the integrand in (8) as piecewise-constant, yielding the piecewise-constant discretization (21) of (8) . A compari-son of these two implies that a basic di ff erence between them lies in the depen-dence of (21) on s κ = sinc( πκ/ n ). As shown by Algorithms TD1 and PCD1,despite this di ff erence, equivalent algorithms based on TD and PCD can neverthe-less be formulated via the fact that s κ is non-zero (and so can be eliminated fromthe algorithm). This is in contrast to the PCD (25) of (8) employed by Eloh et al.(2019) and resulting in Algorithm DGO1. Indeed, the DGO ˇ Γ cH ω from (B.6) in thelatter algorithm depends explicitly on s ν n + ω .All three algorithms TD1, PCD1 and DGO1 are based on Green function precon-ditioning (GFP). In addition, the first two exploit finite di ff erence discretization(FDD) of di ff erential operators, in particular of ∇ and div. In the context of TD(10) and PCD (21) of fields, di ff erent FDDs for these two operators are obtainedin modal form (e.g., (11)-(12) for ∇ ) in terms of e ff ective wavenumbers q a µ (for ∇ a )and q b µ (for div b ). Specific FDDs consider include forward ( a , b = FD), backward( a , b = BD), central ( a , b = CD), and half central ( a , b = hC), di ff erence, respec-tively. Following in particular Willot (2015), given a , one choice for b is basedon conjugacy (18). A second one is based on choosing b in such a way that thediscretized stress divergence is determined in the algorithm at the (displacement)nodes. For example, for the choice a = FD, this results in b = hC, and so the(non-conjugate) combination ( a , b ) = (FD , hC).28omputational comparison of these algorithms is carried out in Section 4 for thematrix-inclusion benchmark case in 1D for material heterogeneity in the form ofboth discontinuous and smooth compliance / sti ff ness distributions. In the formercase, as expected from Fourier series approximation of the analytic solution andseries truncation (Figure 2, left), ˜ E a ( x ) at the sharp MI interface in Figure 3 (left)exhibits Gibbs error independent of, and no convergence with increasing, numer-ical resolution for all choices of ∇ a ( q a µ ) and div b ( q b µ ) in the context of AlgorithmsTD1 and PCD1. On the other hand, ˜ E a ( x ) does converge with increasing numer-ical resolution at the smooth MI interface (Figure 3, right). In addition, for thesharp interface case, comparison of strain results from Algorithms TD1, PCD1and DGO1 in Figure 5 show that the nodal results from PCD1 deviate slightlyfrom those of the anayltic solution and the first two algorithms in the (relativelysoft) matrix, with maximum deviation in the matrix next to the MI interface.Multidimensional versions of the 1D algorithms are formulated in Section 5 viadirect componentwise- and tensor-product-based generalization. In this fashion,3D versions TD2, PCD2 and DGO2, respectively, of the 1D algorithms TD1,PCD1 and DGO1, respectively, follow directly. As in the 1D case, all three arebased on GFP, and the first two on FDDs of ∇ and div, now in 3D. In the contextof either TD (30) or PCD (46) of 3D fields, the latter are formulated in 3D viacomponentwise application of the 1D approaches, resulting in the FDDs ∇ a anddiv b with a = ( a , a , a ) and b = ( b , b , b ). As in 1D, choices for b r given a r are based on conjugacy (45) and the non-conjugate stress divergence criteria, thelatter yielding in particular b r = AhC for a r = AFD and so the "average forward-backward / rotated" FDD, i.e., AFB / R. Computational comparsions of these arepresented in Section 6. Generally speaking, an increase in phase contrast resultsin an increase in the condition number of the algorithmic equation system beingsolved, resulting in slower convergence. As such, preconditioning is clearly es-sential to improved convergence. As shown for example in Willot (2015) and thecurrent work (e.g., Figure 9), the combination of GFP and FDD results in signifi-cant further improvement in convergence rate and behavior over approaches basedon GFP alone such as the DGO of Eloh et al. (2019).In the current work, the focus has been on material inhomogeneity with respectto elastic sti ff ness. Generalization of the current algorithms to the case of suchheterogeneity with respect to both elastic sti ff ness and residual strain (e.g., dueto lattice mismatch between phases) is straightforward. In the 3D case, the corre-sponding generalization T = C [ E − E ∗ ] of the stress-strain relation (with C ( x ) and E ∗ ( x ) known) induces changes in the algorithms. For example, this results in the29hange of T (0) i = C i ¯ E to T (0) i = C i [ ¯ E − E ∗ i ] in the TD-based algorithm TD2. Like-wise, T c(0) i = C c i ¯ E generalizes to T c(0) i = C c i [ ¯ E − E c ∗ i ] in the PCD-based algorithmsPCD2 and DGO2.In the context of strong mechanical equilbrium, a displacement-based approachrelated to the current one has quite recently been developed by Lucarini and Se-gurado (2019). They refer to this approach as displacement-based FFT (DBFFT).As the name implies, they also work with the displacement as the primary disc-retant and the split (29). In particular, their algorithm is based on quasi-staticmechanical equilibrium in Fourier space in the form ˆ A ( q ) ˆ u p = ˆ C [ ¯ E ] q for q (cid:44) in terms of the generalized acoustic tensor ˆ A ( q ) ˆ v : = −F { C [ F − { ˆ v ⊗ q } ] } q . Asnoted by them, Fourier discretization and real-function-based reduction yields adiscrete system for displacement based on a full-rank associated Hermitian ma-trix. In turn, this facilitates use of preconditioners. In the current context, theirchoice corresponds in particular to those C H = ¯ C and ( a r , b r ) = (F , F) in (33).In terms of Fourier transforms, this results in ˆ G − ( q ) a = − ¯ C [ a ⊗ q ] q for q (cid:44) .Algorithmic extension of this to q = then results in a preconditioning operatorˆ M ( q ) which approximately inverts ˆ A ( q ) for use in iterative numerical solution.Traditional (i.e., linear elastic) micromechanics based on a (linear) stress-strainrelation T ( E ) treats the (linear) strain E (and not displacement) as the primaryunknown. In the computational context, this translates into the treatment of E ,or more recently in the geometrically non-linear case, the deformation gradient F , as the primary discretant. In this case, compatibility needs to be enforced viaadditional algorithmic constraints. Alternatively, as discussed in the current work,one can work directly with the displacement or deformation field as the primarydiscretant. This is true in both the current case of strong mechanical equilibriumas well as in the case of weak mechanical equilibrium, as shown by the recentwork of de Geus et al. (2017). To discuss the latter briefly here, consider weakequilibrium (cid:82) T ( x ) · E v ( x ) dv ( x ) = (cid:82) ˆ T ( k ) · ˆ E ∗ v ( k ) dv ( k ) (i.e., via the Rayleigh-Plancherel or power theorem: e.g., Bracewell, 2000, pp. 119-120) for any "virtual"or "test" strain field E v . Given H compatible, ˆ H ∗ v = ˆ u v ⊗ q ∗ holds, again with q = ı k . Then ˆ u v = ˆ H ∗ v q / | q | for q (cid:44) , inducing in turn ˆ T · ˆ E ∗ v = P H ˆ T · ˆ H ∗ v with P H : = I (cid:3) ( q ⊗ q ∗ ) / | q | . Consequently, the projection P H enforces compatibilitywhen H is the primary discretant and unknown. Alternatively, as done in thecurrent work, one can simply work directly with the displacement field or thedeformation field as the primary discretant. Acknowledgements . Financial support of Subproject M03 of the Transregional30ollaborative Research Center SFB / TRR 136 by the German Science Foundation(DFG) is gratefully acknowledged.
References
Anglin, B. S., Lebensohn, R. A., Rollett, A. D., 2014. Validation of a numericalmethod based on fast Fourier transforms for heterogeneous thermoelastic mate-rials by comparison with analytical solutions. Computational Materials Science87, 209–217.Bracewell, R. N., 2000. The Fourier Transform and Its Application, 3rd Edition.McGraw-Hill.Brisard, S., Dormieux, L., 2010. FFT-based methods for the mechanics of com-posites: a general variational framework. Computational Materials Science 49,663–671.Brown, C. M., Dreyer, W., Müller, W., 2002. Discrete fourier transforms and theirapplication to stress-strain problems in composite mechanics: a convergencestudy. Proceedings of the Royal Society of London A 458, 1967–1987.Canuto, C., Hussaini, M. Y., Quateroni, A., Zang, T. A., 1988. Spectral Methodsin Fluid Dynamics. Springer Series in Computational Physics. Springer-Verlag,Berlin.Chen, L.-Q., 2002. Phase-field model for microstructure evolution. Annual Re-view of Material Research 32, 113–140.de Geus, T. W. J., Vondˇrejc, J., Zeman, J., Peerlings, R. H. J., Geers, M. G. D.,2017. Finite strain FFT-based non-linear solvers made simple. Computer Meth-ods in Applied Mechanics and Engineering 318, 412–430.Djaka, K. S., Villani, A., Taupin, V., Capolungo, L., Berbenni, S., 2017. FieldDislocation Mechanics for heterogeneous elastic materials: a numerical spec-tral approach. Computer Methods in Applied Mechanics and Engineering 315,921–942.Dreyer, W., Müller, W., 2000. A study of the coarsening of tin / lead solders. Inter-national Journal of Solids and Structures 37, 3841–3871.31isenlohr, P., Diehl, M., Lebensohn, R. A., Roters, F., 2013. A spectral methodsolution to crystal elasto viscoplasticity at finite strains. International Journal ofPlasticity 46, 37–53.Eloh, K. S., Jacques, A., Berbenni, S., 2019. Development of a new consistent dis-crete Green operator for FFT- based methods to solve heterogeneous problemswith eigenstrains. International Journal of Plasticity 116, 1–23.Eyre, D. J., Milton, G. W., 1999. A fast numerical scheme for computing theresponse of composites using grid refinement. The European Physical Journal6, 41–47.Gottlieb, D., Orszag, S. A., 1977. Numerical Analysis of Spectral Methods: The-ory and Application. SIAM, Philadelphia.Kabel, M., Böhlke, T., Schneider, M., 2014. E ffi cient fixed point and Newton-Krylov solvers for FFT-based homogenization of elasticity at large deforma-tions. Computational Mechanics 54, 1497–1514.Khachaturyan, A. G., 1983. Theory of Structural Transformations in Solids. Wi-ley, New York.Kopriva, D. A., 2009. Implementing Spectral Methods for Partial Di ff erentialEquations. Springer.Lebensohn, R. A., 2001. N-site modeling of a 3d viscoplastic polycrystal usingFast Fourier Transform. Acta Materialia 49, 2723–2737.Lebensohn, R. A., Rollett, A. D., 2020. Spectral methods for full-field microme-chanical modelling of polycrystalline materials. Computational Materials Sci-ence 173, 109336.Liu, S. H., 2011. Numerical Analysis of Partial Di ff erential Equations. Wiley.Lucarini, S., Segurado, J., 2019. DBFFT: a displacement based FFT approach fornon-linear homogenization of the mechanical behavior. International Journal ofEngineering Science 144, 103131.Michel, J. C., Moulinec, H., Suquet, P., 2000. A computational method based onaugmented Lagrangians and Fast Fourier Transforms for composites with highcontrast. Computational Modelling in Engineering Science 1, 79–88.32ichel, J. C., Moulinec, H., Suquet, P., 2001. A computational scheme for linearand non-linear composites with arbitrary phase contrast. International Journalof Numerical Methods in Engineering 52, 139–160.Mishra, N., Vondˇrejc, J., Zeman, J., 2016. A comparative study on low-memoryiterative solvers for FFT-based homogenization of periodic media. Journal ofComputational Physics 321, 151–168.Moulinec, H., Silva, F., 2014. Comparison of three accelerated FFT-basedschemes for computing the mechanical response of composite materials. In-ternational Journal for Numerical Methods in Engineering 97, 960–985.Moulinec, H., Suquet, P., 1994. A fast numerical method for computing the lin-ear and nonlinear mechanical properties of composites. Comptes Rendus del’Académie des Sciences Paris 318, 1417–1423.Mura, T., 1987. Micromechanics of Defects in Solids. Martinus Nijho ff , Dor-drecht.Nemat-Nasser, S., Hori, M., 1999. Micromechanics: Overall Properties of Het-erogeneous Materials. Elsevier.Press, W. H., Teukolsky, S. A., Vetterling, W. T., Flannery, B. P., 2007. NumericalRecipies: The Art of Scientific Computing, Third Edition. Cambridge Univer-sity Press.Schneider, M., Merkert, D., Kabel, M., 2017. FFT-based homogenization for mi-crostructures discretized by linear hexahedral elements. International Journalfor Numerical Methods in Engineering 109, 1461–1489.Shanthraj, P., Eisenlohr, P., Diehl, M., Roters, F., 2015. Numerically robust spec-tral methods for crystal plasticity simulations of heterogeneous materials. Inter-national Journal of Plasticity 66, 31–45.Suquet, P., 1997. Continuum Micromechanics. Vol. 377 of CISM InternationalCenter for Mechanical Sciences. Springer, Berlin.Trefethen, L. N., 2000. Spectral Methods in MATLAB. SIAM.Willot, F., 2015. Fourier-based schemes for computing the mechanical response ofcomposites with accurate local fields. Comptes Rendus Mécanique 343, 232–245. 33 . Trapezoidal approximation / discretization Note that the trapezoidal discretization ˆ f t κ of ˆ f κ in (9) satisfies cardinality (cid:88) m κ = − m e ı k κ x i ˆ f t κ = (cid:88) nj = (cid:88) ( n − / κ = − ( n − / n − e πıκ ( i − j ) / n f ( x j ) = (cid:88) nj = δ i j f ( x j ) = f ( x i ) , (cid:88) m κ = − m + e ı k κ x i ˆ f t κ = (cid:88) nj = (cid:88) n / κ = − n / + n − e πıκ ( i − j ) / n f ( x j ) = (cid:88) nj = δ i j f ( x j ) = f ( x i ) , (cid:88) m − κ = − m e ı k κ x i ˆ f t κ = (cid:88) nj = (cid:88) n / − κ = − n / n − e πıκ ( i − j ) / n f ( x j ) = (cid:88) nj = δ i j f ( x j ) = f ( x i ) , (A.1)(i.e., the interpolation condition) for both n odd ( m = ( n − /
2) and n even ( m = n / I o n f ( x ) : = (cid:88) ( n − / κ = − ( n − / e ı k κ x ˆ f t κ (A.2)for n odd, and I e n f ( x ) : = (cid:88) n / − κ = − n / e ı k κ x ˆ f t κ = (cid:88) n / κ = − n / + e ı k κ x ˆ f t κ (A.3)for n even, interpolate f ( x ). Note that the second form of (A.3) is employed forexample by Willot (2015, Equations (10)-(12)).For Fourier discretization based on (9), the index transformations κ (cid:55)→ µ = κ + + m : {− m , . . . , m } → { , . . . , n } , m = ( n − n odd ,κ (cid:55)→ µ = κ + + m : {− m , . . . , m − } → { , . . . , n } , m = n , n even ,κ (cid:55)→ µ = κ + m : {− m − , . . . , m } → { , . . . , n } , m = n , n even , (A.4)are useful. Restricting attention to the first two, define q µ : = ı k µ − − m , ϕ µ ( x ) : = n − / e q µ x , ˇ f t µ : = n − / (cid:88) ni = e − q µ x i f ( x i ) . (A.5)Based on these, (A.2) and (A.3) take the common form I n f ( x ) = (cid:88) n µ = ϕ µ ( x ) ˇ f t µ , (A.6)with ˇ f t µ = (cid:80) ni = F t − µ i f ( x i ), f ( x i ) = (cid:80) n µ = F t + µ i ˇ f t µ , and F t ± µ i : = n − / e ± q µ x i = n − / e ± πı ( µ − − m )( i − / n , m = ( n − n odd n n even . (A.7)Note that ϕ µ ( x i ) = F + µ i and (cid:80) n µ = F t + µ i F t − µ j = δ ij , consistent with cardinality.34 . Piecewise-constant approximation / discretization The form ˆ f c κ for ˆ f κ from (20) in the case of piecewise-constant f ( x ) determines thecorresponding Fourier series discretization F c n f ( x ) : = (cid:88) n µ = ϕ c µ ( x ) ˇ f c µ , f ( x c i ) = (cid:88) n µ = F c + µ i ˇ f c µ , ˇ f c µ = (cid:88) ni = F c − µ i f ( x c i ) , (B.1)via (A.4) , . Here, ϕ c µ ( x ) : = n − / s c µ e q µ x , s c µ : = s µ − − m and F c ± µ i : = n − / e ± q µ x c i . Since s c µ (cid:44) µ (cid:44) m +
1, note that F c n f ( x c i ) = (cid:88) n µ = ϕ c µ ( x c i ) ˇ f c µ = (cid:88) nj = (cid:88) n µ = s c µ F c + µ i F c − µ j f ( x c j ) ≈ f ( x c i ) (B.2)is only approximately cardinal. In the approach of Eloh et al. (2019) also basedon (20), approximate cardinality takes the form f ( x c i ) ≈ F c p f ( x c i ) = n − / (cid:88) n − ω = e ı k ω x c i (cid:88) m − ν = − m ( − ν s ν n + ω ˇ f c ν n + ω (B.3)via (25) and e − ı k ν n x c i = e − πıν ( i − ) = e πıν e − πıν i = ( − ν ( p = nm even). Comparisonof (B.3) and (26) results inˇ f c ω ≈ (cid:88) m − ν = − m ( − ν s ν n + ω ˇ f c ν n + ω , (cid:88) m − ν = − m ν (cid:44) ( − ν s ν n ˇ f c ν n ≈ , (B.4)via the fact that ¯ f = n − / ˇ f c0 = n − (cid:80) ni = f ( x c i ) from (1) , (8) and (20) . Given(B.4), note that (B.3) reduces to ˜ f ( x c i ) ≈ (cid:80) n − ω = e ı k ω x c i (cid:80) m − ν = − m ( − ν s ν n + ω ˇ f ν n + ω for thefluctuation part ˜ f ( x c i ) of f ( x c i ) = ¯ f + ˜ f ( x c i ) in the context of (1).As done by Eloh et al. (2019), the principle application of these relations and inparticular of (B.4) is to obtain the discretizationˇ Γ cH ω ˇ T c ω = , ω = , . . . , n − , (B.5)of (28) in terms of the "discrete Green operator" (DGO)ˇ Γ cH ω : = (cid:88) m − ν = − m s ν n + ω ˆ Γ H ( k ν n + ω ) , ω = , . . . , n − ..