Elimination of ringing artifacts by finite-element projection in FFT-based homogenization

Micromechanical homogenization is often carried out with Fourier-accelerated methods that are prone to ringing artifacts. We here generalize the compatibility projection introduced by Vond\v{r}ejc, Zeman&Marek [Comput. Math. Appl. 68, 156 (2014)] beyond the Fourier basis. In particular, we formulate the compatibility projection for linear finite elements while maintaining Fourier-acceleration and the fast convergence properties of the original method. We demonstrate that this eliminates ringing artifacts and yields an efficient computational homogenization scheme that is equivalent to canonical finite-element formulations on fully structured grids.


Introduction
Mechanical homogenization seeks the computation of effective mechanical properties, such as homogenized elastic constants or the complete (potentially nonlinear) stress-strain response, given the microstructure of material and constitutive laws for the individual components. Homogenization often operates on periodic "representative" volume elements that are characteristic for the material under investigation [1,2]. The intrinsic periodicity of such volume elements suggests the use of a spectral Fourier-basis [3] for the efficient numerical solution of the homogenization problem.
Since the seminal works by Moulinec and Suquet in 1994 [4] and 1998 [5], (fast) Fourier transform (FFT) based homogenization methods were developed rapidly; see, e.g. Ref. [6] for an authoritative review. The spectral method is here used to solve the partial differential equation (PDE) for a static mechanical equilibrium of a microstructure. The availability of highly optimized FFT implementations (like FFTW [7]) enabled the implementation of efficient and highly parallel spectral solvers which can beat standard implementations of the finite element method (FEM) in speed and accessible system sizes [8].
A major advantage of spectral methods is that it is straightforward to increase the basis set and by this systematically improve the accuracy of the solution. However, because of their global support trigonometric polynomials are not well suited for the solution of PDEs with discontinuous material coefficients or discontinuities in their solution. This leads to a phenomenon known as Gibbs ringing. Gibbs ringing in Fourier spectral methods is well documented, and there are several approaches to suppress the problem [9].
Gibbs ringing persists in FFT-based homogenization methods and was investigated by several authors (see e.g. Refs. [10,11,12,13,14,15,16,17]). A typically mitigation strategy is the introduction of discrete derivatives: for example, Müller [10] used a finite difference discretization or Willot [14] and Schneider et al. [15] describe a central-difference scheme on a staggered grid. However, published works have mainly focused on natural extensions of the Moulinec-Suquet scheme that is formulated in terms of a reference problem in the displacement space, i.e. for small-strain elasticity in form of the Navier-Lamé equations, that explicitly contain the constitutive law. A solution in Fourier space is then only possible for a homogeneous reference medium. Heterogeneity is superimposed onto this reference solution through Lippman-Schwinger-like approaches for capturing it. Convergence properties of such schemes deteriorate with the strength of the heterogeneity.
In this article we discuss the problem of ringing in FFT-based homogenization methods and present a general solution to it for which we combine the advantages of an FEM-type Galerkin approach with the speed-up of a FFT method. We work in the framework of compatibility projection developed by Vondrejc, Zeman, de Geus and coworkers [18,19,20] that considers the deformation gradient as the primary degree of freedom and eliminates the need for an explicit reference medium. Previous work, e.g. Refs. [21,22], has proposed finite differences approaches to reducing ringing artifacts. Different to previous works, we derive a general expression for a discretized operator for the projection on compatible fields and we propose an elimination rather than a mitigation of the problem for certain operator choices. We show that in addition to Gibbs ringing, the solution with projection techniques is prone to ringing artifacts arising from missing degrees of freedom in the formulation of the deformation gradient.
Our suggested improvements to the method are supported by numerical simulations and comparison with analytical results. We suggest a linear finiteelement-based projection operator that eliminates all ringing artifacts and discuss the analogy of the modified FFT-based homogenization method to standard

Compatibility projection
We investigate microstructures by a representative volume element (RVE) in a periodic simulation cell Ω 0 . The microstructure can consist of different phases which are described by arbitrary small-or finite-strain material models. Here and in the following we will denote first-order tensors (vectors) by arrows and second-order tensors by bold symbols.
Any deformation of the simulation cell can be described by a function (the placement map) χ : Ω 0 → Ω, mapping the undeformed grid positions r ∈ Ω 0 into a deformed configuration χ( r) ∈ Ω where Ω is the deformed periodic simulation domain. The overall goal is to solve for the static mechanical equilibrium of the periodic cell for a given deformation. The static mechanic equilibrium is given by [24] ∇ 0 · P T (F( r)) = 0 or ∂ 0,α P iα (F( r)) = 0, with implicit summation over repeated indices (the Einstein summation convention) and a dot product (·) as defined by the second part of Eq. (1).
Greek indices indicate a tensor dimension representing a derivative, whereas Latin indices are used for all other tensor dimensions. P is the first Piola-Kirchhoff stress tensor, in general a non-linear function of the deformation gradient Note that the partial derivatives ∂ 0,α ≡ ∂/∂r α in Eqs. (1) and (2) are with respect to the undeformed configuration Ω 0 of the cell. The most widely used homogeneization schemes combine Eqs. (1) and (2), leading to a set of second-order differential equations in χ. For small strains, those are the well known Navier-Lamé equations that contain spatial derivatives of the elastic constants. The Green's function of this second order differential equation therefore explicitly contains the constitutive law. Solution with Fourier-techniques then requires the introduction of a homogeneous reference material.
In contrast, the formulation employed here [18,19,20] solves the set of first-order differential equations given by Eqs. (1) and (2). This leaves the deformation gradient F as a degree of freedom. Equation (2) can be interpreted as the constraint that F needs to be compatible, i.e. given by the gradient of a respective placement map. The mathematical trick is that this constraint can be implemented by a projection operator G that turns any second-order tensor into a compatible one.
Following Refs. [18,19,20], we reformulate Eq. (1) in the weak (weighted residual) form, where t( r) is an arbitrary periodic (vector-valued) test function, ⊗ the outer product, D the spatial dimension and : the double dot product, a tensor contraction over two indices, A : B = A ij B ji . Surface terms that should appear in Eq. (3) vanish due to periodicity. If we interpret the test function t( r) as a displacement, then δF( r) = ∇ 0 ⊗ t( r) is a suitable set of compatible test gradients. We can therefore write the equilibrium condition as where now δF is an arbitrary (no longer necessarily compatible) test function and G A denotes the application of the self-adjoint operator G to a right-hand side object A. We now discretize the gradients rather than the displacement field, i.e. in the spirit of the Galerkin method, we express the test gradient δF and the deformation gradient F within the same basis set.
Equation (4) no longer contains gradients of the constitutive law, given by P(F). The compatibility operator G is clearly independent of it since it just ensures fulfillment of the compatibility condition, i.e. ∇ 0 × (G : δF) = 0 for finite strain formulations. The numerical scheme can therefore be formulated without the explicit requirement of a homogeneous reference material.
The compatibility operator G is block-diagonal in the Fourier basis, where the blocks are fourth order tensors. This leads to the expression [19] whereĜ is a fourth-order tensor and k the wavevector. We use the hat symbol to denote a Fourier-transformed quantity and F{·}( k) for the explicit Fourier transform of a quantity as given in Appendix A by Eq. (A.3). The numerical solution of Eq. (5) benefits from the fact that for gradient-based optimizers, the steps of the optimization procedure automatically lead to compatible F, see [20]. The operator G projects arbitrary tensor fields onto compatible fields, i.e. fields that can be expressed as the gradient of a lower order tensor. For enforcing the compatibility for second-order tensor fields, the operator is given by [18] where k α is the component α of the wavevector k and k = | k|.

Interpreting the projection operator
The projection operator lends itself to a simple interpretation. Let us assume we have an arbitrary vector field v( r). This field is only a gradient v( r) = ∇ 0 φ( r) of a scalar field φ( r) if its curl vanishes, ∇ 0 × v( r) = 0. If the field is also periodic, we call these fields compatible. In the context of homogenization, v is a row of the deformation gradient F and φ a component of the placement map χ. We are interested in the special case where v( r) is periodic but φ( r) is not. The periodicity of v will be later intrinsically fulfilled through the Fourier transform and we only investigate the gradient property here. For non-compatible fields, we search for the scalar field φ( r) that minimizes the residual vector in a suitable sense (for compatible fields, R ≡ 0). For a minimal residual vector ( R( r) → min) Equation (7) is an equivalent formulation of Eq. (2) and the operator D −1 introduced below is the Green's function of that equation. The canonical requirement is minimization in the least-squares sense, i.e. minimization of We now need to choose a specific basis set for a series expansion of φ( r). In a Fourier basis, and an equivalent expansion holds for v, with v( 0) = v 0 and N the total number of voxels, see Appendix A. Note that in Eq. (9) we have set the (arbitrary) mean value of φ( r) to zero but added a linear function that cannot be represented in the Fourier basis whose derivative gives the mean value of Eq. (10). In terms of homogenization, this mean value describes the affine deformation of the whole RVE and plays the role of the boundary condition of the constituting differential equation.
In the Fourier space the residual becomes for k = 0 whereˆ D( k) = i k is the Fourier representation of the gradient ∇ 0 . For where the star is the complex conjugate. Note that all symbols in Eq. (12) -ˆ v, D andφ -are function that depend explicit on k but that dependence has been omitted here for brevity. Minimization gives the secular equation We can solve this for (note that k = 0) where we interpret the termˆ as the inverse of the derivative, i.e. as some form of "integration". In a next step, we compute the gradient ofφ( k). This yieldŝ The operatorĝ( k) hence projects an arbitrary field on its compatible part in the least squares sense with respect to the integral inner product (L 2 -norm). The full projection operator for the deformation gradient is given by the first part of Eq. (6). Using the Fourier derivativeˆ D( k) = i k yields the specific form of the projection operator given in Eq. (6). Since G contains the Green's function of Eq. (2), but not of Eq. (1), it is independent of the constitutive law and does not require a reference material. The formulation of the projection operator for small-strain elasticity is described in Appendix B.

Discrete projection
Rather than using Eq. (9), we can expand φ( r) in other bases of choice. For example, we will below employ linear finite elements. Other discretizations of the gradient operator ∇ 0 with less suitable properties can be obtained through finite-differences methods. We here assume that the simulation cell is structured in a regular (equally spaced) grid with node positions { r IJ }, where I and J are node indices. We will call the individual grid cell a voxel and develop the theory in two dimensions, but generalization to three dimensions is straightforward. The placement map χ is only known at the nodes which are the corners of the voxels.
The discrete derivative (in some "direction" α) in voxel I, J can generally be written as the convolution with r I+i,J+j = r IJ + r ij and ∆ (α) the voxel size in direction α where the brackets and the upper index indicate that there is no Einstein sum convention applied for the index (α) and periodicity of r in a natural sense is considered. The collection of coefficients s ij is called the stencil of the operation. We will introduce different stencils below and note that derivatives in different directions α require different stencils. In the following we will assume that the deformation gradient is described by a set of d derivatives, but that d is not necessarily equal to the dimension D of the system. This allows the subdivision of voxels into multiple evaluation points, i.e., d may be an integer multiple of D. Using n q evaluation points per voxel leads to d = n q D derivatives per voxel within the framework that we now develop. The subscript () q was chosen because these evaluation points correspond to Gaussian quadrature points in a classic finite element discretization. By expanding the discrete placement map into a Fourier series of the form given by Eq. (A.4), we can write the derivative operation in Fourier space aŝ Again the case k = 0 is special since ij s (α) ij = 0. For k = 0, we can use the same argument as above: A general vector fieldˆ v( k) can be projected (in the least squares sense) onto its compatible partˆ w( k) usinĝ withĝ cf. Eq. (15). We want to mention that by the Helmholtz decomposition δ αβ − g αβ ( k) is a projection to divergence free fields, with applications to error estimation. The projection operator for the deformation gradient is then given byĜ iαβj ( k) = δ ijĝαβ ( k). Note thatĝ ∈ C d×d andĜ ∈ C D×d×d×D . The projection operator is idempotent (i.e. a projection), since, This compatibility operatorĜ is the projection generalized for arbitrary derivative operatorsˆ D.
The case k = 0 contains the boundary condition and is treated as for the Fourier derivative. For multiple elements, each element needs to hold the same gradient for k = 0. We want to emphasize that in this discrete basis, the discrete Fourier transformation has the role of accelerating the convolution rather than providing the basis set for the underlying discretization.

Finite differences and Lanczos-σ correction
Canonical discrete derivative operators are given by finite difference schemes. We here consider the first-order central differences, which can be expressed in Fourier space as follows (see e.g. [22]): Remember that ∆ (α) is as before the grid spacing in direction α and there is no summation over upper indices in parenthesis. We note that this central differences scheme is related to the Lanczos-σ correction for Gibbs ringing in direction α [25], . Expressed using the σ-factor, the derivative operator is given by D cd . Another common first-order finite difference scheme is forward differences, with the Fourier-space representation The full stencil coefficients for these two schemes as applied to two-dimensional problems are shown in Fig. 1a and 1b. Inserting these coefficients into the generic expression Eq. (19) yields the specific Fourier representations given in Eqs. (24) and (26). Figure 1 also has a graphical representation of these discrete derivatives.

Least squares
We now turn to a purely geometric interpretation of the deformation gradient to derive alternative discrete stencils. The deformation maps an infinitesimal fiber vector δ r into δ r = F( r) · δ r. Assuming constant F over a given voxel (light blue rectangle in Fig. 2a), the voxel is affinely deformed by the deformation gradient F (green parallelogram in Fig. 2b) and cannot represent arbitrary displacements of the corners (dark blue trapezoid in Fig. 2b) as long as the (a.1) central differences forward differences least squares linear finite elements deformation gradient F is uniform on that voxel. In order to represent the corner displacements exactly, we can for instance subdivide the voxel, e.g. in 2D into two triangles (see Fig. 2c and Fig. 2d), introducing multiple elements per voxel, with their own uniform deformation gradient per element. We will discuss this decomposition in the next section and will for now continue to work with a uniform deformation gradient per voxel and require matching of the corner displacements in a least-squares sense.
We now regard the displaced corner positions only within a single voxel (see Fig. 2a). The undeformed coordinates of the voxel are r kl with k, l ∈ {0, 1}. The true deformed coordinates are χ kl and the affinely deformed coordinates produced by the deformation gradient F are ψ kl . We require that the affinely deformed rectangle matches the true deformation in a least square sense. This means we are looking for the deformation gradient F that minimizes the residual Since the affinely deformed voxel is spanned by the two vectors F · ∆ r 10 and F · ∆ r 01 with ∆ r kl = ( r kl − r 00 ), we can write the affinely deformed coordinates as We now insert these into Eq. (27) and minimize with respect to F and ψ 00 . This yields the secular equations ∂R or and For rectangular lattices with lattice spacing ∆x and ∆y, this can be solved to give The corresponding stencil coefficients are shown in Fig. 1c.

Linear finite elements
The previous section argued that a uniform deformation gradient per voxel is insufficient to represent the voxel's deformation as measured by the displacement of each corner. This becomes evident from simply counting the degrees of freedom: In two dimensions, the voxel's deformation (and rotation) is given by three vectors (6 degrees of freedom) while the deformation gradient has 4 independent components. In this section, the problem is solved by splitting the voxel into two triangles that are each described by a uniform deformation gradient (see Fig. 2c). The voxel still has 6 degrees of freedom, but now we have 2 deformation gradients with a total of 8 independent components and the constraint that the diagonal boundary between the two triangles remain of the same length and direction (two blocked degrees of freedom). From a simple geometric argument, this triangular decomposition can hence describe arbitrary corner displacements (see Fig. 2d).
This geometric point of view is fully equivalent to a formulation using linear finite elements. Within our rectangular lattice, we use the usual linear shape functions 01 (x, y) = y/∆y (38) where the origin of the coordinate system is at the bottom left of the voxel and an equivalent set of shape functions exists for element (2) in Fig. 2c. The shape function gradient of 01 (x, y) is constant on element (e) and given by the stencil coefficients shown in Fig. 1d for the two deformation gradients. The stencil for element (1) is identical to the forward differences scheme. The stencil for element (2) is a forward differences scheme evaluated on a different set of nodes turning it into a backward differences scheme. Generalizations of this scheme to non-orthogonal voxels and three dimensions are straightforward. This formulation is identical to traditional linear finite elements, but unlike classical finite-elements, the projection formulation yields a condition number that is independent of system size, leading to scale-independent convergence properties [26]. Note that the least square approach described in the previous section simply yields the average of the deformation gradients on the two triangles.
It is important to emphasize that this formulation requires storing two deformation gradients per voxel. In the discrete projection developed above, this means we have n q = 2 evaluation points and d = 4 derivatives for the twodimensional formulation. The deformation gradient F and Piola-Kirchhoff stress P are then both elements of R D×D×nq within the projection scheme. For the evaluation of the constitutive law, both F and P need to be decomposed in their element-wise contributions F (e) and P (e) that are elements of R D×D . We can formally write for this decomposition. We note that finite element discretizations with multiple quadrature points can be described in this discrete projection using similar decompositions. The interpretation put forward at the beginning of this section is that the deformation gradient describes the geometry of the respective triangle (see Fig. 2d). This allows an intuitive interpretation of the projection operator G. We start with a decomposition of some domain into triangles (Fig. 3a). The rotation and shape of each of these triangles is described by a deformation gradient F (e) ∈ R D×D . We now randomly disturb these triangles by adding a random number to the components of their deformation gradients. The resulting structure (see exploded view in Fig. 3b) is clearly no longer compatible. Application of G (Fig. 3c) adjusts the shape of each triangle such that they are again compatible and can be assembled into a continuous deformed structure (Fig. 3d).

Even vs. odd number of grid points
The original formulation of the compatibility projection that employs the Fourier-derivative only works exactly for odd-sized grids, which is a result of the structure of the projection operator. The Fourier derivative is ambiguous at the Nyquist frequency (or the edge of the first Brillouin zone). This ambiguity originates in the freedom of choosing one of two possible equivalent even numbered Fourier grids {k i } from the class k i ∆ ∈ [−π, π) or k i ∆ ∈ (−π, π], where ∆ is the grid spacing. Even sized grids sample the frequency exactly at the Brillouin zone boundary (see Eq. (A.2) for the choice k i ∆ ∈ [−π, π)). As a result the two possible grids differ only in one single grid point, the Nyquist frequency.
The Nyquist frequency is given by k Ny = ±π/∆, where an even sized grid contains either the positive or negative Nyquist frequency. Typically this small difference does not matter for the periodic functionf (k) because the Fourier coefficients are also equivalent in the single point ±k Ny where the two grids differ,f (−k Ny ) from the one grid containing −k Ny has the same value asf (+k Ny ) from the other grid containing +k Ny . Thus, the Fourier series does not depend whether the positive or negative Nyquist frequency is included in the grid and is therefore unambiguous for any (sampled) function f (x). However, when analyzing the Fourier series of D F f (x), where the upper index "F" indicates the explicit use of the Fourier derivativeD F (k) = ik, evaluated on grid points x, we find for the summand S(k Ny , x) at the Nyquist frequency where we used k Ny x ∝ πn, n ∈ Z. SinceD F (k Ny ) = ik Ny = −ik Ny =D F (−k Ny ) the summand S(k Ny , x) = S(−k Ny , x) which gives an ambiguity in the Fourier series of D F f (x) for even sized grids at the Nyquist frequency. The outcome of the Fourier series in Eq. (40) depends on which one of the two possible even sized grids is taken. This ambiguity vanishes for odd sized grids, since the function is never evaluated at the Nyquist frequency and there is only a single possible choice for the Fourier grid. Figure 4 plotsD F (k) at the discrete evaluation points for even and odd-sized grids as an illustration of this problem. This ambiguity is resolved for all the discrete stencils. For example, the forward differences stencils fulfillsD fd (k Ny ) =D fd (−k Ny ) (see Fig. 4). Since we are interested in the compatibility projection operator we additionally require a second conditionD(k Ny )D * (k Ny ) = 0 for a proper derivative operator on even grids to prevent a division by zero in the computation of the compatibility operator by Eq. (21), or in other words we need a properly defined inverse derivative operatorD −1 (k) (see Eq. (15)). An exception from the presented discrete stencils is the central differences stencil, whereD(k Ny )D * (k Ny ) = 0 as a result of a decoupling of the domain into two subgrids. An alternative point of view is that the discrete Fourier-transform is simply employed to carry out a discrete convolution and there is no ambiguity with respect to the Nyquist frequency.

Examples and validation
In the following, we describe four two-dimensional examples to demonstrate the methods developed above with a focus on ringing phenomena. First, we investigate a single soft voxel in a uniform hard matrix under biaxial strain. This minimal example already shows strong ringing artifacts in the xy-component of the stress tensor. Second, we analyse a cell with two pillars separated by one layer of voxels with the Young modulus set to zero. This example shows the ability to handle infinite material contrast with vacuum [27] and simulate a free surface despite the intrinsic periodic boundary conditions 1 . Additionally, one of the pillars contains an inhomogeneity which gives rise to ringing artifacts in the original spectral formulation. Third, we test the numeric correctness of the results by investigating an Eshelby inhomogeneity. We compare the results obtained by the FFT-based method against the analytical Eshelby solution, corrected for periodic boundary conditions. Finally, we demonstrate the feasibility of the method for complex constitutive laws on a damage mechanics problem.

Single voxel inhomogeneity
A classical continuum mechanics problem is the inhomogeneity, an inclusion of a material in a matrix with different material properties. At the boundary of the inhomogeneity, there is a discrete change in the material properties which usually leads to Gibbs ringing artifacts in spectral methods. As minimal example of such an inhomogeneity, we present a single voxel inhomogeneity (in red) placed in the center of a 17 × 17 voxel matrix (in green) as shown in Fig. 5a. The matrix with Young modulus E hard is ten times harder than the inhomogeneity E soft = E hard /10, and both have the same Poisson ratio of ν = 0.33. The material response is described by an isotropic finite strain linear elastic law as described in Refs. [24,20].
We apply a biaxial tensile strain of 10% (F xx = F yy = 0.1 and F xy = F yx = 0) and investigate the shear component P xy of the first Piola-Kirchhoff stress for different implementations of the projection operator. Figure 5 gives an overview of the results where the first row, panels (b.1) to (g.1), show the color coded stress and the second row, panels (b.2) to (g.2), give a more detailed look on the stress along selected rows of the matrix as indicated by the colored arrows in between the subfigures in the first row of Fig. 5. The lines in the panels (b.2) to (g.2) are ordered from top to bottom starting from the center row, i.e. the row with the inhomogeneity in red. Each column represents the result found with a different projection operator: Column (b) for the Fourier-type projection operator as given by Eq. (6), (c) for central differences given by Eq. (24), (d) for forward differences given by Eq. (26), (e) for the least square scheme described in Sec. 2.5, (f) for the Fourier-type projection operator on two evaluation points per voxel (see Appendix C) and (g) for linear finite elements on two elements per voxel as derived in Sec. 2.6.
Panel (b) show the stress field for the original method (e.g. Ref. [20]), with a projection operator based on the Fourier derivative. As expected we observe strong ringing artifacts leading to a checkerboard pattern of the stress field. The stress field and its oscillations are strongest at the inhomogeneity and decay with increasing distance to the discontinuity. However, the ringing does not disappear even at the edges of the cell which was also tested for finer grids, different material contrasts and a slightly inhomogeneous matrix (results not shown). The symmetry of the setup leads to a line of zero stress in the row and column that contains the inhomogeneity (see Fig. 5b.1 and red line in Fig. 5b.2).
Results obtained with the central differences projection operator are shown in Fig. 5c.1 and c.2. The Gibbs ringing artifacts should be strongly suppressed for this method, however we observe a checkerboard pattern of different style compared to (b). This checkerboard pattern originates in the well-known [28,29] (odd-even) decoupling into two subgrids over short distances of the central differences stencil shown in Fig.1a. The decoupling of the two grids is not complete because of an odd-sized grid (17 × 17) and the oscillations decay with increasing distance from the inhomogeneity.
The forward differences stencil, results shown in panels (d), leads to an oscillation-free but slightly asymmetric stress field which originates from the asymmetry of the forward differences scheme (see Fig.1b). This asymmetry is corrected by the least square stencil shown in panels (e). However, the stress has also a checkerboard pattern with checkerboard characteristics of (b.1) and (c.1). We note that the reason for this ringing artifact is neither the Gibbs phenomenon nor the decoupling of two subgrids but the fact that the least-squares derivative cannot represent arbitrary deformations of the voxels as discussed in Sec. 2.5. This discussion, and the outcomes from (b) to (e), indicate that a symmetric and ringing free stress field cannot be obtained by a method with a single deformation gradient per voxel.
Therefore, we also investigated projection operators evaluated on two evaluation points per voxel. For the Fourier type projection operator on two evaluation points per voxel, as described in Appendix C, we still observe ringing, see panel (f). However, ringing is reduced with respect to the Fourier derivative on a single evaluation point. In panel (f.2) the continuous line represents the stress values in the lower triangle (P (1) xy ) and the dashed line represents the stress in the upper triangle (P (2) xy ) (cf. Fig. 2c). In difference to the stress fields computed with a single evaluation point per voxel the two evaluation points per voxel slightly break the symmetry of the problem. This can be seen in the non-zero stress along the row of the inhomogeneity shown by the red line in (f.2).
Finally, we find a ringing free stress field for the discrete projection operator obtained from linear finite elements on two elements per voxel, panel (g). The asymmetry between the two elements of a voxel discussed in the previous paragraph persists for this formulation. The stress field result presented in panel (g) seems to be the most appropriate solution to the problem due to its smooth, ringing free field. This conclusion will be supported by the following three examples. We would like to again emphasize that ringing in this example does not only result from the Gibbs phenomenon, which does not exist for a discrete projection operator (as is evident for the forward differences in panel (d)). A description of local deformation with too few degrees of freedom as shown (e.g., the least square type projection in (e)) also gives rise to ringing.

Two pillars and vacuum
We use a setup consisting of two pillars, as shown in Fig. 6.a, to qualitatively investigate an infinite material contrast and the ability to represent a free surface. This would allow breaking the periodic boundary conditions which are intrinsic to FFT-based methods. We choose a simulation domain of 17 × 17 voxels and subdivide it into two pillars (in green) with Youngs modulus E hard and Poisson number ν = 0.33. The two pillars are separated by a layer consisting of single voxel of a material with zero stiffness (in light blue), i.e. E vac = 0. One Figure 6: The first Piola-Kirchhoff shear stress Pxy of two pillars at 10% strain in y-direction. (a) displays the phase setup. The left pillar (green) has an inhomogeneity (red) of three voxels in its center which are ten times softer than the rest of the pillar. The right pillar (green) is separated by two layers of zero stiffness material (blue) from the left pillar. The first row (x.1) presents the shear stress field for the full grid. The second row (x.2) presents the shear stress as function of x for selected rows of the simulation grid indicated by colored arrows in the first row (x.1). The rest of the figure is illustrated as described for Fig. 5.
can think of this material as "air" or "vacuum". At the surface of the pillars we thus have an infinite material contrast. The left pillar, of width 7 voxels, has additionally an inhomogeneity (in red) of three voxels in its center. The soft inhomogeneity has a Youngs modulus of E soft = E hard /10 and the same Poisson number of ν = 0.33 as the pillar. The inhomogeneity was introduced to generate a non-homogeneous strain field with the ringing artifacts documented in the previous section. The setup is strained by 10% in the y-direction and held at constant size in x-direction. Simulations use the same finite strain model as those of Sec. 3.1. Figure 6 is organized in the same manner as Fig. 5. In the first row, panels (b.1) to (g.1), we show the color coded stress field P xy . In the second row, panels (b.2) to (g.2), we present the stress field as function of χ x at fixed positions χ y , starting from the center of the inhomogeneity (red line) and going row by row down in orange, green and light blue. The gray vertical lines represent zero stress for each curve. The continuous and dashed lines in panels (f.2) and (g.2) represent the stress of the lower and upper triangular element of a voxel, respectively. The columns are ordered as in Fig. 5: (b) Fourier type projection, (c) central differences, (d) forward differences, (e) least squares projection, (f) Fourier-type projection on two evaluation points and (g) linear finite element projection on two elements.
For the Fourier-type projection operator in panel (b), we observe ringing artifacts in the left pillar originating from the inhomogeneity. These artifacts are transmitted through the vacuum region to the right pillar. The shear component of the stress should be zero in the right pillar but shows clear ringing artifacts. The "vacuum" region of zero stiffness is therefor not able to decouple the two pillars. At the symmetry axis in xand y-direction of the inhomogeneity we can again observe a region of zero stress of single voxel thickness.
Central differences, panel (c), lead to a strong decoupling of the grid into two sub grids. The vacuum cannot decouple the strain field in the pillars because the stencil has a range of three voxels. This leads to strong oscillations in xdirection in the two subgrids. In the right pillar, that has a width of an even number of grid points in the x-direction, the decoupling results in almost zero width of the voxels of one sub grid. (Only four voxels are clearly visible in panel (c.1).) Remarkably the left pillar has no ringing in y-direction in the columns of non zero stress. In this example, it becomes very clear that the central differences are sensitive to the setup and number of grid points. For slightly different widths of the pillars or a different mesh grid one can observe these artifacts also in the other pillar (not shown). Panels (d) and (e) (forward differences and least-squares), show the same behavior as discussed for a single voxel inhomogeneity in the previous Sec. 3.1 with an asymmetric, but non-ringing response for the forward differences and a symmetric and ringing response for the least squares scheme. However, the right pillar shows zero shear stress, indicating that for these projection operators the vacuum layer is able to decouple the two pillars. The vacuum layer grows in x-direction to absorb the shrinkage of the two pillars (since ν = 0.33) while the average strain in x-direction is zero. At the surfaces of the left pillar the shear stress field goes to zero as one would expect it for a surface. These discrete derivative schemes decouple the two pillars because their stencil extend only between neighboring nodes.
The Fourier type projection on two evaluation points in panel (f) leads to similar results in the left pillar as for the single voxel inhomogeneity. A ringing artifact from the left pillar is observed in the right pillar which indicates a coupling between the two pillars through the vacuum region. As in the first example, Fig. 5f, the two evaluation points lead to a slight symmetry breaking resulting in non-zero stress at the symmetry axis through the center of the inhomogeneity in the left pillar, i.e. the red line in panel (f.2).
For the linear finite element projection shown in panel (g), we again observe artifact-free results. The left pillar shows a smooth, ringing free and symmetric (besides the previous discussed slight asymmetry between upper and lower triangle) stress field originating from the inhomogeneity. The vacuum regions grow in x-direction to absorb the shrinking of the pillars in that direction. In the right pillar we find no artifacts from the stress field of the left pillar; thus the pillars are fully decoupled.
This example shows that it is possible to simulate infinite material contrast with vacuum as the soft phase. A single layer consisting of material of zero stiffness can decouple different regions in the RVE and thus break the intrinsic periodic boundary conditions of the FFT-based method in one direction. As expected from the theoretical considerations in Sec. 2 the linear finite elements projection operator has the best performance and results in a stress field that is artifact free and qualitatively correct. To further investigate the methods developed here, we continue with a quantitative analysis of an Eshelby inhomogeneity.

Eshelby inhomogeneity
The Eshelby inhomogeneity is similar to the first example of a minimal inhomogeneity consisting of a single voxel. The Eshelby inhomogeneity is an ellipsoidal body inside an infinite elastic medium where the elastic medium differs in its material properties from the ellipsoidal body. The analytical solution to the Eshelby problem is well known [30,31,32,33]. We choose the specific (cylindrical) geometry shown in Fig. 7a with a hard matrix (light green) of Youngs modulus E hard and a soft inhomogeneity (red) of Young's modulus E soft = E hard /10 and zero eigenstrain. The numerical calculations employ a fine mesh of 151 × 151 voxels to properly resolve the cylindrical inhomogeneity with half axes of 10% of the domain edge lengths. The inhomogeneity is placed in the center of the domain and centered on a voxel to retain a symmetric discretized area. For the numerical calculations we use the small-strain formulation (see Appendix B) since the analytical Eshelby expressions are also obtained in this limit.
The results of these calculations are summarized in Figure 7. Rows (c.1) to (i.1) show the full solution of the shear strain ε xy on the 151 × 151 grid. The next row, panels (c.2) to (i.2), show a zoom of the region containing just the inhomogeneity. The three colored lines at the center (dark green), upper half (orange) and lower half (purple) of the inhomogeneity in panel (a) indicate the location of the strain components that are shown in the third row, panels (c.3) to (i.3), for the normal strain ε xx in x-direction at the center and in the fourth row, panels (c.4) to (i.4), for the shear strain ε xy for the upper and lower half of the inhomogeneity. Note that this data is shown only over the zoomed region, not the full calculation. The columns (c) to (h) again represent results obtained for different projection operators as indicated in each column: (c) is the Fourier-type projection, (d) central differences, (e) forward differences, (f) the least square type projection, (g) the Fourier type projection on two evaluation points per voxel and (h) the linear finite elements projection from sec. 2.6.
The additional column (i) represent the analytical results of the Eshelby inhomogeneity. The analytical result is obtained for an inhomogeneity in an infinite media at 1% biaxial strain (ε xx = ε yy = 0.01, ε xy = 0) and is corrected for periodic boundary conditions by summing up the influence of 11 × 11 noninteracting periodic Eshelby inhomogeneities (see panel (b) of Fig. 7). This correction scheme for periodic boundary conditions converges quickly with the number of images. The average strain on the central inhomogeneity is used as the boundary conditions for the periodic numerical computations.
The Fourier-type projection operator gives qualitatively correct results. Panels (c.3) and (c.4) show a quantitative view of two components of the strain tensor that both clearly show oscillations within the inhomogeneity. As expected, we observe strong ringing artifacts which lead to a deviation from the analytical result especially within the inhomogeneity. Note that the normal strain along the center line vanishes due to symmetry reasons and hence does not show ringing.
For the central differences scheme in panels (d.i), we observe a symmetric checkerboard pattern of decreasing amplitude when approaching the center of the inhomogeneity. At the boundary of the inhomogeneity, there is a double ring-like pattern in the strain field originating from the local decoupling into two subgrids by this scheme. The forward differences scheme in column (e) produces ringing-free fields but with the drawback of the already discussed asymmetry, best shown in panel (e.3). However the asymmetry of the field is partly suppressed by working in the small strain limit where ε xy = ε yx . The asymmetry of the strain field is also noticeable in panel (e. 4).
Column (f) shows the results produced by the least squares projection operator. As for the two previous examples we observe a checkerboard like pattern of decreasing intensity when approaching the center of the inhomogeneity (see panel (f.2)).
The Fourier-type projection on two evaluation points per voxel (panel (g)) produces a strain field similar to the standard case with a single point per voxel (panel (c)), however the ringing artifacts appear distributed more homogeneously over the entire inhomogeneity. Panels (g.3) and (g. 4) show for clarity only the strain field for the lower triangular element (element one in Figure 2c). In summary, we find reasonable agreement with the analytical solution for all projection operators. Linear finite elements gives the smoothest curves and results closest to the analytical findings. As in the previous cases, only the forward differences projection and the finite elements projection eliminate the ringing artifact. For a finer grid we observe a convergence towards the analytical result. It is worth noting that the similar behaviour of the finite elements and forward differences projections is easily explained by the fact that the latter corresponds to the finite element projection where only the lower left triangles are considered.

Damage problem
As a drawback of FFT-based solution methods, ringing artifacts can have a drastic effect on the solution of a homogenization problem. Damage mechanics problems are especially vulnerable to fluctuations in the stress field caused by ringing artifacts, since localization is one of the most important characteristics of such problems. A reliable, fast, and ringing-free homogenization method is therefore essential to address damage mechanics problems. In order to illustrate the error introduced by the ringing artifact into a damage problem, we solve a two-dimensional problem representing a concrete microstructure using an Alkali-Silica reaction (ASR) damage model. In this model, the expansion of gel pockets inside aggregates damages the microstructure.
In the modeled concrete microstructure, three phases are considered -a soft matrix with damage (Cement paste), hard inclusions with damage (Aggregates), and gels whose expansion is modeled by a growing spherical eigenstrain. Polygon-shaped aggregates are placed in the RVE using the level set approach (LSA) algorithm [34,35], with an aggregate size distribution chosen to match sieve sizes from real concrete structures. The aggregate size distribution is truncated on the lower end in order to keep the shape of the aggregates physically sound considering the descritisation grid. ASR gel pockets are placed randomly within the aggregate to fill 2% of the cell surface. Concrete paste and aggregate are represented by a linear damage model as their constitutive law [36,37,38]. Due to the brittleness of concrete, the damage part of the bi-linear damage laws is taken steeper than its elastic part.
In the damage phase, the damage surface threshold is defined by the magnitude of strain measured by the L 2 -norm. As long as the damage material's strain is below a determined u , it behaves as a linear elastic isotropic material with Young modulus of E 0 , and afterwards its stress decreases with equivalent stiffness of −αE 0 until the stress becomes zero. From that point on, the material does not carry any stress (complete failure). An eigenstrain with the final amplitude of 20 u was applied on gel voxels, placed as explained beforehand, in 1000 consecutive steps in the carried out simulation. The damage field caused by this loading scenario is depicted in Fig. 8 employing Fourier and FEM projection solvers. As observed in Fig. 8 the damage pattern evolved in the Fourier projection solver solution is checker-boarded and therefore non-physical. While as demonstrated in Fig. 8b, after damage initiation around ASR gel pockets, micro-cracks formed during the damage process coalesce to form cracks with lengths in the range of 0.2 times the RVE length. Comparison of Fig. 8a, b suggests that, in contrast with Fourier projections solvers, ringing-free spectral FEM projection solvers are capable to simulate mechanics damage RVE problems.

Summary & Conclusion
We have extended the compatibility projection of Vondřejc, Zeman, de Geus and coworkers [18,19,20] to basis sets with local support. The nonlocal support of the Fourier basis gives rise to ringing artifacts and makes it impossible to include regions with zero stiffness or vacuum into a calculation. We show that using linear finite elements as a basis set eliminates all of the problems of the Fourier basis but retains the advantages of the compatibility projection, such as the rapid convergence rate. In particular, this formulation allows modeling of regions with zero stiffness, which opens the method to the study of free surfaces or metamaterials. We note that the ideas behind compatibility projection can be exploited to construct preconditioners for standard finite-element formulations beyond the linear elements presented here [39].
Among the projection operators examined in this paper, all but the least squares operator guarantee a compatible strain field. It is useful to think of those compatible projection operators in two categories; discretisations with local support and discretisations with global support, where local support refers to non-overlapping projection stencils. We observe that only two local support projections, the finite element projection and the forward differences projection eliminate ringing, unlike the Fourier and central differences operators (global support and non-local support including neighbouring voxels). This observation allows a well known property of finite element calculations, where discretisation errors are to be expected to be significant at the length scale of the element size (= support size) by St. Venant's principle. This observation allows to interpret the ringing artifact as a discretization error to be expected at the length scale of the support.

Appendix A. Discrete Fourier transformation
The 2D discrete Fourier transformation is defined as follows through out the paper. The generalization to 1D or 3D is straight forward. We divide the simulation domain of edge lengths L x , L y into N x , N y voxels of equal edge lengths ∆ i = L i /N i in each spatial direction i ∈ {x, y}. The lower left corner of voxel n x = I, n y = J is then given by For the discrete Fourier transform of the function f ( r) we use with the corresponding inverse transformation

Appendix B. Small-strain projection
It is often useful to carry out calculations in the small-strain limit. Small strains are special because the strain tensors (that replaces the deformation gradient) has to remain symmetric. For a formulation that involves multiple elements per voxel, we also require multiple symmetric strain tensors per voxel. While in the finite strain formulation, these were absorbed in our derivative indices α, β etc., we need to introduce a specific element index for the smallstrain case and can no longer distinguish between derivatives and Cartesian coordinates because of the symmetry in between derivatives and coordinates. Additionally to the previous introduced indices we will therefore use capital Greek letters (Θ, Λ, Ξ, . . . ) to denote elements. Components of the strain tensor will be denoted by small Latin indices.
We now introduce the strain tensor e in lieu of the deformation gradient, Eq. (2). The strain tensor for element Θ is where ∇ Θ is the (potentially discrete) derivative operator for element Θ and u( r) = χ( r) − r are the displacements from the undeformed positions r. For a given e, we minimize the residual R = Θ k R * Θ ( k) : R Θ ( k) with with respect to u * , where we have transformed into the Fourier-space and introduced the gradient operatorˆ D Θ ( k). The full residual is given by By combining the pairs of indices α = Θ, j and β = Λ, l, we can write this in the same form as the large strain projection,ê iα =Ĝ iαβjêjβ . Note that for a single element Θ = 1, we can write down the expression forĥ analytically, Using this expression and the Fourier derivative forˆ D( k) gives the small-strain projection operator of Ref. [40,Section 6].
Appendix C. Fourier-type projection operator on two evaluation points per voxel The Fourier-type derivative can be extended to several gradient evaluation points per voxel. This is useful to investigate the influence of Gibbs ringing separately from the effect of missing degrees of freedom. The standard Fouriertype derivative is evaluated at the grid points of the Fourier grid, which are the centers of the voxels. Hence, for a triangular mesh we additionally evaluate the Fourier derivative at the geometrical center of each triangle. For a two dimensional rectangular cell of edge lengths ∆ i , as shown in Fig. 2c, that means applying a shift of ±(∆ 1 , ∆ 2 )/6 from the center. The Fourier derivative operator D α ( k) = ik α acquires a phase to yield where we used explicit the indices 1 and 2 to denote the derivative operator in the center of the lower and upper triangle as shown in Fig. 2c.