<figure><img src="image_1.png" /></figure>

# A tetrahedron-based discretization for FFT-based computational homogenization with smooth solution fields

A. Finel ID

Université Paris-Saclay, ONERA, CNRS, Laboratoire d'Etude des Microstructures (LEM), 92322, Châtillon, France

## ARTICLE INFO

**Keywords:**

Homogenization

Fast Fourier transform (FFT)

Lippmann-Schwinger equation

Micromechanics

## ABSTRACT

We propose a new discrete FFT-based method for computational homogenization of micromechanics on a regular grid that is simple, fast and robust. The discretization scheme is based on a tetrahedral stencil that generates two independent mechanical problems on distinct face centred cubic subgrids. The resulting equilibrium problems may be merged into a single one that is solved by a gradient-based method equivalent to the original Moulinec-Suquet basic scheme.

The tetrahedral stencil displays three crucial properties. First, and most importantly, the Fourier representation of the associated Green operator is defined for any finite q-vector generated by the periodic boundary conditions and that does not belong to the reciprocal lattice of the discrete grids. As shown in the paper, this property guarantees that, for any elastic contrast, even infinite, mechanical equilibrium is always mathematically stable, i.e. free of any unphysical patterns, such as oscillations, ringing or checkerboarding, a property which is not shared by the original Moulinec-Suquet discretization (Moulinec and Suquet, 1994;1998) nor by the rotated scheme proposed by Willot (2015). Second, the components of tensorial quantities are all defined on the same location, which permits the use of any elastic anisotropy and any spatial variation of the material fields. Third, convergence to equilibrium using the simplest iterative scheme, i.e. the basic scheme, is fast and the number of iterates stabilizes at high contrast, so that infinite contrast is obtained without additional computational cost.

The efficiency and accuracy of the proposed tetrahedral stencil are assessed through several case studies, including microstructures that consist of a single inclusion with various elastic properties (finite eigenstrain, void, hard precipitate), a microstructure that contains simultaneously soft and stiff inclusions and complex 3D microstructures with random spherical voids. Convergence of effective elastic properties with grid refinement is investigated. An implementation of the tetrahedral stencil within the polarization-based fixed point algorithm of Monchiet-Bonnet is also proposed.

## 1. Introduction

Computational micromechanics and homogenization require the solution of the mechanical equilibrium of complex microstructures. Closed-form estimates of average properties, often based on the Eshelby inclusion problem, are of little use because real microstructures depart a lot from a single inclusion embedded in an infinite matrix. The alternative is to resort to full-field simulations of a complete and realistic microstructure. Among the available numerical techniques, Fourier-based methods have, during the last decades, gained an increasing popularity for several reasons. A full finite-element discretization of the microstructure,
which is a very complex operation in highly heterogeneous systems, is not required. The calculation is simply based on a regular grid which, in turn, may be directly identified with the voxels of an image obtained by 3D imaging techniques. Also, the underlying regular grid permits the use of very efficient libraries that implement Fast Fourier Transform (FFT) algorithms. Finally, the method is easy to implement and produces impressive results for large-scale microstructures in a very short time.

Fourier-based methods have been introduced by Moulinec-Suquet [1,2] in the context of computational homogenization (see [3–5] for recent reviews). They have also long been popular in the field of structural transformations since the pioneering work of Khachaturyan [6] on the elastic theory of solids and its extension to elastically inhomogeneous systems [7] (see an earlier review in [8]).

In fact, the original Moulinec-Suquet [1,2] and Khachaturyan [7] methods for dealing with elastically inhomogeneous solids are based on the same foundations, consisting of a continuous formulation of mechanical equilibrium in a periodic unit cell, which is then solved in Fourier space using a simple iterative algorithm. Indeed, as noted for example in [9], the so-called "basic scheme" introduced by Moulinec-Suquet may be interpreted as a Neumann expansion of the inhomogeneous elastic Green function associated to the microstructure. This is precisely the starting point of the perturbation method introduced by Khachaturyan, which can be recast as an iterative algorithm, as shown in [10], which in turn happens to be equivalent to the basic scheme.

It has rapidly been observed that the original method, in its most basic form, suffers from two difficulties. First, oscillations are often observed in the computed stress and strain fields in the presence of sharp material discontinuities. This is due to the numerical implementation of the method. Indeed, the mechanical equilibrium of the periodic unit cell and its Lippmann-Schwinger formulation are initially derived in a continuous setting. However, for the numerical implementation, a discretization of the real space must be considered. This discretization is associated, in the Fourier space, to a high frequency cut-off that generates oscillations and ringing patterns in real space. Second, convergence to equilibrium is often very slow. Indeed, as shown in [11] (see also [12,13]), when the elastic contrast $C$ is large but finite, convergence to equilibrium is linear but the number of iterates to reach equilibrium is proportional to $C$, which may lead to the conclusion that the basic scheme does not converge when $C$ is infinite, such as in porous materials. In fact, regarding this matter, a distinction must be made between the iterative scheme itself, to which we refer below as the basic scheme, and the discretization method used to approximate the continuum. It has been shown [14] that, in the continuous setting and under suitable conditions on the elastic energy density, the basic scheme for porous materials converges logarithmically for any distribution of the porous phase and that the convergence becomes linear under stronger assumptions on the distribution of voids (as example, convex voids). However, linear convergence may be lost when the continuum is discretized, which is the case with the Moulinec-Suquet discretization.

Accelerated algorithms have been proposed to overcome the convergence problems encountered with the basic scheme, such as the method proposed by Eyre and Milton [9] in the context of electrical conductivity and further extended to elasticity [11,15], the augmented Lagrangian method proposed by Michel et al. [11,16] and the polarization method proposed by Monchiet and Bonnet [13] which, under specific choices of parameters, reduces to each of the previous ones [12]. As shown in [12,17], these methods converge also linearly if the material contrast $C$ is finite but the convergence rate is proportional to the square root of $C$, which leads to a significant improvement over the basic scheme. However, the simplicity of the original method is lost and, with the Moulinec-Suquet discretization and infinite contrast, linear convergence is again lost and the number of iterates to reach equilibrium becomes exceedingly high (see for example [18]).

Another FFT methodology, based on the Hashin-Shtrikman variational principle, has been proposed by Brisard and Dormieux [19, 20]. The main ingredient of this method is to consider piecewise constant approximations of the polarization field and ask the corresponding discrete variables to fulfil an energy principle instead of solving a discrete Lippmann-Schwinger equation, as it is done within the original Moulinec-Suquet method. As the associated algorithm iterates on successive approximations of the polarization field, the method is also referred to as a polarization-based method.

In its spirit, this polarization-based method is similar to the method introduced earlier by Rodney et al. [21] for modelling dislocations within a phase field approach. In this latter method, dislocation loops are decomposed into elementary loops, called *loopons*, within which the associated phase field is piecewise constant. The elastic interactions between the loopons is then given by *decorated* interactions matrices whose Fourier transforms are given by infinite sums over the discrete sites of the reciprocal lattice in Fourier space. The resulting Green operator is exactly the so-called *periodized* Green operator of the polarization-based method proposed by Brisard and Dormieux [19,20]. Thanks to its (piecewise) continuous character, the polarization-based method is automatically free of oscillatory artefacts and can handle infinite stiffness contrasts. However, as noted in [20], the infinite series involved in the calculation of the periodized Green operator converge very slowly, making their numerical implementation rather involved, especially for 3D cases.

More recently, Eloh et al. [22] proposed a method which is also based on a piecewise constant approximation of local fields, but applied at the level of the Lippmann-Schwinger equation. The method therefore also leads to a redefinition of the Green operator in the form of infinite series in reciprocal space and, consequently, suffers from the same weakness as the polarization-based method of Brisard and Dormieux.

As alternatives to the original fixed-point basic scheme and its subsequent evolutions, numerous gradient-based iterative solvers have been proposed, such as Krylov Conjugate Gradient methods for linear problems [20,23,24], extensions of these Krylov methods to Newton-Krylov solvers for nonlinear situations [25,26], fast gradient methods [27–29] and Quasi-Newton methods [30].

All the FFT-based methods mentioned above have in common the formulation of the initial mechanical problem in a continuous setting, followed by a specific discretization strategy to generate a discrete Green operator. An alternative to this route is to readily start with a finite difference discretization of the relation between the strain and displacement fields, both defined on specific discrete
grids. Besides the centred discretization scheme, which is known to generate strong oscillations, various discretization methods have been proposed.

The standard staggered grid scheme, initially developed for fluid mechanics problems [31], has been used in a Fourier-based setting by Geslin and Finel [32] and Schneider et al. [33]. Its main advantage is that it is free from any spurious oscillations. It has also been observed that the staggered grid discretization associated with the basic scheme leads to a linear convergence in the presence of spherical pores [14]. However, the price to pay for these benefits is that different components of the displacement field and also of the strain field are stored on different locations, which generates difficulties in presence of sharp interfaces and in situations in which the constitutive law links together tensor components that do not sit on the same location. The discretization scheme proposed by Willot [34], based on a rotated staggered grid initially introduced for simulating the propagation of elastic waves [35], avoids these difficulties. As it presents several positive properties, in particular the possibility to simulate porous materials and the reduction of the oscillations observed in the Moulinec–Suquet discretization, the rotated scheme has become the standard discretization method for FFT-based numerical homogenization. However, as will be shown below, it is not exempt of spurious oscillations or checkerboard patterning effects. Finally, the discretization scheme proposed by Ruffini and Finel [36] avoids the weaknesses of the two previous schemes: the components of a given tensor are all stored on the same location and it is free from any non-physical instability. However, this comes at the expense of introducing a subgrid discretization in which, if $D$ is the spatial dimension, each initial cubic voxel, to which is assigned a single displacement vector, is subdivided into $2^D$ subvoxels that each receives a discrete strain tensor. Even though the method operates on the displacement field only, the $2^D$ subdivisions increases in the same proportion the computational load.

More recently, Schneider et al. [37] introduced another FFT-based method that relies on a finite-element discretization with trilinear hexahedral elements and which shares similarities with the previous method. The displacement field is discretized on the sites of a single cubic grid whereas, inside each cubic voxel, the strain field is discretized on $2^D$ Gauss points. This scheme is also devoid of any unphysical oscillations but, as for the previous method, this comes at the expense of the computational effort due to the doubly fine grid and the corresponding $2^D$ integration points per cubic voxel. This shortcoming was addressed by Leuschner and Fritzen [38], who derived a formulation of the method that operates on the displacement field only, but at the expense of a high memory demand associated with the assembly of the finite-element stiffness matrix, which penalizes the method when dealing with non-linear problems. This formulation was subsequently recast by Ladecký et al. [39] within a matrix-free iterative method. Interestingly, it has been shown in [37] that the rotated finite difference scheme proposed by Willot can also be interpreted as a finite-element discretization with trilinear hexahedral elements, but with only one Gauss point per cubic voxel instead of eight for full integration. This analogy shed a new light on the deficiencies of the rotated scheme. It is argued that the reduced integration is at the origin of the hourglass instabilities mentioned above, which subsequently led Schneider [40] to propose a method to control hourglassing.

In the present work, we propose a new discretization method, based on a tetrahedral stencil between the strain and displacement fields. This stencil avoids the difficulties encountered with the discrete schemes discussed above. It is mathematically free from any non-physical instabilities, does not resort on subvoxels defined on a finer scale and the components of a given tensor are all defined on the same location. Finally, convergence to mechanical equilibrium using the basic scheme is fast, even in the presence of voids.

The formulation of the tetrahedral discretization is based on the division of two initial staggered simple cubic grids, on which we define respectively the strain and displacement fields, into disjoint face-centred cubic subgrids. The procedure leads to two independent discretization schemes. We show that the two resulting mechanical equilibria can be merged into a single gradient-based iterative scheme which operates on the full initial displacement grid, and which can be interpreted as the adaptation of the Moulinec–Suquet basic scheme to the current discretization. This reduction is made possible by the introduction of fields, called $T_1$ and $T_2$ fields, defined on the full initial strain grid and which appropriately mix the two discretizations.

We also propose an adaptation of the polarization method of Monchiet–Bonnet that can be useful when the basic scheme is not efficient enough, which can be the case for example when considering a microstructure containing simultaneously soft and hard inclusions.

The paper is organized as follows. In Section 2, we develop the model, its implementation within a gradient-based method and we present an adaptation of the polarization method of Monchiet–Bonnet. In Section 3, we analyse the performance of the tetrahedral stencil method for different applications, including comparisons with the original Moulinec–Suquet method, the rotated finite difference scheme and the standard staggered grid method. Convergence of effective elastic properties under grid refinement is also investigated. We conclude in Section 4.

The following notation is used throughout the paper. $a$ refers to a scalar quantity whereas bold symbols $\underline{a}$, $\underline{a}$, $\underline{a}$ and $\underline{a}$ denote first-, second-, third- and fourth-order tensor quantities, respectively. $\bar{A}$ is the complex conjugate of $A$, $\hat{f}(q)$ represents the Fourier transform of $f(r)$, $\langle X \rangle$ is the spatial average of $X$ and, as usual, single dot (.), double-dot (:) and the symbol $\otimes_s$ represent contraction, double contraction and the symmetrized tensor product, respectively.

## 2. Method : a tetrahedral stencil

The problem at hand is to solve for static mechanical equilibrium in a 3D cell containing a representative volume of a microstructure composed of one or more phases embedded into a matrix and subjected to an applied overall strain or external stress. We use linear elasticity at small strains in its standard form as our point of departure for discussing the discretization schemes and the solution methods. However, the discretization method proposed below is general and may be used in more complex situations, such as nonlinear elasticity or plastic flow. The microstructure is described by a stiffness distribution $\lambda(r)$ and by an eigenstrain
distribution $\epsilon^0(r)$ that measures a free-stress deformation field with respect to the matrix. The local stress $\sigma(r)$ is therefore related to the local strain $\epsilon(r)$ by

$$ \sigma(r) = \lambda(r) : (\epsilon(r) - \epsilon^0(r)). \qquad (2.1) $$

The central point of our treatment is a proper choice of a discretization scheme that relates the strain and displacement fields. We proceed as follows.

## 2.1. Discretization scheme

We consider a 3D periodic rectangular cell discretized into $N_1 \times N_2 \times N_3 = N$ cubic voxels. For a reason that will be given later, the integers $N_i, i = 1, 2, 3$, are even. The edges of the cell are parallel to the cartesian coordinates and the voxel edge length is equal to $d$ (see Fig. 1). We refer to voxel $(l_1, l_2, l_3)$ as the voxel whose centre sits at position $\underline{r}$ given (in units of $d$) by

$$ \begin{gathered} \underline{r} = l_1 \underline{a}_1 + l_2 \underline{a}_2 + l_3 \underline{a}_3, \\ l_i \in \{0, \dots, N_i - 1\} \quad (i = 1, 2, 3), \end{gathered} \qquad (2.2) $$

where $\underline{a}_i, i = 1, 2, 3$, are unit vectors along the cartesian axes. The centres of the voxels form a simple cubic (SC) lattice with a grid spacing equal to $d$. The eight corners of cell $(l_1, l_2, l_3)$ sit on sites whose coordinates are $(l_1 \pm 1/2, l_2 \pm 1/2, l_3 \pm 1/2)$. Altogether, the corners of all cells form another SC lattice.

The first step is to choose the locations where the different tensors are defined. In order to handle any elastic anisotropy and any spatial variation of the material fields, the discrete version of the constitutive relation that link stresses and strains must involve tensorial quantities all defined on the same location. Here, all the 2nd-order and 4th-order tensors, i.e. the strain and stress fields as well as the material stiffness and eigenstrain, will be defined on the SC voxel grid, hereafter referred to as the SC strain grid, whereas the displacement field will live on the SC corner grid, hereafter noted as the SC displacement grid. In particular, all the components (diagonal and non-diagonal) of a tensor live on the same location, which avoids the difficulty that arises when the discretization of a constitutive laws relates a component of a tensor to components of another one defined on different locations, as it is the case in the usual SC staggered scheme mentioned above (see for example [33]).

The key point now is to select a judicious stencil that relates the strain field to the displacement field. In order to automatically avoid the instabilities often associated to finite difference schemes, we propose a tetrahedral stencil, in which the strain at the centre of a voxel is linked to the displacement vectors that sit at the fourth corners of one of the two regular tetrahedra that can be constructed out of the eight corners of that voxel (see Fig. 1). In other words, we split the SC displacement grid into two face centred cubic (FCC) subgrids $Y_1$ and $Y_2$, whose sites in Fig. 1 are identified with labels 1 and 2, respectively. Considering one of the FCC displacement subgrid, say subgrid $Y_1$, we observe that any site of the SC strain grid sits at the centre of one of the tetrahedra formed by the nodes of this subgrid. This offers the possibility to define two tetrahedral discretization schemes, labelled 1 and 2, associated to displacement subgrids $Y_1$ and $Y_2$, respectively (see Fig. 2).

The discretization scheme associated to one of the FCC displacement subgrids does not involve the fluctuating displacements on sites that belongs to the other subgrid. As a result, if we select only one of the two FCC discretization schemes, some symmetry elements on the SC strain grid that would be preserved if we had used a full SC displacement grid are now lost. In order to recover the full cubic symmetry on the SC strain grid, we consider simultaneously the two discretization schemes. The final strain field will be simply defined as the average of the strain fields generated by the two discretization schemes. However, since this averaging is “on site”, it will not generate a strain field less accurate than the two initial fields, as would be observed if some spatial averaging process had been used (see for example [41]). Also, although one might naively expect that using two discretization schemes would double the computational effort in terms of memory and CPU time, we will see below that this is not the case, as the two mechanical equilibria that emerge from the two discretization schemes can be collapsed into a single one.

Next, considering one of these FCC discretization schemes, say discretization scheme 1, we observe that the translations that relate the four vertices of a tetrahedron, which sit on subgrid $Y_1$, to its centre, which sits on the SC strain grid, alternate between two different sets, represented in red and blue, respectively, in Figs. 1 and 2. In other words, the SC strain grid must also be split into two FCC subgrids $Y_A$ and $Y_B$, whose sites, in Figs. 1 and 2, are identified with the labels A and B, respectively.

We are now in position to define the finite difference schemes that relate the strains defined on the FCC strain subgrids to the displacements defined on the FCC displacement subgrids. We first split the strains into their spatial averages and fluctuating terms, i.e. for discretization schemes $k = 1, 2$ and for strain subgrids $Y_K$, $K = A, B$, we write:

$$ \text{if } \underline{r} \in Y_K, \quad \epsilon^{(Kk)}(\underline{r}) = \bar{\epsilon} + \delta\epsilon^{(Kk)}(\underline{r}), \qquad (2.3) $$

where the notation $A^{(Kk)}(\underline{r})$ refers to a quantity that lives on a site $\underline{r}$ that belongs to strain subgrid $Y_K$ and that is defined using the discretization scheme associated to displacement subgrid $Y_k$. For strain subgrids $Y_K, K = A, B$ and discretization schemes $Y_k, k = 1, 2$, the fluctuating strain fields $\delta\epsilon^{(Kk)}(\underline{r})$, whose spatial averages are equal to zero, are linked to the fluctuating displacement field $\underline{u}(\underline{r})$ by:

$$ \text{if } \underline{r} \in Y_K, \quad \delta\epsilon_{ij}^{(Kk)}(\underline{r}) = \frac{1}{2} \{ D_i^{(Kk)}[u_j](\underline{r}) + D_j^{(Kk)}[u_i](\underline{r}) \}, \qquad (2.4) $$

where the dimensionless quantities $u_i, i = 1, 2, 3$, represent the coordinates of the displacement field in units of the cell length $d$. Following the lines of the tetrahedral stencil described above, inspection of Fig. 1 shows that the finite difference schemes that connect strain subgrid $Y_A$ to displacement subgrid $Y_1$ on the one hand, and strain subgrid $Y_B$ to displacement subgrid $Y_2$ on the
<figure><img src="image_2.png" /><figcaption>Fig. 1. Tetrahedral stencil for a simple cubic discretization. Only two voxels along axis (100) are shown. The simple cubic displacement grid, whose sites are located at the corners of the voxels, is split into two FCC subgrids whose sites are labelled by the indices 1 and 2. Similarly, the simple cubic strain grid, whose sites are located at the centres of the voxels, is split into two FCC subgrids whose sites are labelled by the indices A and B. The finite difference schemes that connect strain subgrid $Y_A$ to displacement subgrid $Y_1$ and strain subgrid $Y_B$ to displacement subgrid $Y_2$ are geometrically associated with tetrahedra (in blue) that, up to translations, are identical to a single one (tetrahedron $T_1$), which makes the finite difference operators $D^{(A1)}$ and $D^{(B2)}$ equal. Similarly, the finite difference schemes $D^{(A2)}$ and $D^{(B1)}$ are equal to each other, because they are associated to tetrahedra (in red) also identical to a single one (tetrahedron $T_2$). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2. FCC displacements subgrids 1 and 2 and their links to the FCC strain subgrids A and B. Left: within each initial cubic voxel, the four sites that belongs to the FCC displacement subgrid 1 are related to the strain site that sits on its centre either by a blue tetrahedron (strain site of type A) or a red one (strain site of type B). Right: conversely, again within each cubic voxel, the four sites that belongs to the FCC displacement subgrid 2 are related to the strain site that sits on its centre either by a red tetrahedron (strain site of type A) or a blue one (strain site of type B). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

other, are geometrically associated with tetrahedra (in blue in Fig. 1) that, up to translations, are identical to a single one, henceforth referred to as tetrahedron $T_1$. This implies that the finite difference operators $D^{(A1)}$ and $D^{(B2)}$ are equal. Similarly, the operators $D^{(A2)}$ and $D^{(B1)}$, also associated to tetrahedra (in red in Fig. 1) identical to a single one referred to as tetrahedron $T_2$, are equal. Next, we observe that the finite difference scheme associated to $T_2$ is nothing but the negative of the image of the finite difference scheme $T_1$ through a central symmetry around its centre. In all, this leads to, for $i = 1, 2, 3$,

$$D_i^{(A1)}[f](lmn) = D_i[f](lmn), \qquad (2.5a)$$

$$D_i^{(B2)}[f](lmn) = D_i[f](lmn), \qquad (2.5b)$$

$$D_i^{(A2)}[f](lmn) = -C_{lmn} \{ D_i[f](lmn) \}, \qquad (2.5c)$$

$$D_i^{(B1)}[f](lmn) = -C_{lmn} \{ D_i[f](lmn) \}, \qquad (2.5d)$$
where the operator $D_i$, $i = 1, 2, 3$, act on a discrete scalar field $f$ according to

$$ 
\begin{aligned}
D_1[f](lmn) &= \frac{1}{2} \left\{ f \left(lmn + \frac{111}{2}\right) + f \left(lmn + \frac{\bar{1}\bar{1}}{2}\right) - f \left(lmn + \frac{\bar{1}\bar{1}1}{2}\right) - f \left(lmn + \frac{\bar{1}\bar{1}\bar{1}}{2}\right) \right\}, \\
D_2[f](lmn) &= \frac{1}{2} \left\{ f \left(lmn + \frac{111}{2}\right) + f \left(lmn + \frac{\bar{1}\bar{1}\bar{1}}{2}\right) - f \left(lmn + \frac{\bar{1}\bar{1}1}{2}\right) - f \left(lmn + \frac{\bar{1}\bar{1}\bar{1}}{2}\right) \right\}, \\
D_3[f](lmn) &= \frac{1}{2} \left\{ f \left(lmn + \frac{111}{2}\right) + f \left(lmn + \frac{\bar{1}\bar{1}\bar{1}}{2}\right) - f \left(lmn + \frac{\bar{1}\bar{1}\bar{1}}{2}\right) - f \left(lmn + \frac{\bar{1}\bar{1}\bar{1}\bar{1}}{2}\right) \right\}.
\end{aligned} 
\quad (2.6) $$

The short hand notation $(lmn)$ refers to a site that belongs to the strain grid and whose coordinates are $(l, m, n)$, $(lmn + \frac{abc}{2})$ refers to a site that sits on the displacement grid and whose coordinates are $(l + \frac{a}{2}, m + \frac{b}{2}, n + \frac{c}{2})$ and $\bar{1}$ is equal to $-1$. The operator $C_{lmn}$ in Eqs. (2.5c) and (2.5d) represents the central symmetry around site $(lmn)$, i.e. $C_{lmn}(f(lmn + \frac{abc}{2})) = f(lmn - \frac{abc}{2})$. Identification of the strain components using (2.4) and the finite difference operators $D_i^{(Kk)}$ defined in (2.5) is illustrated in Appendix A. Finally, we impose periodic boundary conditions (PBC) on the microstructure, i.e. on the stiffness $\lambda_{\tilde{\omega}}(\mathbf{r})$ and eigenstrain $\epsilon^0(\mathbf{r})$ which are defined on the SC strain grid, and on the fluctuating displacement field $\underline{u}(\mathbf{r})$ which is defined on the SC displacement grid,

$$ \underline{r} \in Y_1 \cup Y_2, \quad \underline{u}(\underline{r} + N_1 \underline{a}_1) = \underline{u}(\underline{r} + N_2 \underline{a}_2) = \underline{u}(\underline{r} + N_3 \underline{a}_3) = \underline{u}(\underline{r}), \quad (2.7) $$

and similar conditions on $\tilde{\omega}(\underline{r})$ and $\tilde{\epsilon}^0(\underline{r})$. This automatically guarantees that the averages of the fluctuating strain fields $\delta\tilde{\epsilon}^{(Kk)}(\underline{r})$ defined in Eq. (2.4) are equal to zero. Also, the PBC are transmitted to $\delta\tilde{\epsilon}^{(Kk)}(\underline{r})$, and therefore to $\tilde{\epsilon}^{(Kk)}(\underline{r})$, through Eq. (2.4) and, consequently, to the stress field $\tilde{\sigma}^{(Kk)}(\underline{r})$ through the constitutive law, here given by Eq. (2.1). Note that within a given FCC subgrid, the coordinates of two consecutive sites along a cubic direction differ by two units (see Figs. 1 or 2). Therefore, the integers $N_i$, $i = 1, 2, 3$ must be even.

Next, we define the discrete Fourier transforms associated to our FCC subgrids.

## 2.2. Discrete fourier transforms

Fourier transforms on FCC displacement subgrids $Y_k$, $k = 1, 2$ and FCC strain subgrids $Y_K$, $K = A, B$, are defined by

$$ \hat{\underline{u}}^{(k)}(\underline{q}) = \frac{1}{N} \sum_{\underline{r} \in Y_k} \underline{u}(\underline{r}) \exp(-i\underline{q} \cdot \underline{r}), \quad (2.8a) $$

$$ \hat{\underline{A}}^{(Kk)}(\underline{q}) = \frac{1}{N} \sum_{\underline{r} \in Y_K} \underline{A}^{(Kk)}(\underline{r}) \exp(-i\underline{q} \cdot \underline{r}), \quad (2.8b) $$

where $\underline{A}^{(Kk)}(\underline{r})$ stands for any 2nd-order tensor field defined on subgrid $Y_K$ and associated to discretization scheme $k$. These Fourier transforms are defined for all q-vectors that account for the PBC on the SC grids defined in (2.7), i.e for the $N = N_1 \times N_2 \times N_3$ q-vectors defined (in units of $1/d$) by

$$ 
\begin{aligned}
\underline{q} &= q_1 \underline{a}_1^* + q_2 \underline{a}_2^* + q_3 \underline{a}_3^*, \\
q_i &= 2\pi \frac{h_i}{N_i} \quad \text{and} \quad h_i = 0, \dots, N_i - 1, \quad (i = 1, 2, 3),
\end{aligned} 
\quad (2.9) $$

where $(\underline{a}_1^*, \underline{a}_2^*, \underline{a}_3^*)$ is the reciprocal basis associated to $(\underline{a}_1, \underline{a}_2, \underline{a}_3)^*$.<fnref n="1" /> In the following, the short hand notation $(q_1 q_2 q_3)^*$ will refer to vector $\underline{q} = q_1 \underline{a}_1^* + q_2 \underline{a}_2^* + q_3 \underline{a}_3^*$.

The inverse Fourier transforms associated to the direct Fourier transforms (2.8a) and (2.8b) are given by

$$ \underline{r} \in Y_k, \quad \underline{u}(\underline{r}) = \sum_{\underline{q}} \hat{\underline{u}}^{(k)}(\underline{q}) \exp(i\underline{q} \cdot \underline{r}), \quad (2.10a) $$

$$ \underline{r} \in Y_K, \quad \underline{A}^{(Kk)}(\underline{r}) = \sum_{\underline{q}} \hat{\underline{A}}^{(Kk)}(\underline{q}) \exp(i\underline{q} \cdot \underline{r}), \quad (2.10b) $$

where the sums run on the q-vectors defined in (2.9). Note that these inverse Fourier transforms are defined only for the $N/2$ sites that belong to the relevant FCC subgrid. In other words, an inverse Fourier transform computed on a site $\underline{r}$ that does not belong to the relevant FCC subgrid is equal to zero:

$$ \underline{r} \notin Y_k : \quad \sum_{\underline{q}} \hat{\underline{u}}^{(k)}(\underline{q}) \exp(i\underline{q} \cdot \underline{r}) = 0, \quad (2.11a) $$

$$ \underline{r} \notin Y_K : \quad \sum_{\underline{q}} \hat{\underline{A}}^{(Kk)}(\underline{q}) \exp(i\underline{q} \cdot \underline{r}) = 0. \quad (2.11b) $$

We may now write the Fourier transform versions of the stencil (2.4). We obtain, for $k = 1, 2$ and $K = A, B$,

$$ \hat{\delta\tilde{\epsilon}}^{(Kk)}(\underline{q}) = \underline{D}^{(Kk)}(\underline{q}) \otimes_3 \hat{\underline{u}}^{(k)}(\underline{q}), \quad (2.12) $$

<footnote n="1">¹ Basis $(\underline{a}_1^*, \underline{a}_2^*, \underline{a}_3^*)$ and $(\underline{a}_1, \underline{a}_2, \underline{a}_3)$ are linked by $\underline{a}_j^* \cdot \underline{a}_j = \delta_{ij}$ where $\delta_{ij}$ is the Kronecker symbol. In the present situation, as $(\underline{a}_1, \underline{a}_2, \underline{a}_3)$ is orthonormal, we have $\underline{a}_i^* = \underline{a}_i$, $i = 1, 2, 3$.</footnote>
where the functions $\underline{D}^{(Kk)}(\underline{q})$ are the Fourier transforms of the discrete operators defined in (2.5). They are given by, for $i = 1, 2, 3$,

$$D_i^{(A1)}(\underline{q}) = D_i(\underline{q}), \qquad (2.13a)$$

$$D_i^{(B2)}(\underline{q}) = D_i(\underline{q}), \qquad (2.13b)$$

$$D_i^{(A2)}(\underline{q}) = -\overline{D_i(\underline{q})}, \qquad (2.13c)$$

$$D_i^{(B1)}(\underline{q}) = -\overline{D_i(\underline{q})}, \qquad (2.13d)$$

where the operators $D_i(\underline{q})$, $i = 1, 2, 3$, are the Fourier transforms of the operators given in (2.6),

$$D_1(\underline{q}) = \frac{1}{2} \left\{ \exp\left(\frac{i}{2}(q_1 + q_2 + q_3)\right) + \exp\left(\frac{i}{2}(q_1 - q_2 - q_3)\right) - \exp\left(\frac{i}{2}(-q_1 + q_2 - q_3)\right) - \exp\left(\frac{i}{2}(-q_1 - q_2 + q_3)\right) \right\}, \\ D_2(\underline{q}) = \frac{1}{2} \left\{ \exp\left(\frac{i}{2}(q_1 + q_2 + q_3)\right) - \exp\left(\frac{i}{2}(q_1 - q_2 - q_3)\right) + \exp\left(\frac{i}{2}(-q_1 + q_2 - q_3)\right) - \exp\left(\frac{i}{2}(-q_1 - q_2 + q_3)\right) \right\}, \qquad (2.14) \\ D_3(\underline{q}) = \frac{1}{2} \left\{ \exp\left(\frac{i}{2}(q_1 + q_2 + q_3)\right) - \exp\left(\frac{i}{2}(q_1 - q_2 - q_3)\right) - \exp\left(\frac{i}{2}(-q_1 + q_2 - q_3)\right) + \exp\left(\frac{i}{2}(-q_1 - q_2 + q_3)\right) \right\},$$

for all vectors $\underline{q} = q_1 \underline{a}_1^* + q_2 \underline{a}_2^* + q_3 \underline{a}_3^*$ defined in Eq. (2.9).

We note that, as expected, the operators $D_i(\underline{q})$, $i = 1, 2, 3$, are equal to zero for $\underline{q} = (000)^*$. This results from the fact that mode $\underline{q} = (000)^*$ corresponds to a uniform translation in real space and thus does not generate any strain. Therefore, the Fourier transform of any differential operator, whether discrete or continuous, vanishes at the origin of Fourier space.

We also observe that, in the present situation, the operators $D_i(\underline{q})$, $i = 1, 2, 3$, are equal to zero for $\underline{q} = (\pi\pi\pi)^*$, which implies that, within our tetrahedral stencil, mode $\underline{q} = (\pi\pi\pi)^*$ also cannot generate any strain. However, this is not problematic because this mode belongs to the reciprocal lattice of the FCC subgrids. Indeed, it is easily seen that we have, for $k = 1, 2$,

$$\text{if } \underline{r} \text{ and } \underline{r}' \in Y_k \text{ and } \underline{q} = (\pi\pi\pi)^*, \quad \exp(i\underline{q}\cdot\underline{r}) = \exp(i\underline{q}\cdot\underline{r}'). \qquad (2.15)$$

This implies that the contribution of the Fourier component $\hat{u}^{(k)}((\pi\pi\pi)^*)$ to the displacement field defined on any site of subgrid $Y_k$ does not depend on the site. In other words, the mode $\underline{q} = (\pi\pi\pi)^*$ corresponds also to uniform translations of the discrete displacements subgrids and, therefore, does not generate any strain.

In brief, the stencil used here has the important property that the Fourier transforms of the associated differential operators vanish only when they should, i.e. for q-modes that corresponds to uniform translations.

## 2.3. Mechanical equilibrium and a gradient-descent method

Formally, the discretization scheme $k$ is associated to the elastic energy (per unit volume)

$$E^{(k)} = \frac{1}{2N} \sum_{K=\Lambda,B} \sum_{\underline{r} \in Y_k} (\underline{\epsilon}^{(Kk)}(\underline{r}) - \underline{\epsilon}^0(\underline{r})) : \underline{\lambda}(\underline{r}) : (\underline{\epsilon}^{(Kk)}(\underline{r}) - \underline{\epsilon}^0(\underline{r})), \qquad (2.16)$$

which corresponds to the initial constitutive law considered here, i.e.

$$\underline{\sigma}^{(Kk)}(\underline{r}) = N \frac{\partial E^{(k)}}{\partial \underline{\epsilon}^{(Kk)}(\underline{r})} = \underline{\lambda}(\underline{r}) : (\underline{\epsilon}^{(Kk)}(\underline{r}) - \underline{\epsilon}^0(\underline{r})). \qquad (2.17)$$

Under the condition of an overall applied strain, the average strain $\underline{\epsilon}$ introduced in the decomposition (2.3) is simply equal to the applied strain $\underline{\epsilon}^a$. Within this boundary condition, mechanical equilibrium corresponds solely to minimization of elastic energy with respect to the fluctuating components of the displacement field. Boundary conditions corresponding to an applied stress will be discussed later. Until otherwise stated, we consider here that the system is subjected to an applied strain.

As mentioned above, we consider simultaneously the two discretization schemes $k = 1$ and $k = 2$. Therefore, we introduce the energy functional $E$ defined as the average of the individual elastic energies,

$$E = \frac{1}{2} (E^{(1)} + E^{(2)}). \qquad (2.18)$$

As $E^{(k)}$ depends only on the displacements $\underline{u}(\underline{r})$ attached to the sites $\underline{r}$ that belong to the displacement subgrid $Y_k$, minimization of the average energy $E$ with respect to the full displacement field $\underline{u}$ is equivalent to the simultaneous minimizations of the individual energies $E^{(1)}$ and $E^{(2)}$. Consequently, the two mechanical equilibria are simultaneously resolved by

$$\forall \underline{r} \in Y_1 \cup Y_2, \quad \frac{\partial E}{\partial \underline{u}(\underline{r})} = 0, \qquad (2.19)$$

or equivalently by

$$\forall \underline{q}, \quad \frac{\partial E}{\partial \underline{q}} = 0, \qquad (2.20)$$

where $\underline{u}(\underline{q})$ is the Fourier transform of the full displacement field, i.e.

$$\begin{aligned} \hat{\underline{u}}(\underline{q}) &= \frac{1}{N} \sum_{\underline{r} \in Y_1 \cup Y_2} \underline{u}(\underline{r}) \exp(-i\underline{q} \cdot \underline{r}), \\ &= \hat{\underline{u}}^{(1)}(\underline{q}) + \hat{\underline{u}}^{(2)}(\underline{q}), \end{aligned} \qquad (2.21)$$
and

$$ \underline{u}(\underline{r}) = \sum_{\underline{q}} \hat{u}(\underline{q}) \exp(i\underline{q} \cdot \underline{r}), \qquad (2.22) $$

which, due to property (2.11a), is compatible with (2.10a). A straightforward calculation (see Appendix C) shows that the derivative of energy $E^{(k)}$ with respect to Fourier component $\hat{u}^{(k)}(\underline{q})$ is given by

$$ \frac{\partial E^{(k)}}{\partial \hat{u}^{(k)}(\underline{q})} = \sum_{K=A,B} \hat{\sigma}^{(Kk)}(\underline{q}) \cdot \underline{D}^{(Kk)}(\underline{q}), \quad k = 1, 2, \qquad (2.23) $$

where $\hat{\sigma}^{(Kk)}(\underline{q})$ is the Fourier transform of the stress field $\sigma^{(Kk)}(\underline{r})$. Therefore, from

$$ \frac{\partial E}{\partial \hat{u}(\underline{q})} = \frac{1}{2} \left\{ \frac{\partial E^{(1)}}{\partial \hat{u}^{(1)}(\underline{q})} + \frac{\partial E^{(2)}}{\partial \hat{u}^{(2)}(\underline{q})} \right\}, \qquad (2.24) $$

one gets

$$ \frac{\partial E}{\partial \hat{u}(\underline{q})} = \frac{1}{2} \left\{ \overline{\hat{\sigma}^{(T_1)}(\underline{q}) \cdot \underline{D}(\underline{q})} - \overline{\hat{\sigma}^{(T_2)}(\underline{q}) \cdot \underline{D}(\underline{q})} \right\}, \qquad (2.25) $$

where the quantities $\hat{\sigma}^{(T_1)}(\underline{q})$ and $\hat{\sigma}^{(T_2)}(\underline{q})$ are defined by

$$ \begin{aligned} \hat{\sigma}^{(T_1)}(\underline{q}) &= \hat{\sigma}^{(A1)}(\underline{q}) + \hat{\sigma}^{(B2)}(\underline{q}), \\ \hat{\sigma}^{(T_2)}(\underline{q}) &= \hat{\sigma}^{(A2)}(\underline{q}) + \hat{\sigma}^{(B1)}(\underline{q}). \end{aligned} \qquad (2.26) $$

A central point here is the emergence of these two combinations, that we have labelled respectively with the superscripts $T_1$ and $T_2$ initially introduced in Section 2.1. This emergence is due to the equality of the operators $\underline{D}^{(A1)}(\underline{q})$ and $\underline{D}^{(B2)}(\underline{q})$ (respectively, $\underline{D}^{(A2)}(\underline{q})$ and $\underline{D}^{(B1)}(\underline{q})$) observed in Eqs. (2.13), which in turn is due to the fact that the corresponding finite difference schemes are based on the same tetrahedron, namely tetrahedron $T_1$ (respectively, tetrahedron $T_2$), as explained in Section 2.1. Note that each of the two combinations labelled with the superscripts $T_1$ and $T_2$ mixes the two discretizations $Y_1$ and $Y_2$ and, therefore, none of them refers to a specific discretization. Without risk of confusion, we will hereinafter refer to fields $T_1$ and $T_2$ as fields defined by combinations such as those presented in (2.26).

We now formulate a minimization scheme in the lines of the basic scheme. As first shown in [26] (see also [14,27]), the Moulinec-Suquet basic scheme can be interpreted as a gradient method with a fixed step size if the gradient of the elastic energy is defined according to a suitable scalar product. This gradient-based interpretation can be adapted to the present discretization scheme in a way that leads to a single iterative scheme that incorporates simultaneously the two discretization schemes associated to the FCC displacement subgrids $Y_1$ and $Y_2$. This is done as follows.

First, in a similar way to the introduction of the stress fields $\hat{\sigma}^{(T_1)}(\underline{q})$ and $\hat{\sigma}^{(T_2)}(\underline{q})$, we introduce the q-space fluctuating strain fields

$$ \begin{aligned} \hat{\delta}_{\underline{q}}^{(T_1)}(\underline{q}) &= \hat{\delta}_{\underline{q}}^{(A1)}(\underline{q}) + \hat{\delta}_{\underline{q}}^{(B2)}(\underline{q}), \\ \hat{\delta}_{\underline{q}}^{(T_2)}(\underline{q}) &= \hat{\delta}_{\underline{q}}^{(A2)}(\underline{q}) + \hat{\delta}_{\underline{q}}^{(B1)}(\underline{q}). \end{aligned} \qquad (2.27) $$

Using Eq. (2.12) that relates the subgrid strains to the subgrid displacement fields, together with Eqs. (2.13) that gives the differential operators and the relation $\hat{u}(\underline{q}) = \hat{u}^{(1)}(\underline{q}) + \hat{u}^{(2)}(\underline{q})$, one gets

$$ \begin{aligned} \hat{\delta}_{\underline{q}}^{(T_1)}(\underline{q}) &= \underline{D}(\underline{q}) \otimes_s \hat{u}(\underline{q}), \\ \hat{\delta}_{\underline{q}}^{(T_2)}(\underline{q}) &= -\underline{D}(\underline{q}) \otimes_s \hat{u}(\underline{q}), \end{aligned} \qquad (2.28) $$

where we again observe that each of the fields $T_1$ and $T_2$ is linked to a specific finite difference operator. Also, as a $T_1$ or $T_2$ field is the combination of fields defined on the two subgrids $Y_A$ and $Y_B$, its inverse Fourier transform is defined on each site of the full SC strain grid, i.e., $\underline{X}^{(T_1)}$ being one of the fields that appear in (2.26) and (2.27), we have

$$ \underline{r} \in Y_A \cup Y_B, \quad \underline{X}^{(T_1)}(\underline{r}) = \sum_{\underline{q}} \hat{\underline{X}}^{(T_1)}(\underline{q}) \exp(i\underline{q} \cdot \underline{r}), \qquad (2.29) $$

and

$$ \hat{\underline{X}}^{(T_1)}(\underline{q}) = \frac{1}{N} \sum_{\underline{r} \in Y_A \cup Y_B} \underline{X}^{(T_1)}(\underline{r}) \exp(-i\underline{q} \cdot \underline{r}), \qquad (2.30) $$

which, due to property (2.11b), is compatible with the initial definitions of $\hat{\underline{X}}^{(T_1)}(\underline{q})$ given in (2.26) and (2.27).

Next, we define the space $H^1_{\text{per}}$ of periodic displacements fields with vanishing Fourier components for the modes $\underline{q} = (000)^*$ and $\underline{q} = (\pi\pi)^*$. The motivation for excluding these two modes is that, as explained at the end of the previous section, they both correspond to uniform translations of the discretization subgrids $Y_1$ and $Y_2$ on which the tetrahedral stencil is based. Therefore, within this stencil, these modes do not contribute to the elastic energy and can be excluded from the minimization scheme with
no prejudice. Let $\hat{\lambda}_\sim^0$ denotes a uniform elastic stiffness tensor, i.e. a positive definite fourth-order tensor with major and minor symmetries. We endow the space $H_{per}^1$ with an Hilbertian structure through the scalar product defined by

$$ 
\begin{aligned}
(\underline{u}, \underline{v})_{H_{per}^1} &= \frac{1}{2N} \sum_{\underline{r} \in Y_1 \cup Y_2} \left( \delta_{\underline{r}}^{(T_1)}(\underline{r}) : \hat{\lambda}_\sim^0 : \delta_{\underline{r}}^{(T_1)}(\underline{r}) + \delta_{\underline{r}}^{(T_2)}(\underline{r}) : \hat{\lambda}_\sim^0 : \delta_{\underline{r}}^{(T_2)}(\underline{r}) \right), \\
&= \frac{1}{2} \sum_{\underline{q}} \left( \delta_{\underline{q}}^{(T_1)}(\underline{q}) : \hat{\lambda}_\sim^0 : \delta_{\underline{q}}^{(T_1)}(\underline{q}) + \delta_{\underline{q}}^{(T_2)}(\underline{q}) : \hat{\lambda}_\sim^0 : \delta_{\underline{q}}^{(T_2)}(\underline{q}) \right),
\end{aligned}
\quad (2.31) $$

where, in q-space, the fluctuating strain fields $\hat{\delta}_{\underline{q}}^{(T_1)}$, $\hat{\delta}_{\underline{q}}^{(T_2)}$, $\hat{\delta}_{\underline{q}}^{(T_1)}$ and $\hat{\delta}_{\underline{q}}^{(T_2)}$ are related to the fluctuating displacements $\hat{u}$ and $\hat{v}$ through (2.28), i.e.

$$ 
\left\{ 
\begin{aligned}
\hat{\delta}_{\underline{q}}^{(T_1)}(\underline{q}) &= \underline{D}(\underline{q}) \otimes_s \hat{u}(\underline{q}), \\
\hat{\delta}_{\underline{q}}^{(T_2)}(\underline{q}) &= -\underline{D}(\underline{q}) \otimes_s \hat{u}(\underline{q}),
\end{aligned}
\right. 
\quad 
\left\{ 
\begin{aligned}
\hat{\delta}_{\underline{q}}^{(T_1)}(\underline{q}) &= \underline{D}(\underline{q}) \otimes_s \hat{v}(\underline{q}), \\
\hat{\delta}_{\underline{q}}^{(T_2)}(\underline{q}) &= -\underline{D}(\underline{q}) \otimes_s \hat{v}(\underline{q}).
\end{aligned}
\right. 
\quad (2.32) $$

We easily see that the form (2.31) is indeed symmetric and, thanks to the exclusion of the modes $\underline{q} = (000)^*$ and $\underline{q} = (\pi\pi\pi)^*$, positive definite on $H_{per}^1$. Inserting relations (2.32) into (2.31) leads to

$$ (\underline{u}, \underline{v})_{H_{per}^1} = \sum_{\underline{q}} \hat{u}(\underline{q}) \cdot \hat{G}^0(\underline{q})^{-1} \cdot \hat{v}(\underline{q}), \quad (2.33) $$

where $\hat{G}^0(\underline{q})^{-1}$, the inverse of the displacement Green operator, is given by

$$ \hat{G}^0(\underline{q})^{-1} = \frac{1}{2} \{ \underline{D}(\underline{q}) \cdot \hat{\lambda}_\sim^0 \cdot \underline{D}(\underline{q}) + \underline{D}(\underline{q}) \cdot \hat{\lambda}_\sim^0 \cdot \underline{D}(\underline{q}) \}. \quad (2.34) $$

According to the Riesz representation theorem, there exists an isomorphic map $R$ between $H_{per}^1$ and its dual space $H_{per}^{-1}$ such that, for any linear form $\underline{d}$ in $H_{per}^{-1}$,

$$ (\underline{d}, \underline{v})_{H_{per}^{-1} \times H_{per}^{1}} = (R \underline{d}, \underline{v})_{H_{per}^{1}}, \quad \forall \underline{v} \in H_{per}^{1}, \quad (2.35) $$

where the duality product is given by

$$ (\underline{d}, \underline{v})_{H_{per}^{-1} \times H_{per}^{1}} \equiv \frac{1}{N} \sum_{\underline{r} \in Y_1 \cup Y_2} \underline{d}(\underline{r}) \cdot \underline{v}(\underline{r}) = \sum_{\underline{q}} \hat{\underline{d}}(\underline{q}) \cdot \hat{\underline{v}}(\underline{q}). \quad (2.36) $$

With the scalar product (2.33), equality (2.35) leads to

$$ \hat{G}^0(\underline{q})^{-1} \cdot \hat{R} \underline{d}(\underline{q}) = \underline{d}(\underline{q}). \quad (2.37) $$

where we used the fact that $\hat{G}^0(\underline{q})^{-1}$ is symmetric. The operator $\underline{D}(\underline{q})$ is equal to zero for $\underline{q} = (000)^*$ and $\underline{q} = (\pi\pi\pi)^*$, which makes the inverse of the Green functions also equal to zero, leaving the Green function $\hat{G}^0(\underline{q})$ undefined. However, these modes have been excluded in the definition of the space $H_{per}^1$. Hence, we can safely set $\hat{G}^0(\underline{q}) = \underline{\alpha}$ for $\underline{q} = (000)^*$ and $\underline{q} = (\pi\pi\pi)^*$ and write (2.37) in the form

$$ \forall \underline{q}, \quad \hat{R} \underline{d}(\underline{q}) = \hat{G}^0(\underline{q}) \cdot \underline{d}(\underline{q}), \quad (2.38) $$

which shows that the Riesz map $R$ is well defined and equal to the Green operator $\hat{G}^0$. Note that, as seen in (2.34), $\hat{G}^0$ fulfils the following symmetry properties:

$$ \hat{G}_{ij}^0(\underline{q}) = \overline{\hat{G}_{ij}^0(\underline{q})} = \hat{G}_{ji}^0(\underline{q}). \quad (2.39) $$

The Hilbertian setting and associated scalar product defined above permit to define a gradient of the elastic energy, that lives in $H_{per}^1$, through

$$ \forall \underline{v} \in H_{per}^1, \quad (\underline{v} \cdot \underline{E}[\underline{u}], \underline{v})_{H_{per}^1} = \underline{E}[\underline{u}](\underline{v}), \quad (2.40) $$

where $\underline{E}[\underline{u}](\underline{v})$ is the directional derivative of the elastic energy $E$ at $\underline{u}$ in the direction of the field $\underline{v}$, i.e.

$$ 
\begin{aligned}
\underline{E}[\underline{u}](\underline{v}) &\equiv \sum_{\underline{r} \in Y_1 \cup Y_2} \frac{\partial \underline{E}}{\partial \underline{u}(\underline{r})} \cdot \underline{v}(\underline{r}), \\
&= \left( N \frac{\partial \underline{E}}{\partial \underline{u}}, \underline{v} \right)_{H_{per}^{-1} \times H_{per}^{1}},
\end{aligned}
\quad (2.41) $$

where $N \partial \underline{E}/\partial \underline{u}$ is regarded as a linear form that lives in $H_{per}^{-1}$. According to the Riesz map identified above, one gets

$$ 
\begin{aligned}
\underline{\nabla} \underline{E}[\underline{u}](\underline{q}) &= \hat{G}^0(\underline{q}) \cdot N \frac{\partial \underline{E}}{\partial \underline{u}}(\underline{q}), \\
&= \hat{G}^0(\underline{q}) \cdot \frac{\partial \underline{E}}{\partial \hat{u}(\underline{q})},
\end{aligned}
\quad (2.42) $$
and, with (2.25),

$$ \underline{\nabla E[\underline{u}]}(\underline{q}) = \hat{\mathcal{G}}^0(\underline{q}). \frac{1}{2} \{ \hat{\mathcal{G}}^{(T_1)}(\underline{q}).\underline{\mathcal{D}}(\underline{q}) - \hat{\mathcal{G}}^{(T_2)}(\underline{q}).\underline{\mathcal{D}}(\underline{q}) \}. \quad (2.43) $$

Finally, with a fixed step size equal to one, the gradient descent method generated by $\underline{\nabla E[\underline{u}]}$ reads, in q-space,

$$ \begin{aligned} \forall \underline{q}. \quad \hat{\underline{u}}_{(n+1)}(\underline{q}) &= \hat{\underline{u}}_{(n)}(\underline{q}) - \underline{\nabla E[\underline{u}]}_{(n)}(\underline{q}), \\ &= \hat{\underline{u}}_{(n)}(\underline{q}) - \hat{\mathcal{G}}^0(\underline{q}). \frac{1}{2} \{ \hat{\mathcal{G}}^{(T_1)}_{(n)}(\underline{q}).\underline{\mathcal{D}}(\underline{q}) - \hat{\mathcal{G}}^{(T_2)}_{(n)}(\underline{q}).\underline{\mathcal{D}}(\underline{q}) \}, \end{aligned} \quad (2.44) $$

where we recognize the structure of the Moulinec–Suquet basic scheme, here adapted to the tetrahedral stencil. We note that this equation explicitly uses, in the last term on the right-hand side, a quantity which, measuring in q-space the deviation from mechanical equilibrium (see (2.20) and (2.25)), is needed to test the convergence of the iterative scheme.

A comment on the Green operator $\hat{\mathcal{G}}^0$ is here in order. We can easily show that, for any vector $\underline{q}_1$ and if $\underline{q}_1 + (\pi\pi\pi)^*$ belongs to the set (2.9), we have

$$ \hat{\mathcal{G}}^0(\underline{q}_1 + (\pi\pi\pi)^*) = \hat{\mathcal{G}}^0(\underline{q}_1), \quad (2.45) $$

which means that $\hat{\mathcal{G}}^0(\underline{q})$ is invariant through the translation $(\pi\pi\pi)^*$. This reinforces the equivalence, within the tetrahedral stencil, between the modes $(000)^*$ and $(\pi\pi\pi)^*$ in that, beyond the fact that they both correspond to uniform translations of the subgrids $Y_1$ and $Y_2$, they confer to the displacement Green function (and to the associated strain-based Green functions defined below in Section 2.6) the same analytic behaviour in their respective neighbourhoods. In brief, mode $\underline{q} = (\pi\pi\pi)^*$ is strictly equivalent to mode $\underline{q} = (000)^*$, which is the only singular point that one should expect for a well-behaved Green operator. This is not the case with the rotated stencil proposed by Willot [34], in which the cancellation of the differential operator $\underline{\mathcal{D}}(\underline{q})$ at finite q-points generates singularities on the associated displacement and strain Green functions that are not equivalent to those present at $\underline{q} = (000)^*$. The resulting non-analytical behaviour generates, in real space, slowly decaying unphysical oscillations that cannot be suppressed by any recipe for assigning an arbitrary value to the Green functions when they are not defined, because, as shown in Appendix D, analytic continuation to these problematic q-points is impossible. In short, the tetrahedral stencil fulfils the important property that the Fourier representation of its Green's operator is defined and analytic everywhere in q-space except at the origin $(000)^*$ and at a q-point which is equivalent to the origin, i.e. at a q-point that belongs to the reciprocal lattice of the subgrids on which the stencil is defined. This property is not shared by the rotated stencil, which inevitably generates unphysical instabilities in real space, such as oscillations and checkerboard effects on the strain and stress fields. In contrast, the tetrahedral scheme proposed here, which guarantees that the operators $D_i(\underline{q})$ do not vanish for any q-vector that does not lie on the reciprocal lattice of the real-space subgrids, is automatically free from these artefacts.

We now comment further on the fields $T_1$ and $T_2$ introduced in (2.26) and (2.27). Let $\hat{\mathcal{X}}^{(T_i)}$ be one of the 2nd-order tensors defined in these equations. According to these definitions and using (2.10b) and (2.11b), the inverse Fourier transforms (2.29) of $\hat{\mathcal{X}}^{(T_i)}(\underline{q})$, $i = 1, 2$, fulfil

$$ \begin{aligned} \underline{\mathcal{X}}^{(T_1)}(\underline{r}) &= \underline{\mathcal{X}}^{(A1)}(\underline{r}) \quad \text{if } \underline{r} \in Y_A, \qquad \underline{\mathcal{X}}^{(T_1)}(\underline{r}) = \underline{\mathcal{X}}^{(B2)}(\underline{r}) \quad \text{if } \underline{r} \in Y_B, \\ \underline{\mathcal{X}}^{(T_2)}(\underline{r}) &= \underline{\mathcal{X}}^{(A2)}(\underline{r}) \quad \text{if } \underline{r} \in Y_A, \qquad \underline{\mathcal{X}}^{(T_2)}(\underline{r}) = \underline{\mathcal{X}}^{(B1)}(\underline{r}) \quad \text{if } \underline{r} \in Y_B. \end{aligned} \quad (2.46) $$

Hence, a stress or strain field $T_1$ or $T_2$ is, in real space, defined on each site of the SC strain grid. Also, according to (2.28), the strain fields $T_1$ and $T_2$ are directly linked to the displacement field $\underline{u}(\underline{r})$ defined on the full SC displacement grid. Therefore, our tetrahedral stencil can be interpreted without recourse to the FCC displacements subgrids $Y_1$, $Y_2$ and FCC strain subgrids $Y_A$, $Y_B$, even if its mathematical formulation must go through these lattice subdivisions, notably to show that, thanks to its tetrahedral character, the stencil is mathematically stable, i.e. automatically free from any oscillation or checkerboard effect that would be observed if the associated q-space operators were to vanish on any q-point not belonging to the relevant reciprocal lattice.

In summary, even-though it relies initially on two independent FCC discretization, the final formulation of the tetrahedral stencil is condensed into a single iterative algorithm that involves quantities defined on full SC grids. The underlying FCC character of the stencil is reflected in the appearance of two strain and stress fields, labelled above with the indices $T_1$ and $T_2$. In case of a linear constitutive law, these fields are of course linked by the relations

$$ \begin{aligned} \underline{\sigma}^{(T_1)}(\underline{r}) &= \underline{\lambda}_{\underline{\sigma}}(\underline{r}) : (\underline{\varepsilon}^{(T_1)}(\underline{r}) - \underline{\varepsilon}^0(\underline{r})), \\ \underline{\sigma}^{(T_2)}(\underline{r}) &= \underline{\lambda}_{\underline{\sigma}}(\underline{r}) : (\underline{\varepsilon}^{(T_2)}(\underline{r}) - \underline{\varepsilon}^0(\underline{r})), \end{aligned} \quad (2.47) $$

which are inherited from the constitutive law (2.17) through relations (2.46). Also, note that $\underline{\varepsilon}^{(T_1)}(\underline{r})$ and $\underline{\varepsilon}^{(T_2)}(\underline{r})$ are related to their corresponding fluctuating strain fields through the same decomposition as in (2.3):

$$ \begin{aligned} \underline{\varepsilon}^{(T_1)}(\underline{r}) &= \tilde{\underline{\varepsilon}} + \delta\underline{\varepsilon}^{(T_1)}(\underline{r}), \\ \underline{\varepsilon}^{(T_2)}(\underline{r}) &= \tilde{\underline{\varepsilon}} + \delta\underline{\varepsilon}^{(T_2)}(\underline{r}). \end{aligned} \quad (2.48) $$

Once the equilibrium has been solved, the final stress and strain fields are simply defined as the arithmetic averages of the corresponding fields,

$$ \underline{\sigma}(\underline{r}) = \frac{1}{2} \{ \underline{\sigma}^{(T_1)}(\underline{r}) + \underline{\sigma}^{(T_2)}(\underline{r}) \}, \quad (2.49a) $$

$$ \underline{\varepsilon}(\underline{r}) = \frac{1}{2} \{ \underline{\varepsilon}^{(T_1)}(\underline{r}) + \underline{\varepsilon}^{(T_2)}(\underline{r}) \}, \quad (2.49b) $$
which, again in case of a linear constitutive law, are linked by

$$ \underline{\sigma}(\underline{r}) = \underline{\lambda}(\underline{r}) : (\underline{\epsilon}(\underline{r}) - \underline{\epsilon}^0(\underline{r})). \qquad (2.50) $$

Finally, we note that, since the definitions (2.49) only involve quantities defined on the same site $\underline{r}$, the averaged fields $\underline{\sigma}(\underline{r})$ and $\underline{\epsilon}(\underline{r})$ are no less accurate than the initial ones, as might be observed if some spatial averaging process had been used (see for example [41]).

## 2.4. Iterative loop on the average strain when the system is subjected to an applied stress

The boundary conditions considered above correspond to an applied strain, which is identified with the average strain $\bar{\epsilon}$. We consider here the opposite simple situation of an applied external stress. More complex boundary conditions could also be addressed, such as mixed stress and strain loading [42,43], but are not considered here.

When the boundary conditions consists in applying an external stress, the elastic energy given in Eq. (2.18) must be supplemented by an external contribution,

$$ E = \frac{1}{2} (E^{(1)} + E^{(2)}) - \underline{\sigma}^a : \bar{\underline{\epsilon}}, \qquad (2.51) $$

where $\underline{\sigma}^a$ is the applied stress. Minimization of this energy with respect to $\bar{\underline{\epsilon}}$ leads to

$$ \langle \underline{\sigma}(\underline{r}) \rangle = \underline{\sigma}^a, \qquad (2.52) $$

where $\underline{\sigma}(\underline{r})$ is the local stress defined in Eq. (2.49a). This equation may be easily recast in a form that allows it to be solved using an iterative algorithm. This can be done as follows. Using successively Eqs. (2.50), (2.49b) and (2.48), Eq. (2.52) leads to

$$ \langle \underline{\lambda}(\underline{r}) \rangle : \bar{\underline{\epsilon}} = \underline{\sigma}^a - \langle \underline{\lambda}(\underline{r}) : \underline{\delta\epsilon}(\underline{r}) \rangle + \langle \underline{\lambda}(\underline{r}) : \underline{\epsilon}^0(\underline{r}) \rangle, \qquad (2.53) $$

where $\underline{\delta\epsilon}(\underline{r}) = \frac{1}{2}(\underline{\delta\epsilon}^{(T_1)}(\underline{r}) + \underline{\delta\epsilon}^{(T_2)}(\underline{r}))$. Shifting the left member of this equation to the right and adding the product $\underline{\lambda}^0 : \bar{\underline{\epsilon}}$ to each member gives

$$ \underline{\lambda}^0 : \bar{\underline{\epsilon}} = \underline{\sigma}^a + \langle \underline{\lambda}^0 - \langle \underline{\lambda}(\underline{r}) \rangle \rangle : \bar{\underline{\epsilon}} - \langle \underline{\lambda}(\underline{r}) : \underline{\delta\epsilon}(\underline{r}) \rangle + \langle \underline{\lambda}(\underline{r}) : \underline{\epsilon}^0(\underline{r}) \rangle, \qquad (2.54) $$

which suggests the following iterative scheme

$$ \bar{\underline{\epsilon}}_{(n+1)} = (\underline{\lambda}^0)^{-1} : \{ \underline{\sigma}^a + \langle \underline{\lambda}^0 - \langle \underline{\lambda}(\underline{r}) \rangle \rangle : \bar{\underline{\epsilon}}_{(n)} - \langle \underline{\lambda}(\underline{r}) : \underline{\delta\epsilon}_{(n)}(\underline{r}) \rangle + \langle \underline{\lambda}(\underline{r}) : \underline{\epsilon}^0(\underline{r}) \rangle \}, \qquad (2.55) $$

where $\bar{\underline{\epsilon}}_{(n)}$ and $\underline{\delta\epsilon}_{(n)}(\underline{r})$ are the average and fluctuating strain fields at iteration $n$.

## 2.5. Summary: the tetrahedral stencil within the displacement-based basic scheme

To highlight the simplicity of the overall method, we show in Algorithm 1 the algorithm which implements the displacement-based iterative procedure given in Eq. (2.44). The operator $\underline{D}(\underline{q})$ that appears in steps 1 and 5 is given in Eq. (2.14) and the Green function $\hat{\underline{G}}^0(\underline{q})$ is defined, through its inverse, in Eq. (2.34). The symbols FT and FT$^{-1}$ denote Fourier transform and Inverse Fourier transform, respectively. The iterative scheme is stopped when a convergence criterion, such as one of those discussed in Section 2.8, decreases below a prescribed tolerance.

## 2.6. The tetrahedral stencil and the strain-based basic scheme

For the sake of completeness, we translate now the previous algorithm, that iterates on the displacement field, into an algorithm that iterates on strains. The resulting scheme will be of course strictly equivalent, at each iteration step, to the initial one. From the iterative Eq. (2.44) and relations (2.28), which relate the fluctuating deformation fields $\hat{\underline{\delta\epsilon}}^{(T_1)}(\underline{q})$ and $\hat{\underline{\delta\epsilon}}^{(T_2)}(\underline{q})$ to the displacement field $\underline{\hat{u}}(\underline{q})$, we obtain, after straightforward algebra,

$$ \begin{aligned} \hat{\underline{\delta\epsilon}}^{(T_1)}(\underline{q}) &= \hat{\underline{\delta\epsilon}}^{(T_1)}(\underline{q}) - \frac{1}{2} \{ \hat{\underline{\Gamma}}^{0,d}_{\underline{\epsilon}}(\underline{q}) : \hat{\underline{\sigma}}^{(T_1)}(\underline{q}) + \hat{\underline{\Gamma}}^{0,nd}_{\underline{\epsilon}}(\underline{q}) : \hat{\underline{\sigma}}^{(T_2)}(\underline{q}) \}, \\ \hat{\underline{\delta\epsilon}}^{(T_2)}(\underline{q}) &= \hat{\underline{\delta\epsilon}}^{(T_2)}(\underline{q}) - \frac{1}{2} \{ \hat{\underline{\Gamma}}^{0,nd}_{\underline{\epsilon}}(\underline{q}) : \hat{\underline{\sigma}}^{(T_1)}(\underline{q}) + \hat{\underline{\Gamma}}^{0,d}_{\underline{\epsilon}}(\underline{q}) : \hat{\underline{\sigma}}^{(T_2)}(\underline{q}) \}. \end{aligned} \qquad (2.56) $$

The diagonal and non-diagonal Green operators $\hat{\underline{\Gamma}}^{0,d}_{\underline{\epsilon}}(\underline{q})$ and $\hat{\underline{\Gamma}}^{0,nd}_{\underline{\epsilon}}(\underline{q})$ are defined by

$$ \begin{aligned} \hat{\underline{\Gamma}}^{0,d}_{ijkl}(\underline{q}) &= \frac{1}{4} \{ D_l(\underline{q}) \hat{\underline{G}}^0_{jk}(\underline{q}) D_l(\underline{q}) + D_j(\underline{q}) \hat{\underline{G}}^0_{lk}(\underline{q}) D_l(\underline{q}) + D_l(\underline{q}) \hat{\underline{G}}^0_{jl}(\underline{q}) D_k(\underline{q}) + D_j(\underline{q}) \hat{\underline{G}}^0_{il}(\underline{q}) D_k(\underline{q}) \}, \\ \hat{\underline{\Gamma}}^{0,nd}_{ijkl}(\underline{q}) &= -\frac{1}{4} \{ D_l(\underline{q}) \hat{\underline{G}}^0_{jk}(\underline{q}) D_l(\underline{q}) + D_j(\underline{q}) \hat{\underline{G}}^0_{lk}(\underline{q}) D_l(\underline{q}) + D_l(\underline{q}) \hat{\underline{G}}^0_{jl}(\underline{q}) D_k(\underline{q}) + D_j(\underline{q}) \hat{\underline{G}}^0_{il}(\underline{q}) D_k(\underline{q}) \}, \end{aligned} $$

or, in tensor notation,

$$ \begin{aligned} \hat{\underline{\Gamma}}^{0,d}_{\underline{\epsilon}}(\underline{q}) &= \underline{D}(\underline{q}) \otimes_s \hat{\underline{G}}^0(\underline{q}) \otimes_s \underline{D}(\underline{q}), \\ \hat{\underline{\Gamma}}^{0,nd}_{\underline{\epsilon}}(\underline{q}) &= -\underline{D}(\underline{q}) \otimes_s \hat{\underline{G}}^0(\underline{q}) \otimes_s \underline{D}(\underline{q}). \end{aligned} \qquad (2.57) $$
<caption>**Algorithm 1**: displacement-based basic scheme for the tetrahedral stencil on a simple cubic grid</caption>

initial guess for $\hat{u}(q)$ (e.g. $\hat{u}(q) = 0$)
initial guess for $\bar{\epsilon}$ (e.g. $\bar{\epsilon} = (\lambda^0)^{-1} : \sigma^a + \langle \epsilon^0(r) \rangle$ if applied stress, $\bar{\epsilon} = \bar{\epsilon}^a$ if applied strain)

$$1: \quad \hat{\delta}\epsilon^{(T_1)}(\underline{q}) = \underline{D}(q) \otimes_s \hat{u}(q) \\ \hat{\delta}\epsilon^{(T_2)}(\underline{q}) = -\overline{\underline{D}(q)} \otimes_s \hat{u}(\underline{q})$$

$$2: \quad \delta\epsilon^{(T_1)}(\underline{r}) = \text{FT}^{-1}\{\hat{\delta}\epsilon^{(T_1)}(\underline{q})\} \\ \delta\epsilon^{(T_2)}(\underline{r}) = \text{FT}^{-1}\{\hat{\delta}\epsilon^{(T_2)}(\underline{q})\}$$

$$3: \quad \sigma^{(T_1)}(\underline{r}) = \lambda_{\underline{r}} : (\bar{\epsilon} + \delta\epsilon^{(T_1)}(\underline{r}) - \epsilon^0(\underline{r})) \\ \sigma^{(T_2)}(\underline{r}) = \lambda_{\underline{r}} : (\bar{\epsilon} + \delta\epsilon^{(T_2)}(\underline{r}) - \epsilon^0(\underline{r}))$$

$$4: \quad \hat{\sigma}^{(T_1)}(\underline{q}) = \text{FT}\{\sigma^{(T_1)}(\underline{r})\} \\ \hat{\sigma}^{(T_2)}(\underline{q}) = \text{FT}\{\sigma^{(T_2)}(\underline{r})\}$$

$$5: \quad \hat{u}(\underline{q}) \leftarrow \hat{u}(\underline{q}) - \hat{G}^0(\underline{q}) \cdot \{\hat{\sigma}^{(T_1)}(\underline{q}) \cdot \overline{\underline{D}(\underline{q})} - \hat{\sigma}^{(T_2)}(\underline{q}) \cdot \underline{D}(\underline{q})\} / 2$$

$$6: \quad \text{if applied stress:} \\ \delta\epsilon(\underline{r}) = \frac{1}{2}(\delta\epsilon^{(T_1)}(\underline{r}) + \delta\epsilon^{(T_2)}(\underline{r})) \\ \bar{\epsilon} = (\lambda^0)^{-1} \cdot \{\sigma^a + (\lambda^0 - \langle \lambda(\underline{r}) \rangle) : \bar{\epsilon} = \langle \lambda(\underline{r}) : \delta\epsilon(\underline{r}) \rangle + \langle \lambda^0(\underline{r}) : \epsilon^0(\underline{r}) \rangle\}$$

7: return to 1 until Convergence

We note that, as a consequence of property (2.39), these Green operators display the following symmetry relations:

$$\begin{gathered} \hat{\Gamma}_{ijkl}^{0,d}(\underline{q}) = \hat{\Gamma}_{jikl}^{0,d}(\underline{q}) = \hat{\Gamma}_{ijkl}^{0,d}(\underline{q}) = \overline{\hat{\Gamma}_{klij}^{0,d}(\underline{q})}, \\ \hat{\Gamma}_{ijkl}^{0,nd}(\underline{q}) = \hat{\Gamma}_{jikl}^{0,nd}(\underline{q}) = \hat{\Gamma}_{ijkl}^{0,nd}(\underline{q}) = \overline{\hat{\Gamma}_{klij}^{0,nd}(\underline{q})}. \end{gathered} \tag{2.58}$$

Using the relations,

$$\begin{gathered} \frac{1}{2} \{ \hat{\Gamma}_{\underline{\epsilon}}^{0,d}(\underline{q}) : \lambda_{\underline{\epsilon}}^0 : \hat{\delta}\epsilon^{(T_1)}(\underline{q}) + \hat{\Gamma}_{\underline{\epsilon}}^{0,nd}(\underline{q}) : \lambda_{\underline{\epsilon}}^0 : \hat{\delta}\epsilon^{(T_2)}(\underline{q}) \} = \hat{\delta}\epsilon^{(T_1)}(\underline{q}), \\ \frac{1}{2} \{ \hat{\Gamma}_{\underline{\epsilon}}^{0,nd}(\underline{q}) : \lambda_{\underline{\epsilon}}^0 : \hat{\delta}\epsilon^{(T_1)}(\underline{q}) + \hat{\Gamma}_{\underline{\epsilon}}^{0,d}(\underline{q}) : \lambda_{\underline{\epsilon}}^0 : \hat{\delta}\epsilon^{(T_2)}(\underline{q}) \} = \hat{\delta}\epsilon^{(T_2)}(\underline{q}), \end{gathered} \tag{2.59}$$

which, because $\hat{\delta}\epsilon^{(T_1)}(\underline{q})$ and $\hat{\delta}\epsilon^{(T_2)}(\underline{q})$ are equal to zero when $\underline{q} = 0$, are valid for any q-vector, the iterative scheme (2.56) can be written into a commonly used Lippmann-Schwinger form,

$$\begin{aligned} \hat{\delta}\epsilon^{(T_1)}(\underline{q}) &= -\frac{1}{2} \left\{ \hat{\Gamma}_{\underline{\epsilon}}^{0,d}(\underline{q}) : \hat{\zeta}_{(n)}^{(T_1)}(\underline{q}) + \hat{\Gamma}_{\underline{\epsilon}}^{0,nd}(\underline{q}) : \hat{\zeta}_{(n)}^{(T_2)}(\underline{q}) \right\}, \\ \hat{\delta}\epsilon^{(T_2)}(\underline{q}) &= -\frac{1}{2} \left\{ \overline{\hat{\Gamma}_{\underline{\epsilon}}^{0,nd}(\underline{q})} : \hat{\zeta}_{(n)}^{(T_1)}(\underline{q}) + \overline{\hat{\Gamma}_{\underline{\epsilon}}^{0,d}(\underline{q})} : \hat{\zeta}_{(n)}^{(T_2)}(\underline{q}) \right\}, \end{aligned} \tag{2.60}$$

where the polarization tensors $\hat{\sigma}^{(T_1)}(\underline{q})$ and $\hat{\sigma}^{(T_2)}(\underline{q})$ are defined by

$$\hat{\zeta}^{(T_1)}(\underline{q}) = \hat{\sigma}^{(T_1)}(\underline{q}) - \lambda_{\underline{\epsilon}}^0 : \hat{\epsilon}^{(T_1)}(\underline{q}), \quad i = 1, 2. \tag{2.61}$$

Finally, the iterative scheme (2.56) may also be formulated as

$$\begin{aligned} \hat{\delta}\epsilon_{(n+1)}^{(T_1)}(\underline{q}) &= \hat{\delta}\epsilon_{(n)}^{(T_1)}(\underline{q}) - \hat{\Xi}_{\underline{\epsilon}}^{0}(\underline{q}) \cdot \frac{1}{2} \{ \hat{\sigma}_{(n)}^{(T_1)}(\underline{q}) \cdot \overline{\underline{D}(\underline{q})} - \hat{\sigma}_{(n)}^{(T_2)}(\underline{q}) \cdot \underline{D}(\underline{q}) \}, \\ \hat{\delta}\epsilon_{(n+1)}^{(T_2)}(\underline{q}) &= \hat{\delta}\epsilon_{(n)}^{(T_2)}(\underline{q}) + \hat{\Xi}_{\underline{\epsilon}}^{0}(\underline{q}) \cdot \frac{1}{2} \{ \hat{\sigma}_{(n)}^{(T_1)}(\underline{q}) \cdot \overline{\underline{D}(\underline{q})} - \hat{\sigma}_{(n)}^{(T_2)}(\underline{q}) \cdot \underline{D}(\underline{q}) \}, \end{aligned} \tag{2.62}$$

where the third-order tensor $\hat{\Xi}_{\underline{\epsilon}}^0(\underline{q})$ is defined by

$$\hat{\Xi}_{ijk}^{0}(\underline{q}) = \frac{1}{2} \{ D_{l}(\underline{q}) \hat{G}_{jk}^{0}(\underline{q}) + D_{j}(\underline{q}) \hat{G}_{lk}^{0}(\underline{q}) \}. \tag{2.63}$$

Note that this tensor is symmetric with respect to the first two indices. The advantage of algorithm (2.62) over (2.56) and (2.60) is that only one Green operator is needed and that it explicitly uses the quantity which, measuring in q-space the deviation from mechanical equilibrium (see Eq. (2.25)), is needed to test the convergence of the iterative scheme. Finally, the strain-based algorithm that results from (2.62) is summarized in Algorithm 2. If formulation (2.56) is preferred, equations in step (4) are replaced by (2.56).
<caption>**Algorithm 2**: strain-based basic scheme for the tetrahedral stencil on a simple cubic grid</caption>

initial guess for $\hat{\delta}_{\tilde{\epsilon}}^{(T_1)}(\underline{q})$ and $\hat{\delta}_{\tilde{\epsilon}}^{(T_2)}(\underline{q})$ (e.g. $\hat{\delta}_{\tilde{\epsilon}}^{(T_1)}(\underline{q}) = 0$, $\hat{\delta}_{\tilde{\epsilon}}^{(T_2)}(\underline{q}) = 0$)
initial guess for $\tilde{\epsilon}$ (e.g. $\tilde{\epsilon} = (\tilde{\lambda}^0)^{-1} : \sigma^a + \langle \epsilon^0(\underline{r}) \rangle$ if applied stress, $\tilde{\epsilon} = \tilde{\epsilon}^a$ if applied strain)

$$1: \delta\tilde{\epsilon}^{(T_1)}(\underline{r}) = \text{FT}^{-1} \left\{ \hat{\delta}_{\tilde{\epsilon}}^{(T_1)}(\underline{q}) \right\} \\ \hat{\delta}_{\tilde{\epsilon}}^{(T_2)}(\underline{r}) = \text{FT}^{-1} \left\{ \hat{\delta}_{\tilde{\epsilon}}^{(T_2)}(\underline{q}) \right\}$$

$$2: \sigma^{(T_1)}(\underline{r}) = \tilde{\lambda}(\underline{r}) : (\tilde{\epsilon} + \delta\epsilon^{(T_1)}(\underline{r}) - \epsilon^0(\underline{r})) \\ \sigma^{(T_2)}(\underline{r}) = \tilde{\lambda}(\underline{r}) : (\tilde{\epsilon} + \delta\epsilon^{(T_2)}(\underline{r}) - \epsilon^0(\underline{r}))$$

$$3: \hat{\sigma}^{(T_1)}(\underline{q}) = \text{FT} \{ \sigma^{(T_1)}(\underline{r}) \} \\ \hat{\sigma}^{(T_2)}(\underline{q}) = \text{FT} \{ \sigma^{(T_2)}(\underline{r}) \}$$

$$4: \hat{\delta}_{\tilde{\epsilon}}^{(T_1)}(\underline{q}) \leftarrow \hat{\delta}_{\tilde{\epsilon}}^{(T_1)}(\underline{q}) - \underline{\hat{\delta}}_{\tilde{\epsilon}}^{0}(\underline{q}). \left\{ \hat{\sigma}_{\tilde{\epsilon}}^{(T_1)}(\underline{q}) \underline{\cdot} \underline{\underline{D}}(\underline{q}) - \hat{\sigma}_{\tilde{\epsilon}}^{(T_2)}(\underline{q}) \underline{\cdot} \underline{\underline{D}}(\underline{q}) \right\} /2 \\ \hat{\delta}_{\tilde{\epsilon}}^{(T_2)}(\underline{q}) \leftarrow \hat{\delta}_{\tilde{\epsilon}}^{(T_2)}(\underline{q}) + \underline{\hat{\delta}}_{\tilde{\epsilon}}^{0}(\underline{q}). \left\{ \hat{\sigma}_{\tilde{\epsilon}}^{(T_1)}(\underline{q}) \underline{\cdot} \underline{\underline{D}}(\underline{q}) - \hat{\sigma}_{\tilde{\epsilon}}^{(T_2)}(\underline{q}) \underline{\cdot} \underline{\underline{D}}(\underline{q}) \right\} /2$$

5: if applied stress :

$$\delta\tilde{\epsilon}(\underline{r}) = \frac{1}{2}(\delta\epsilon^{(T_1)}(\underline{r}) + \delta\epsilon^{(T_2)}(\underline{r})) \\ \tilde{\epsilon} \leftarrow (\tilde{\lambda}^0)^{-1} : \left\{ \sigma^a + (\tilde{\lambda}^0 - \langle \tilde{\lambda}(\underline{r}) \rangle) \right\} : \tilde{\epsilon} - \langle \tilde{\lambda}(\underline{r}) : \delta\epsilon(\underline{r}) \rangle + \langle \tilde{\lambda}(\underline{r}) : \epsilon^0(\underline{r}) \rangle$$

6: return to 1 until Convergence

## 2.7. The tetrahedral stencil and a polarization based method

In this section, we propose an adaptation to the tetrahedral stencil of the Monchiet-Bonnet polarization method [13] based on its reformulation by Sab et al. [44]. This adaptation can be useful when the basic scheme is not efficient enough, which can be the case for example when considering a microstructure containing simultaneously soft and hard inclusions with high contrasts between these inclusions and the matrix. This will be the case in Section 3.5 in which we consider large double contrasts.

Adapted to the tetrahedral stencil, the polarization-based method leads to the following iterative scheme:

$$\delta_{\tilde{\epsilon}}^{(T_1)}(\underline{r}) = \delta_{\tilde{\epsilon}}^{(T_1)}(\underline{r}) + (1 - \gamma) \alpha_{\tilde{\epsilon}}(\underline{r}) : \underline{X}_{\tilde{\epsilon}}^{(T_1)}(\underline{r}), \qquad (2.64)$$

$$\delta_{\tilde{\epsilon}}^{(T_2)}(\underline{r}) = \delta_{\tilde{\epsilon}}^{(T_2)}(\underline{r}) + (1 - \gamma) \alpha_{\tilde{\epsilon}}(\underline{r}) : \underline{X}_{\tilde{\epsilon}}^{(T_2)}(\underline{r}).$$

The so-called residuals $\underline{X}_{\tilde{\epsilon}}^{(T_1)}(\underline{r})$ and $\underline{X}_{\tilde{\epsilon}}^{(T_2)}(\underline{r})$ are given by

$$\underline{X}_{\tilde{\epsilon}}^{(T_1)}(\underline{r}) = -\delta_{\tilde{\epsilon}}^{(T_1)}(\underline{r}) - \text{FT}^{-1} \left\{ \hat{\underline{X}}_{\tilde{\epsilon}}^{0,d}(\underline{q}) : \hat{\underline{\epsilon}}^{(T_1)}(\underline{q}) + \hat{\underline{X}}_{\tilde{\epsilon}}^{0,nd}(\underline{q}) : \hat{\underline{\epsilon}}^{(T_2)}(\underline{q}) \right\} /2, \qquad (2.65)$$

$$\underline{X}_{\tilde{\epsilon}}^{(T_2)}(\underline{r}) = -\delta_{\tilde{\epsilon}}^{(T_2)}(\underline{r}) - \text{FT}^{-1} \left\{ \hat{\underline{X}}_{\tilde{\epsilon}}^{0,nd}(\underline{q}) : \hat{\underline{\epsilon}}^{(T_1)}(\underline{q}) + \hat{\underline{X}}_{\tilde{\epsilon}}^{0,d}(\underline{q}) : \hat{\underline{\epsilon}}^{(T_2)}(\underline{q}) \right\} /2,$$

where the diagonal and non-diagonal Green operators $\hat{\underline{X}}_{\tilde{\epsilon}}^{0,d}(\underline{q})$ and $\hat{\underline{X}}_{\tilde{\epsilon}}^{0,nd}(\underline{q})$ have been introduced in (2.57), the polarization tensors $\hat{\underline{\epsilon}}^{(T_1)}(\underline{q})$ and $\hat{\underline{\epsilon}}^{(T_2)}(\underline{q})$ have been defined in (2.61) and the tensor $\underline{\alpha}_{\tilde{\epsilon}}(\underline{r})$ is linked to the elastic tensor $\underline{\lambda}_{\tilde{\epsilon}}(\underline{r})$ by

$$\underline{\alpha}_{\tilde{\epsilon}}(\underline{r}) = 2 (\tilde{\lambda}_{\tilde{\epsilon}}(\underline{r}) + \tilde{\lambda}_{0}(\underline{r}))^{-1} : \underline{\lambda}_{0}(\underline{r}). \qquad (2.66)$$

The corresponding algorithm is shown in Algorithm 3. As pointed out in [12,13,45], the polarization scheme contains the Eyre-Milton and augmented Lagrangian methods, which are recovered for $\gamma = 0$ and $\gamma = 1/2$, respectively, while, based on numerical experiences, Monchiet-Bonnet [13] recommend the value $\gamma = 1/4$.

## 2.8. Convergence criteria

Convergence to equilibrium can be evaluated using the $L^2$-norm of the quantity that measures deviation from equilibrium, i.e.

$$L_2 = \| N \frac{\partial E}{\partial \underline{u}} \|_{L^2} = \sqrt{ \frac{1}{N} \sum_{\underline{r} \in Y_1 \cup Y_2} \left| N \frac{\partial E}{\partial \underline{u}(\underline{r})} \right|^2 }, \qquad (2.67)$$

where the quantity $N \partial E / \partial \underline{u}(\underline{r})$ represents, in units of $1/d$, the local force density that acts on site $\underline{r}$ and $|\cdot|$ denotes the euclidean norm of a vector. Using Parseval's theorem, we get

$$L_2 = \sqrt{\sum_{\underline{q}} \left| \frac{\partial E}{\partial \hat{\underline{u}}(\underline{q})} \right|^2}. \qquad (2.68)$$
**Algorithm 3**: strain-based polarization method for the tetrahedral stencil

initial guess for $\delta\epsilon^{(T_1)}(\underline{\boldsymbol{r}})$ and $\delta\epsilon^{(T_2)}(\underline{\boldsymbol{r}})$ (e.g. $\delta\epsilon^{(T_1)}(\underline{\boldsymbol{r}}) = 0$, $\delta\epsilon^{(T_2)}(\underline{\boldsymbol{r}}) = 0$)

initial guess for $\tilde{\epsilon}$ (e.g. $\tilde{\epsilon} = (\tilde{\lambda}^0)^{-1} : \sigma^a + \langle \epsilon^0(\underline{\boldsymbol{r}}) \rangle$ if applied stress, $\tilde{\epsilon} = \tilde{\epsilon}^a$ if applied strain)

$$ 
\begin{aligned}
1: \quad & \sigma^{(T_1)}(\underline{\boldsymbol{r}}) = \tilde{\lambda}^{(T_1)}(\underline{\boldsymbol{r}}) : (\tilde{\epsilon} + \delta\epsilon^{(T_1)}(\underline{\boldsymbol{r}}) - \epsilon^0(\underline{\boldsymbol{r}})) \\
& \tilde{\sigma}^{(T_2)}(\underline{\boldsymbol{r}}) = \tilde{\lambda}^{(T_2)}(\underline{\boldsymbol{r}}) : (\tilde{\epsilon} + \delta\epsilon^{(T_2)}(\underline{\boldsymbol{r}}) - \epsilon^0(\underline{\boldsymbol{r}}))
\end{aligned}
$$

$$ 
\begin{aligned}
2: \quad & \tilde{\boldsymbol{\tau}}^{(T_1)}(\underline{\boldsymbol{r}}) = \sigma^{(T_1)}(\underline{\boldsymbol{r}}) - \tilde{\lambda}^{(T_1)}: (\tilde{\epsilon} + \delta\epsilon^{(T_1)}(\underline{\boldsymbol{r}})) \\
& \tilde{\boldsymbol{\tau}}^{(T_2)}(\underline{\boldsymbol{r}}) = \sigma^{(T_2)}(\underline{\boldsymbol{r}}) - \tilde{\lambda}^{(T_2)}: (\tilde{\epsilon} + \delta\epsilon^{(T_2)}(\underline{\boldsymbol{r}}))
\end{aligned}
$$

$$ 
\begin{aligned}
3: \quad & \hat{\tilde{\boldsymbol{\epsilon}}}^{(T_1)}(\underline{\boldsymbol{q}}) = \text{FT} \{ \tilde{\boldsymbol{\epsilon}}^{(T_1)}(\underline{\boldsymbol{r}}) \} \\
& \hat{\tilde{\boldsymbol{\epsilon}}}^{(T_2)}(\underline{\boldsymbol{q}}) = \text{FT} \{ \tilde{\boldsymbol{\epsilon}}^{(T_2)}(\underline{\boldsymbol{r}}) \}
\end{aligned}
$$

$$ 
\begin{aligned}
4: \quad & \tilde{\boldsymbol{X}}^{(T_1)}(\underline{\boldsymbol{r}}) = -\delta\tilde{\epsilon}^{(T_1)}(\underline{\boldsymbol{r}}) - \text{FT}^{-1} \left\{ \hat{\tilde{\boldsymbol{\Gamma}}}^{0,d}_{\tilde{\sigma}}(\underline{\boldsymbol{q}}) : \hat{\tilde{\boldsymbol{\epsilon}}}^{(T_1)}(\underline{\boldsymbol{q}}) + \hat{\tilde{\boldsymbol{\Gamma}}}^{0,nd}_{\tilde{\sigma}}(\underline{\boldsymbol{q}}) : \hat{\tilde{\boldsymbol{\epsilon}}}^{(T_2)}(\underline{\boldsymbol{q}}) \right\} / 2 \\
& \tilde{\boldsymbol{X}}^{(T_2)}(\underline{\boldsymbol{r}}) = -\delta\tilde{\epsilon}^{(T_2)}(\underline{\boldsymbol{r}}) - \text{FT}^{-1} \left\{ \hat{\tilde{\boldsymbol{\Gamma}}}^{0,nd}_{\tilde{\sigma}}(\underline{\boldsymbol{q}}) : \hat{\tilde{\boldsymbol{\epsilon}}}^{(T_1)}(\underline{\boldsymbol{q}}) + \hat{\tilde{\boldsymbol{\Gamma}}}^{0,d}_{\tilde{\sigma}}(\underline{\boldsymbol{q}}) : \hat{\tilde{\boldsymbol{\epsilon}}}^{(T_2)}(\underline{\boldsymbol{q}}) \right\} / 2
\end{aligned}
$$

$$ 
\begin{aligned}
5: \quad & \delta\tilde{\epsilon}^{(T_1)}(\underline{\boldsymbol{r}}) \leftarrow \delta\tilde{\epsilon}^{(T_1)}(\underline{\boldsymbol{r}}) + (1-\gamma) \underline{\boldsymbol{\alpha}}^{(T_1)}(\underline{\boldsymbol{r}}) : \tilde{\boldsymbol{X}}^{(T_1)}(\underline{\boldsymbol{r}}) \\
& \delta\tilde{\epsilon}^{(T_2)}(\underline{\boldsymbol{r}}) \leftarrow \delta\tilde{\epsilon}^{(T_2)}(\underline{\boldsymbol{r}}) + (1-\gamma) \underline{\boldsymbol{\alpha}}^{(T_2)}(\underline{\boldsymbol{r}}) : \tilde{\boldsymbol{X}}^{(T_2)}(\underline{\boldsymbol{r}})
\end{aligned}
$$

6: if applied stress :

$$ 
\begin{aligned}
\delta\tilde{\epsilon}(\underline{\boldsymbol{r}}) &= \frac{1}{2} (\delta\epsilon^{(T_1)}(\underline{\boldsymbol{r}}) + \delta\epsilon^{(T_2)}(\underline{\boldsymbol{r}})) \\
\tilde{\epsilon} &\leftarrow (\tilde{\lambda}^0)^{-1} : \left\{ \tilde{\sigma}^a + \left( \tilde{\lambda}^0 - \langle \tilde{\lambda}(\underline{\boldsymbol{r}}) \rangle \right) : \tilde{\epsilon} - \langle \tilde{\lambda}(\underline{\boldsymbol{r}}) : \delta\epsilon(\underline{\boldsymbol{r}}) \rangle + \langle \tilde{\lambda}(\underline{\boldsymbol{r}}) : \epsilon^0(\underline{\boldsymbol{r}}) \rangle \right\}
\end{aligned}
$$

7: return to 1 until Convergence

In Appendix B, we show that, provided the divergence operator is defined in a proper way, the derivative $\partial E / \partial \hat{\underline{\boldsymbol{u}}}(\underline{\boldsymbol{q}})$ is related to the divergence of the average stress $\sigma$ defined in (2.49a) by

$$ \widehat{\text{div}} \tilde{\sigma}(\underline{\boldsymbol{q}}) = - \frac{\partial E}{\partial \hat{\underline{\boldsymbol{u}}}(\underline{\boldsymbol{q}})}. \qquad (2.69) $$

Therefore, as one could have anticipated, the present criterion amounts simply to the $L^2$-norm of the stress divergence,

$$ L_2 = \sqrt{\sum_{\underline{\boldsymbol{q}}} |\widehat{\text{div}} \tilde{\sigma}(\underline{\boldsymbol{q}})|^2} = \|\text{div} \sigma\|_{L^2}. \qquad (2.70) $$

Finally, to make it dimensionless and insensitive to the amplitude of the stress, we define the convergence criterion associated to the $L^2$-norm of the stress divergence by

$$ \epsilon_{L^2} = \frac{\|\text{div} \sigma\|_{L^2}}{\|\langle \sigma \rangle\|}, \qquad (2.71) $$

where $\|\langle \sigma \rangle\| = \sqrt{\sum_{i,j} |\langle \sigma_{ij} \rangle|^2}$ is the Frobenius norm of the average stress. According to (2.69) and (2.25), $\epsilon_{L^2}$ can be numerically computed at iteration step n as

$$ \epsilon_{L^2}^{(n)} = \frac{\cdot \sqrt{\frac{1}{4} \sum_q |\hat{\tilde{\sigma}}_{(n)}^{(T_1)}(\underline{\boldsymbol{q}}) \cdot \overline{\underline{\boldsymbol{D}}}(\underline{\boldsymbol{q}}) - \hat{\tilde{\sigma}}_{(n)}^{(T_2)}(\underline{\boldsymbol{q}}) \cdot \underline{\boldsymbol{D}}}(\underline{\boldsymbol{q}})|^2}{\| \langle \sigma_{(n)} \rangle \|} \qquad (2.72) $$

Other convergence criteria have been discussed in the literature, often in link with the method used to solve the equilibrium equations (see for example [28,44–46]). As shown in Section 2.3, the basic scheme implementation of the tetrahedral stencil can be formulated as a gradient descent method within the appropriate Hilbertian structure. In this context, a sensible convergence criterion may be given by the norm of the energy gradient $\|\underline{\boldsymbol{\nabla}} E[\underline{\boldsymbol{u}}]\|_{H_{\text{per}}^1}$, where $\|\cdot\|_{H_{\text{per}}^1}$ is the norm induced by the scalar product defined in (2.31), i.e.

$$ \|\underline{\boldsymbol{\nabla}} E[\underline{\boldsymbol{u}}]\|_{H_{\text{per}}^1} = \sqrt{\langle \underline{\boldsymbol{\nabla}} E[\underline{\boldsymbol{u}}], \underline{\boldsymbol{\nabla}} E[\underline{\boldsymbol{u}}] \rangle_{H_{\text{per}}^1}}. \qquad (2.73) $$

According to the Riesz mapping from the dual $H_{\text{per}}^{-1}$ to $H_{\text{per}}^1$, which was identified above as being equal to $\tilde{\sigma}^0$, we have

$$ \|\underline{\boldsymbol{\nabla}} E[\underline{\boldsymbol{u}}]\|_{H_{\text{per}}^1} = \|\tilde{\boldsymbol{G}}^{0-1} \underline{\boldsymbol{\nabla}} E[\underline{\boldsymbol{u}}]\|_{H_{\text{per}}^{-1}}, \qquad (2.74) $$
where $||.||_{H_{per}^{-1}}$ is the norm in $H_{per}^{-1}$ induced by the norm in $H_{per}^{1}$ and $G^{0-1}$ acts as a convolution operator. Using (2.42) and (2.69), one gets

$$ \underline{\nabla} E[\underline{u}](\underline{r}) = -\underline{G}^0 \operatorname{div} \underline{\sigma}(\underline{r}), \qquad (2.75) $$

and therefore

$$ \| \underline{\nabla} E[\underline{u}] \|_{H_{per}^{1}} = \| \operatorname{div} \underline{\sigma} \|_{H_{per}^{-1}}, \qquad (2.76) $$

which shows that the $H_{per}^{1}$-norm of the energy gradient $\underline{\nabla} E[\underline{u}]$ is equal to the $H_{per}^{-1}$-norm of the divergence of the stress when $\operatorname{div} \underline{\sigma}$ is regarded as an operator that lives in $H_{per}^{-1}$. Finally, to make it dimensionless and insensitive to the amplitude of the average stress, we define the convergence criterion associated to the $H_{per}^{1}$-norm of the energy gradient, or equivalently to the $H_{per}^{-1}$-norm of the stress divergence, by

$$ \epsilon_{H_{per}^{-1}} = \frac{\|\operatorname{div} \underline{\sigma}\|_{H_{per}^{-1}}}{\sqrt{\langle \underline{\sigma} \rangle : \underline{\lambda}^{0-1} : \langle \underline{\sigma} \rangle}}. \qquad (2.77) $$

Inserting the Fourier representation (2.43) of the energy gradient into the q-space formulation (2.33) of the scalar product, one gets that $\epsilon_{H_{per}^{-1}}$ can be numerically computed at iteration step n as

$$ \epsilon_{H_{per}^{-1}}^{(n)} = \frac{\sqrt{\frac{1}{4} \sum_q \{\hat{\underline{\sigma}}_{(n)}^{(T_1)}(\underline{q}) \cdot \overline{\underline{D}(\underline{q})} - \hat{\underline{\sigma}}_{(n)}^{(T_2)}(\underline{q}) \cdot \underline{D}(\underline{q})\} \cdot \hat{\underline{G}}^0(\underline{q}) \cdot \{\hat{\underline{\sigma}}_{(n)}^{(T_1)}(\underline{q}) \cdot \underline{D}(\underline{q}) - \hat{\underline{\sigma}}_{(n)}^{(T_2)}(\underline{q}) \cdot \underline{D}(\underline{q})\}}}{\sqrt{\langle \underline{\sigma}_{(n)} \rangle : \underline{\lambda}^{0-1} : \langle \underline{\sigma}_{(n)} \rangle}}, \qquad (2.78) $$

or, because within the basic scheme, the increment on the displacement field is equal to the negative of the energy gradient, as

$$ \epsilon_{H_{per}^{-1}}^{(n)} = \frac{\sqrt{\sum_q \{\hat{\underline{u}}_{(n+1)}(\underline{q}) - \hat{\underline{u}}_{(n)}(\underline{q})\} \cdot \overline{\hat{\underline{G}}^0(\underline{q})^{-1}} \cdot \{\hat{\underline{u}}_{(n+1)}(\underline{q}) - \hat{\underline{u}}_{(n)}(\underline{q})\}}}{\sqrt{\langle \underline{\sigma}_{(n)} \rangle : \underline{\lambda}^{0-1} : \langle \underline{\sigma}_{(n)} \rangle}}. \qquad (2.79) $$

If the strain-based basic scheme proposed in Section 2.6 is preferred, the previous expression becomes, using the initial definition (2.31) of the scalar product,

$$ \epsilon_{H_{per}^{-1}}^{(n)} = \frac{\sqrt{\frac{1}{2} \sum_q \left( \overline{\Delta \delta \epsilon_{(n)}^{(T_1)}(\underline{q})} : \underline{\lambda}^0 : \overline{\Delta \delta \epsilon_{(n)}^{(T_1)}(\underline{q})} + \overline{\Delta \delta \epsilon_{(n)}^{(T_2)}(\underline{q})} : \underline{\lambda}^0 : \overline{\Delta \delta \epsilon_{(n)}^{(T_2)}(\underline{q})} \right)}}{\sqrt{\langle \underline{\sigma}_{(n)} \rangle : \underline{\lambda}^{0-1} : \langle \underline{\sigma}_{(n)} \rangle}}, \qquad (2.80) $$

where $\overline{\Delta \delta \epsilon_{(n)}^{(T_1)}(\underline{q})}$ and $\overline{\Delta \delta \epsilon_{(n)}^{(T_2)}(\underline{q})}$ are the increments at step n on the strain fields $\overline{\delta \epsilon_{(n)}^{(T_1)}(\underline{q})}$ and $\overline{\delta \epsilon_{(n)}^{(T_2)}(\underline{q})}$.

We will present in Section 3.6 a comparison of the evolution of the convergence criteria $\epsilon_{L^2}$ and $\epsilon_{H_{per}^{-1}}$ as a function of the iterate index. We will observe that, for the situation considered in this section at least, they behave in a very similar way.

Finally, for the polarization method presented in Section 2.7, we propose to use

$$ \epsilon_{pol} = \frac{\sqrt{\langle \underline{\mathbf{X}}^{(T_1)} : \underline{\lambda}^0 : \underline{\mathbf{X}}^{(T_1)} \rangle + \sqrt{\langle \underline{\mathbf{X}}^{(T_2)} : \underline{\lambda}^0 : \underline{\mathbf{X}}^{(T_2)} \rangle}}}{\sqrt{\overline{\epsilon} : \underline{\lambda}^0 : \overline{\epsilon}}}, \qquad (2.81) $$

which is the adaptation of the $\underline{\lambda}^0$-norm criterion proposed by Sab et al. [44] to the present situation, in which we iterate over the two strain fields $\overline{\delta \epsilon_{(n)}^{(T_1)}}$ and $\overline{\delta \epsilon_{(n)}^{(T_2)}}$ which both converge if their respective residuals $\underline{\mathbf{X}}^{(T_1)}$ and $\underline{\mathbf{X}}^{(T_2)}$ vanish.

### 3. Numerical experiments

In this section, we analyse the performance of the tetrahedral stencil approach for different applications, including comparisons with the original Moulinec–Suquet discretization [1,2] the rotated finite difference scheme proposed by Willot [34] and the staggered grid discretization. Our aim is to test the accuracy and efficiency of the method and also its ability to treat, in term of convergence rate, extreme situations where the elastic contrast is large or infinite.

All the results presented below, with the exception of those in Section 3.4 that concern accelerated schemes and that have been taken from the literature and those in Section 3.5 where we used the polarization method presented in Section 2.7, were obtained with the basic scheme which, for our tetrahedral stencil, was implemented in the form described in Algorithm 1 and the FFT implementation relies upon the FFTW library [47].
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** A cubic inclusion with eigenstrain $\epsilon_{33}^0$ in an isotropic medium with a zero applied stress: internal stresses $\sigma_{11}$, $\sigma_{22}$ and $\sigma_{33}$, in units of $\mu\epsilon_{33}^0$, along the line parallel to direction (100) and that goes through the centre of the cell. The numerical results obtained with the tetrahedral stencil (red symbols) are compared to the analytic solution (blue lines); the right panel shows a zoom on the small zone marked in the left panel. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

### 3.1. A cubic inclusion with eigenstrain in an elastically isotropic homogeneous material

For our first application, following [22], we consider a homogeneous, elastically isotropic medium containing a cubic inclusion characterized by a finite eigenstrain $\epsilon^0$. The interest of this situation is that mechanical equilibrium is known analytically (see for example [48]), which permits to test the accuracy of our numerical method.

In the present situation, we use the same material parameters as in [22]. The media is elastically homogeneous and its shear modulus and Poisson's ratio are given by $\mu = 132.30$ GPa and $\nu = 0.26$. The eigenstrain in the inclusion is such that $\epsilon_{ij}^0 = 0$ if $(ij) \neq (33)$ whereas $\epsilon_{33}^0$ is finite. The periodic medium is discretized with 256³ cells and the cubic inclusion consists of 31³ cells. The boundary condition used here is a zero external stress. The results obtained with the tetrahedral stencil are compared with the exact solution (obtained from [48]) in Fig. 3, where we report the stress components $\sigma_{11}$, $\sigma_{22}$, $\sigma_{33}$ along a line parallel to the cubic direction (100) and that goes through the centre of the cubic inclusion. The comparison with the analytical solution shows that the tetrahedral stencil leads to very accurate numerical results, as observed in particular in the right panel of Fig. 3, where we display a zoom on a small region inside the inclusion.

### 3.2. A cubic void within an elastically isotropic homogeneous matrix

Next, we consider a cubic void embedded into an elastically isotropic matrix. Our aim is to test the behaviour of the tetrahedral stencil in the presence of infinite elastic contrast, a situation which is known to be challenging, both in term of convergence of the iterative algorithms and, if convergence is achieved, in term of accuracy. Indeed, it has long been observed that, in the original Moulinec-Suquet method, the high frequency cut-off generated by the discretization of the Lippmann-Schwinger equation generates oscillations whose amplitude increases with the elastic contrast. It has also been observed that the iterative basic scheme, associated with the Moulinec-Suquet discretization, converges slowly for large contrast and does not converge if the contrast is infinite.

We consider a material with the same elastic parameters as in the previous case. The unit cell is discretized into 64³ voxels and we insert a cubic cavity which consists of 31³ cells and whose faces are perpendicular to the cartesian axis (this relatively small size has been selected since it may typically correspond to linear dimensions of heterogeneities in a complex microstructure). The overall system is subjected to a hydrostatic pressure $P = 300$ MPa.

In the following, we will compare the results obtained with the tetrahedral stencil to those obtained with the Moulinec-Suquet discretization and the rotated scheme used by Willot. We mention that, for these two methods and when needed, we use the regularization recipes proposed by the corresponding authors to fix problems encountered with the behaviour of the respective Green operators. Indeed, within the Moulinec-Suquet discretization, if the number of voxels along cartesian axis $i$ is even, conjugation relations on the Fourier representation of the Green function are lost when the component $q_i$ of a q-vector reaches the highest frequency, i.e. when $q_i = \pi$, which makes the imaginary part of backward Fourier transforms of the strain field non zero. Following the recipe proposed by Moulinec-Suquet [2], we then fix the strain Green function equal to the inverse of the elastic tensor $\lambda^0$. As for the rotated scheme, if the number of voxels along any two of the cartesian axis, say $i$ and $j$, are even, all the components of the q-space differential operator $\mathcal{D}(q)$ are equal to zero when components $q_i$ and $q_j$ of the q-vector reach the highest frequency, i.e. when $q_i = q_j = \pi$, which makes the Green function undefined at the corresponding q-lines. Following the recipe proposed by Willot, we then fix the displacement (or strain) Green function to zero. We recall that this problem does not exist with our discretization scheme, because the tetrahedral stencil has the property that the q-space differential operator $\mathcal{D}(q)$ never vanishes for $q \neq 0$, except at a q-point equivalent to the origin of Fourier space (see discussion at the end of Section 2.3).

First, we report in Fig. 4 the evolution of the error on equilibrium $c_{L2}$, defined in Eq. (2.71), as a function of the iterate index for the tetrahedral stencil and, for comparison, for the rotated and Moulinec-Suquet discretizations. The stiffness tensor $\lambda_{\approx}^0$ of the reference medium has not been optimized and, for the three methods, has been simply set to $\lambda_{\approx}^0 = 0.8 \lambda_{\approx}^m$ where $\lambda_{\approx}^m$ is the stiffness of
<figure><img src="image_5.png" /><figcaption>Fig. 4. A cubic void within an isotropic medium subjected to a hydrostatic pressure $P = 300$ MPa: error to equilibrium $ϵ_{L2}$ as a function of iterates within the basic scheme. The dotted line indicates the precision required to reach, with the tetrahedral stencil, an accuracy of 1 MPa on the stress field.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. A cubic void within an isotropic matrix: error $ϵ_{L2}$ vs. iteration within the basic scheme for different values of the reference stiffness tensor $\lambda_{\alpha}^0 = \alpha \lambda_{\alpha}^m$.</figcaption></figure>

the matrix. We observe that the tetrahedral and rotated schemes converge very rapidly, with the highest precision investigated here, $ϵ_{L2} = 10^{-10}$, being reached with less than 100 iterates, with a slight advantage for the tetrahedral stencil, whereas the error with the Moulinec–Suquet discretization decreases very slowly and, in fact, does not appear to converge to zero, a fact that has already been observed when the elastic contrast is infinite, as in the present case. The highest precision investigated here is much more than necessary in common applications. In Fig. 4, the dotted line indicates the precision required to reach, with the tetrahedral stencil, an accuracy of 1 MPa on the stress field, accuracy been defined as the maximum difference observed between the calculated stress field and that obtained in a high-precision reference calculation ($ϵ_{L2} = 10^{-10}$). We observe that only approximately 20 iterations are needed to achieve the required accuracy.

We also analysed the influence of the reference stiffness tensor $\lambda_{\alpha}^0$ on the convergence to equilibrium. Setting $\lambda_{\alpha}^0 = \alpha \lambda_{\alpha}^m$, we report in Fig. 5 the evolution of the error as a function of the iterate for different values of $\alpha$. We observe that the precision $ϵ_{L2} = 10^{-10}$ is always reached for $\alpha \ge 0.5$ and that convergence is then linear. For $\alpha < 0.5$, the lowest achievable accuracy increases rapidly with decreasing $\alpha$. Overall, the optimal convergence rate is obtained for $\alpha = 0.5$, i.e. when the reference stiffness tensor is equal to the arithmetic average of the matrix and void stiffness tensors, as expected for isotropic materials [11,12]. We also note that, for $\alpha \ge 0.5$, the number of iterates to reach equilibrium varies linearly with $\alpha$ and we mention that, within a precision smaller than $10^{-6}$ MPa, the stress field at equilibrium is insensitive to the reference stiffness tensors $\lambda_{\alpha}^0$.

Next, we display in Fig. 6 the stress components $σ_{11}$, $σ_{22}$ and the von Mises stress along a line parallel to direction (100) and passing through the centre of the cubic void. Results obtained with the tetrahedral stencil are compared to those obtained with the Moulinec–Suquet and rotated discretizations. For the last two methods, computations have also been performed with a computational box whose linear dimensions along the cartesian axis are odd, namely a box that consists of 65³ voxels, for which no regularization recipes are needed. For all the computations, the stiffness tensor $\lambda_{\alpha}^0$ of the reference homogeneous medium has been chosen equal to 80% of that of the matrix. The results reported for the tetrahedral and rotated stencils are those obtained when the precision $ϵ_{L2}$ has reached $10^{-10}$, whereas, because of the slow convergence, those obtained with the Moulinec–Suquet discretization corresponds
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** A cubic void within an isotropic medium subjected to a hydrostatic pressure $P = 300$ MPa : internal stresses $\sigma_{11}$, $\sigma_{22}$ and von Mises stress (MPa) along the line parallel to direction (100) and passing through the centre of the void. The numerical results obtained with the tetrahedral stencil (red lines) are compared to the results obtained with the rotated (green lines) and Moulinec–Suquet (blue lines) discretization; the cubic cavity consists of $31^3$ voxels; full lines represent results obtained with a unit cell containing $64^3$ voxels and dotted lines those obtained with a unit cell containing $65^3$ voxels; the middle and bottom panels show zooms on the zones marked 1 and 2, respectively, in the top panel. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

to a precision equal to $10^{-4}$. As seen in Fig. 6, the results obtained with the tetrahedral stencil are exempt of any artefact, whereas significant oscillations are observed with the Moulinec–Suquet discretization and also with the rotated scheme, whether the linear dimensions of the computational box are even or odd.

Clearly, the recipe proposed by Willot to correct for the lack of definition of the Green function when the box dimensions are even is unsatisfactory, as it leaves the Green function singular at q-points where it should be analytic, which in turn generates oscillations in real space. As seen in Fig. 6, these oscillations are hardly suppressed when the box dimensions are odd, even though the Green function in that case is defined for any finite q-vector. The reason why these oscillations are always present, whether the linear dimensions are even or odd, is explained in Appendix D.

Finally, we show in Fig. 7 the map of the von Mises stress obtained with the tetrahedral stencil and, for comparison, those obtained with Moulinec–Suquet and Willot discretizations. We clearly observed oscillations with the Moulinec–Suquet discretization, as well as checkerboard effects with the rotated stencil, whereas the stress map obtained with the tetrahedral stencil is again exempt of any artefact.
<figure><img src="image_8.png" /><figcaption>Fig. 7. A cubic void within an isotropic medium subjected to a hydrostatic pressure $P = 300$ MPa : von Mises stress maps along the (001) plane passing through the centre of computational box; comparison between the results obtained with the Moulinec–Suquet discretization (top left panel), the rotated discretization (top right panel) and the tetrahedral stencil (bottom panel) (units: MPa). The computational box is discretized into $64^3$ cells and the cubic void consists of $31^3$ voxels.</figcaption></figure>

### 3.3. A hard spherical inclusion within a soft matrix

We consider here a situation opposite to the previous one, i.e. a soft matrix reinforced by hard precipitates. The microstructure consists of a single glass inclusion, with Young's modulus and Poisson coefficient given by $E = 72$ GPa and $v = 0.22$, embedded into a polyamide matrix characterized by $E = 2.1$ GPa and $v = 0.3$. The periodic medium consists of $128^3$ cells and the radius of the spherical inclusion is, in cell units, $R = 40$. The overall system is subjected to the uniaxial external strain $\bar{e}_{11} = 0.05$.

We show in Fig. 8 results concerning the von Mises stress obtained with the tetrahedral stencil and, for comparison, the stresses obtained with Moulinec–Suquet and rotated discretizations. The stiffness tensor $\lambda_{\bar{m}}^0$ of the reference medium has been set to $\lambda_{\bar{m}}^0 = 0.5 (\lambda_{\bar{m}}^m + \lambda_{\bar{m}}^p)$ where $\lambda_{\bar{m}}^m$ and $\lambda_{\bar{m}}^p$ are the stiffness matrices of matrix and precipitate, respectively. The results reported for the tetrahedral and rotated stencils are those obtained when the precision $\epsilon_L^2$ defined in (2.71) has reached $10^{-10}$, whereas, because of the slow convergence, those obtained with the Moulinec–Suquet discretization corresponds to a precision equal to $10^{-8}$. Strong oscillations within the inclusion are observed with the Moulinec–Suquet method, which are hardly reduced with the rotated scheme, whereas the results obtained with the tetrahedral stencil are again exempt of any artefacts. The same conclusion are obtained when we compare 2D maps, as shown in Fig. 9, where we observe ringing artefacts with the Moulinec–Suquet discretization and checkerboarding with the rotated scheme, whereas the tetrahedral stencil map is free of these artefacts.

Finally, for the interested reader, we compare in Fig. 10 the fluctuating displacement fields $u(r)$ obtained with the tetrahedral and rotated stencils. Maps are displayed along a (100) plane perpendicular to the applied strain $\bar{e}_{11}$ and close to the centre of the inclusion. For the component $u_1(r)$, the rotated stencil result displays strong checkerboarding within the matrix near but also relatively
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** A hard spherical inclusion within a soft isotropic matrix subjected to the uniaxial external strain $\bar{\epsilon}_{11} = 0.05$ : internal von Mises stress (MPa) along the line parallel to direction (100) and passing through the centre of the inclusion. Numerical results obtained with the tetrahedral stencil (red lines) are compared to those obtained with the rotated (green lines) and Moulinec-Suquet (blue lines) discretizations. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** A hard spherical inclusion within a soft isotropic matrix subjected to the uniaxial external strain $\bar{\epsilon}_{11} = 0.05$ : von Mises stress maps along the (001) plane that goes through the centre of computational box; comparison between the results obtained with the Moulinec-Suquet discretization (top left panel), the rotated scheme (top, right panel) and the tetrahedral stencil (bottom panel) (units: MPa). The computational box consists of 64³ voxels.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** A hard spherical inclusion within a soft isotropic matrix subjected to the uniaxial external strain $\bar{\epsilon}_{11} = 0.05$: fluctuating displacement field maps along the (100) plane whose centre is offset by half a unit cell below the centre of the spherical hard inclusion; comparison between the results obtained with the tetrahedral stencil (left panels) and the rotated discretization (right panels) for the component $u_1(\mathbf{r})$ (top row) and $u_2(\mathbf{r})$ (bottom row) (units: cubic voxel length $d$). Remark: the displacement field $\underline{u}(\mathbf{r})$ is defined on sites $\mathbf{r}$ whose coordinates are half-integers (see Fig. 1); therefore, since Fourier transforms in FFT libraries involve integers, the computation of the displacement field in real space by inverse Fourier transform is carried out as follows: $\underline{u}(lmn - \frac{111}{2}) = \sum_q \exp(-i(q_1 + q_2 + q_3)/2) \hat{u}(\underline{q}) \exp(i\underline{q} \cdot (lmn))$.</figcaption></figure>

far from the interfaces with the hard inclusion, whereas the solution field predicted by the tetrahedral stencil is smooth. For the components $u_2(\mathbf{r})$ and $u_3(\mathbf{r})$ (which is not shown because, to a rotation of $\pi/2$ around the X-axis, it is identical to $u_2(\mathbf{r})$), the maps generated by the tetrahedral and rotated stencil are similar.

### 3.4. Rate of convergence vs. contrast

We present here a brief analysis of the convergence rate of the tetrahedral solver as a function of the elastic contrast. Our aim is to test the performance of the tetrahedral stencil, associated to the basic iterative scheme, in presence of large contrasts.

The microstructure consists of a single spherical inclusion whose radius, in cell units, is $R = 17.6$, embedded into a periodic medium which consists of $44^3$ cells. Matrix and inclusion are elastically isotropic with the same Poisson coefficient $\nu = 0.125$ and a large range of elastic contrasts $C = k_{inclusion}/k_{matrix}$ between the bulk moduli of inclusion and matrix is considered. The overall system is subjected to the uniaxial external strain $\bar{\epsilon}_{11} = 1$.

We report in Fig. 11 the number of iterates needed to reach equilibrium as a function of the contrast $C$. Precision on equilibrium has been fixed to $e_{L2} = 10^{-10}$. Results obtained with the tetrahedral stencil are compared to those obtained with the rotated discretization, the original Moulinec–Suquet method and three accelerated schemes that have been proposed to overcome the slow convergence observed with the Moulinec–Suquet basic scheme, namely the augmented Lagrangian method, the Eyre–Milton approach and the polarization-based method proposed by Monchiet and Bonnet. The results concerning the last three methods have been taken from Ref. [12]. For the tetrahedral stencil, the rotated discretization and Moulinec–Suquet discretization, we used an isotropic reference medium that has the same Poisson coefficient as the inclusion and matrix and a bulk modulus given by $k_0 = \alpha k_{matrix} + (1 - \alpha) k_{inclusion}$ with $\alpha = 0.55$ for $C < 1$ and $\alpha = 0.45$ for $C > 1$.
<figure><img src="image_12.png" /><figcaption>**Fig. 11.** An isotropic spherical inclusion embedded into an isotropic matrix: number of iterations at convergence (tolerance $\epsilon_L^{2} = 10^{-10}$) for different contrasts between inclusion and matrix bulk moduli. Results obtained with the basic scheme associated with the tetrahedral stencil, the rotated discretization and the Moulinec-Suquet (MS) discretization have been obtained with a reference medium that has the same Poisson coefficient as the inclusion and matrix and a bulk modulus given by $k_0 = \alpha k_{\text{matrix}} + (1-\alpha) k_{\text{inclusion}}$ with $\alpha = 0.55$ for $C < 1$ and $\alpha = 0.45$ for $C > 1$. Results for the MS method accelerated by the augmented Lagrangian scheme, Eyre-Milton scheme and the Monchiet-Bonnet polarization method have been taken from Ref. [12].</figcaption></figure>

The main conclusion is that, whereas the number of iterates increases linearly with $C$ (or with $1/C$ if $C < 1$) for the Moulinec-Suquet discretization combined with the basic scheme and displays square-root dependencies for the three accelerated iterative schemes, the tetrahedral stencil together with the rotated discretization, for which the number of iterates is however slightly higher, converge rapidly and, most importantly, require a number of iterates that essentially does not depend on the contrast. Therefore, microstructures with voids or infinitely rigid inclusions can be handled as easily as situations where the contrast is finite, even though the iterative scheme used here follows the simple basic scheme.

### 3.5. A 3-phases 3D medium with double contrast

We consider now a microstructure that contains simultaneously soft and hard inclusions. Our purpose is to test the ability of the tetrahedral stencil to treat this “double contrast” situation which, as discussed recently by Sab et al. [44], is particularly challenging, both in terms of discretization scheme and iterative solver, when the contrast between the hard and soft inclusions is large and eventually infinite.

The periodic medium consists of $64^3$ cells. We insert a hard spherical inclusion at the centre of the simulation box and soft spherical inclusions centred on its height vertices. The radius of the spheres is, in cell units, $R = 16$. The matrix, the hard and the soft inclusions are elastically isotropic and have the same Poisson’s ratio $\nu = 0.3$ but different shear moduli given by $\mu_1, \chi\mu_1$ and $\mu_1/\chi$, respectively, where $\chi$ denotes the double contrast ratio between the hard inclusions and the matrix and between the matrix and the soft inclusions.

In the situations analysed in the previous paragraph, where single contrasts were considered, we observed that, for the tetrahedral stencil associated to the basic scheme, equilibrium with a precision $\epsilon_L^2$ equal to $10^{-10}$ was achieved with a moderate number of iterations (around 300) that essentially does not depend on the contrast up to the smallest ($10^{-5}$) and largest ($10^5$) contrasts considered. In the present situation, we observed that, with the same precision, the number of iterates $N$ to reach equilibrium, which is still relatively moderate for a double contrast $\chi$ equal to $10$ ($N = 600$), becomes prohibitively too large (larger than $10^4$) for $\chi \ge 100$. Therefore, microstructures with large double contrasts require the use of a faster solver, such as the Adaptive Eyre-Milton (AEM) scheme recently proposed by Sab et al. [44] which permits to consider easily large and even infinite double contrasts.

Here, we confine our numerical investigations to large but finite double contrast $\chi$. We therefore use the polarization method, whose adaptation to the tetrahedral discretization has been presented in Section 2.7, Algorithm 3. We set $\gamma = 0$, which means that we use the Eyre-Milton method. The reference stiffness tensor $J_s^0$ is isotropic and, as recommended [11], its bulk and shear moduli have been set to the geometric means of the hard and soft particle moduli, which here amounts to those of the matrix. Convergence to equilibrium was monitored with the criterion $\epsilon_{pol}$ defined at the end of Section 2.8 and iterations were stopped when $\epsilon_{pol}$ reaches $10^{-10}$.

We report in Fig. 12 results for a double contrast $\chi = 10^5$ when the overall system is subjected to the uniaxial unit strain $\bar{\epsilon}_{ij} = \delta_{ij} \delta_{ji}$. The von Mises stress map obtained with the tetrahedral stencil is compared to those obtained with the rotated scheme and the standard staggered grid [32,33]. We observed typical checkerboard effects generated by the rotated stencil within the hard
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** A 3-phase 3D microstructure with double contrast $\chi = 10^5$ subjected to a uniaxial strain $\bar{\epsilon}_{11} = 1$: von Mises stress maps along the (110) plane that goes through the centre of computational box; comparison between the rotated discretization (left), the standard staggered grid (middle) and the tetrahedral stencil (right) (units: matrix shear modulus $\mu_1$). The solver used here is the polarization-based method shown in Algorithm 3 with $\gamma = 0$. The computational box consists of $64^3$ voxels.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 13.** A 3D microstructure with a double contrast $\chi = 10^5$ subjected to a uniaxial strain $\bar{\epsilon}_{11} = 1$: errors to equilibrium $\epsilon_{L2}$ and $\epsilon_{pol}$ as a function of iterates with the tetrahedral stencil and the Eyre-Milton polarization-based method.</figcaption></figure>

spherical inclusion and in the neighbourhood of its interface with the matrix. The solution corresponding to the staggered grid is devoid of oscillations but is non-symmetric with respect to the vertical and horizontal middle lines, albeit the microstructure does obey these symmetries. This artefact results from the staggered placement of the deformation components assigned to a given cubic voxel on four sites that do not all belong to the grid to which the corresponding stiffness tensor is assigned, with the result that the deformation components are not mathematically bound to follow any local symmetry displayed by the stiffness tensor or by the geometry of the microstructure. As discussed in [33], the remedy to correct for this artefact is to consider a doubly-fine grid that involves height sub-voxels for each initial cubic voxel, at the expense of the computational effort. Finally, we note that the solution corresponding to the tetrahedral stencil is exempt of the artefacts displayed by the two other stencils.

We close this section by a brief analysis of the convergence rate of the polarization-based method adapted to the tetrahedral stencil. We show in Fig. 13 the evolution of the errors to equilibrium as a function of the iterate for the double contrast $\chi = 10^5$. The error $\epsilon_{pol}$, which is based on the $\lambda^0$-norm of the residuals (see (2.81)) is compared to the error $\epsilon_{L2}$, which is based on the $L^2$-norm of the stress divergence (see (2.71)). We note that, after a finite transient regime, the two errors appear to be almost identical and
<caption>**Table 1**</caption>

3D medium with double contrast : convergence rate of the polarization-based algorithm associated with the tetrahedral stencil (Algorithm 3 with $\gamma = 0$); number of iterates as a function of the double contrast $\chi$ to reach equilibrium with the precision $\epsilon_{pol} = 10^{-10}$.

<table><thead><tr><th>χ</th><th>10</th><th>10<sup>2</sup></th><th>10<sup>3</sup></th><th>10<sup>4</sup></th><th>10<sup>5</sup></th><th>10<sup>6</sup></th></tr></thead><tbody><tr><td>Iterates</td><td>67</td><td>237</td><td>339</td><td>386</td><td>425</td><td>462</td></tr></tbody></table>

<figure><img src="image_15.png" /><figcaption>**Fig. 14.** 3D microstructure with randomly distributed non-overlapping spherical voids : von Mises stress maps along a (001) plane; comparison between the results obtained with the tetrahedral stencil (left panel) and with the staggered grid discretization (right panel) (units: MPa). The computational box consists of 128³ voxels.</figcaption></figure>

display both a linear convergence. Next, we report in Table 1 the number of iterates needed to reach equilibrium as a function of the double contrast $\chi$ with a tolerance $\epsilon_{pol}$ set equal to $10^{-10}$. Strain fields $\epsilon^{(T_1)}$ and $\epsilon^{(T_2)}$ are initially set to the applied strain $\bar{\epsilon}_{ij} = \delta_{ij}\delta_{j1}$. We observe that, up to the high double contrast $\chi = 10^6$, the number of iterates increases slowly with $\chi$ and stays moderate, which, in light of the simplicity of Algorithm 3, is satisfactory.

### 3.6. A 3D microstructure with non-overlapping voids

We consider here a more realistic heterogeneous materials that consists of non-overlapping spherical voids randomly distributed into an elastically isotropic matrix. The periodic medium is discretized into 128³ cells, the volume fraction of the void phase is 30% and the radius of the voids is, in cell units, $R = 13$. The matrix Young's modulus and Poisson coefficient are given by $E = 66.9$ GPa and $\nu = 0.25$ and the overall system is subjected to the uniaxial external strain $\bar{\epsilon}_{22} = 0.05$.

First, we compare in Fig. 14 von Mises maps obtained with the tetrahedral and staggered grid discretizations. We observe that the two maps are similar and devoid of any unphysical oscillations. The two stress fields differ in the close neighbourhood of the interfaces between voids and matrix. A close examination shows that these differences come from the artefacts in the staggered grid map which are generated by the non-symmetrical treatment of interfaces discussed in the previous example.

The stress field $\sigma$ that results from the tetrahedral discretization is the arithmetic average of the stress fields $\sigma^{(T_1)}$ and $\sigma^{(T_2)}$ (see Eq. (2.49)). We report in Fig. 15 von Mises stress maps of $\sigma^{(T_1)}$ and $\sigma^{(T_2)}$ (top row) together with the von Mises maps of their arithmetic average $\sigma$ and (half) difference $(\sigma^{(T_1)} - \sigma^{(T_2)})/2$ (bottom row). We observe that the stress fields $\sigma^{(T_1)}$, $\sigma^{(T_2)}$ and $\sigma$ are almost identical, except in the immediate neighbourhood of the interfaces between voids and matrix. This results from the FCC nature of the tetrahedral stencil. In fact, the geometry of the initial cubic grid is recovered only if we consider the arithmetic average of the two FCC discretizations. Indeed, looking at Fig. 1 and considering a specific strain site, say a strain site of type A, we observe that the action of the displacement field on this strain will be consistent with the full simple cubic geometry only if we consider simultaneously the FCC displacement subgrids 1 and 2, which, using the notations defined in Section 2.1, generate the stress components $\sigma^{(A1)}$ and $\sigma^{(A2)}$, respectively. Therefore, we must consider their arithmetic average $(\sigma^{(A1)} + \sigma^{(A2)})/2$. The same applies to a strain site of type B. These are the reasons why only the average $(\sigma^{(T_1)} + \sigma^{(T_2)})/2$ is meaningful, as it realizes simultaneously the two averages $(\sigma^{(A1)} + \sigma^{(A2)})/2$ and $(\sigma^{(B1)} + \sigma^{(B2)})/2$, as can be seen from Eqs. (2.46) and (2.49). Of course, the difference between $\sigma^{(T_1)}$ and $\sigma^{(T_2)}$ is accentuated in regions where the stress is highly inhomogeneous at the scale of the voxel length, such as close to interfaces, but this difference only reflects the underlying local heterogeneity and has no other significance.

Next, we show in Fig. 16 the evolution of the errors to equilibrium as a function of the iterate index for the tetrahedral and staggered grid discretizations within the basic scheme. For comparison, we report the error $\epsilon_{L^2}$, defined in Eq. (2.71), together
<figure><img src="image_16.png" /><figcaption>Fig. 15. 3D microstructure with randomly distributed non-overlapping spherical voids: von Mises stress maps along a (100) plane with the tetrahedral stencil; $\sigma(T_1)$ (top left), $\sigma(T_2)$ (top right), $\frac{1}{2}(\sigma(T_1) + \sigma(T_2))$ (bottom left), $\frac{1}{2}(\sigma(T_1) - \sigma(T_2))$ (bottom right) (units: MPa).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 16. 3D microstructure with randomly distributed non-overlapping spherical voids: errors to equilibrium $\epsilon_L$ and $\epsilon_{H_{per}^{-1}}$ as a function of iterates with the basic scheme; comparison between the tetrahedral stencil (red) and the standard staggered grid discretization (green). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Fig. 17. (a) 3D microstructure with randomly distributed spherical voids with two different radii and N = 256 (see text). (b) Convergence of the effective Young's modulus with respect to the number N³ of voxels with the tetrahedral, rotated and standard staggered grid methods (error bars reflect the dispersion of the Young's modulus measures along the three cartesian axis); in black, Hashin-Shtrikman theoretical upper bounds. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<table><caption>**Table 2**

Tensile test under stress control : equilibrium average strain components (×10³) with applied stress $\sigma_{ij}^a = \sigma_0\delta_{il}\delta_{j1}$ and N = 256 voxels along each cartesian axis ($\sigma_0$ = 1 MPa).</caption><thead><tr><th>ε̃<sub>11</sub></th><th>ε̃<sub>22</sub></th><th>ε̃<sub>33</sub></th><th>ε̃<sub>23</sub></th><th>ε̃<sub>13</sub></th><th>ε̃<sub>12</sub></th></tr></thead><tbody><tr><td>1.484</td><td>-0.462</td><td>-0.460</td><td>-0.001</td><td>0.003</td><td>0.009</td></tr></tbody></table>

with the error $\epsilon_{H_{per}^{-1}}$, defined in Eq. (2.77), which are based respectively on the L²-norm and $H_{per}^{-1}$-norm of the (discrete equivalent) divergence of the stress. Errors $\epsilon_{L^2}$ and $\epsilon_{H_{per}^{-1}}$ have been numerically computed according to (2.72) and (2.78), respectively. First, we note that, for each discretization scheme, $\epsilon_{L^2}$ and $\epsilon_{H_{per}^{-1}}$ are close to each other and therefore display identical evolutions. Second, we observe that, for the two discretization stencils, the convergence to equilibrium is linear, with a sensible advantage for the tetrahedral discretization as compared to the staggered one.

## 3.7. Homogenization and convergence with grid refinement

For our last numerical experiments, we present a brief analysis of the convergence of the effective properties of a heterogeneous microstructure as a function of the refinement of the grid used to discretize a representative cell. The material consists of a random repartition of two non-overlapping spherical void assemblies that contain respectively 155 and 63 voids of radius R = 8N/128 and R = 10N/128 (in voxel units) where N is the number of voxels along each cartesian axis (see Fig. 17(a) for N = 256). The volume fraction of the void phase is equal to 0.283 for N = 256, with slight variations with N due to the voxellation of the interfaces of the spherical voids. The matrix is elastically isotropic with a Young's modulus and a Poisson coefficient given by E = 1.25 GPa and ν = 0.35. We perform a simple tensile experiment under stress control with a tensile force applied along one of the cartesian axis, e.g. $\sigma_{ij}^a = \sigma_0\delta_{il}\delta_{j1}$. The overall effective Young's modulus $E^{eff}$ is computed through the ratio $\sigma_0/\bar{\epsilon}_{il}$ where $\bar{\epsilon}$ is the obtained average strain and i the direction of the applied tensile force.

If the representative volume is sufficiently large, we expect the microstructure to be elastically isotropic. This is confirmed in Table 2 where we show, as an example, the calculated average strain obtained in one of our simulations. We observe that the results are consistent, within small statistical errors, with what we would expect from a stress control tensile experiment on an isotropic material, i.e. equal contractions along the directions perpendicular to the load and zero shear strains.

We present in Fig. 17(b) the evolution of $E^{eff}$ as a function of the number N³ of D.O.F. used to discretize the microstructure. Results obtained with the tetrahedral stencil are compared to those obtained with the rotated and standard staggered grid methods. We also report the theoretical upper bound on the effective Young's modulus that results from the Hashin-Shtrikman variational inequalities [49,50]. Small variations of this upper bound with the number of voxels are due to fluctuations in the void volume generated by the voxellation of the interfaces, while the theoretical lower bound vanishes because the inclusion phase is made up of voids. First, we observe that, as expected, the computed effective Young's modulus are all below the upper bounds. Next, we note that, in contrast to the rotated and staggered grid methods, the tetrahedral stencil leads to an effective Young's modulus that decreases monotonically with increasing the grid refinement. Finally, we note that the tetrahedral stencil predicts a (slightly) stiffer behaviour than the other two methods.
## 4. Conclusions

We have introduced a new FFT-based method for the computation of the macroscopic mechanical properties of heterogeneous microstructures. The discretization scheme relies on a tetrahedral stencil. Its formulation is based on the division of the initial strain and displacement simple cubic grids into face-centred cubic subgrids. The two resulting minimization problems can be reduced to a single gradient-based iterative scheme that operates on the full initial displacement grid and that can be interpreted as the adaptation of the Moulinec-Suquet basic scheme to the tetrahedral discretization. This reduction is made possible by the introduction of stress and strain fields, called $T_1$ and $T_2$ fields, defined on the full initial strain grid and that appropriately mix the two discretizations, and to the definition of a judiciously chosen scalar product that operates on $T_1$ and $T_2$ strain fields. Implementation of the corresponding algorithm is simple and can operate on the displacement field, as it is derived and which is preferred, or on the strain fields. Also, an implementation within the general polarization-based method of Monchiet-Bonnet has been proposed, which, if needed, can also be easily implemented.

To summarize, the main characters and advantages of the method are the following:

* the discretization scheme, which is based on a tetrahedral stencil, has been chosen in such a way that mechanical equilibrium is mathematically stable, i.e. automatically free from any unphysical oscillations or checkerboard effects, such as those observed with the original Moulinec-Suquet discretization or the discrete rotated stencil;

* all the components (diagonal and non-diagonal) of a tensorial quantity live on the same location, which avoids the numerical artefacts that we may observe when, as within the standard staggered grid method, the components of a given tensor are staggered on different grids, which may generate non-symmetrical treatments of sharp interfaces; it avoids also the difficulty that arises when using a constitutive law that links together tensor components that do not sit on the same location;

* convergence to equilibrium with the fixed point algorithm of the original basic scheme is fast and requires a number of iterates that, for large contrasts, is essentially independent on the contrast, so that infinite contrast is obtained without additional computational cost;

* if needed, an implementation within the more efficient polarization-based fixed point algorithm of Monchiet-Bonnet is proposed.

* convergence of effective elastic properties with grid refinement is monotonic and fast.

To conclude, the focus in this work has been on providing a discretization scheme that excludes unphysical responses in mechanical fields. The proposed stencil may be however easily transposed to other fields, such as electrokinetics or multiferroics. Lastly, in order to better discretize interfaces that do not conform to a regular grid, an interesting improvement could be to use the proposed stencil on non-uniform grids with a mesh adaptive strategy, such as recently proposed in [51,52].

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgements

The author thanks Yann Le Bouar, Antoine Ruffini and Mathieu Fèvre for numerous stimulating discussions and the three anonymous reviewers whose comments and suggestions helped improve and clarify the manuscript.

## Appendix A. Action of the finite difference operators $D_i^{(Kk)}$

We illustrate the use of the operators $D_i^{(Kk)}$, defined in Eq. (2.5), to obtain the explicit expressions of the fluctuating strain fields defined in Eq. (2.4).

Consider for example the fluctuating strain components $\delta\epsilon_{11}^{(A1)}$ and $\delta\epsilon_{12}^{(A1)}$ that sit on site (000). Application of the operator $\underline{D}^{A1}$, defined in (2.5a), leads to

$$ \delta\epsilon_{11}^{(A1)}(000) = \frac{1}{2} \left\{ u_1 \left(000 + \frac{111}{2}\right) + u_1 \left(000 + \frac{\overline{111}}{2}\right) - u_1 \left(000 + \frac{\overline{111}}{2}\right) - u_1 \left(000 + \frac{\overline{111}}{2}\right) \right\}, $$

and

$$ \delta\epsilon_{12}^{(A1)}(000) = \frac{1}{4} \left\{ \begin{aligned} & u_2 \left(000 + \frac{111}{2}\right) + u_2 \left(000 + \frac{\overline{111}}{2}\right) - u_2 \left(000 + \frac{\overline{111}}{2}\right) - u_2 \left(000 + \frac{\overline{111}}{2}\right) \\ & + \frac{1}{4} \left\{ u_1 \left(000 + \frac{111}{2}\right) + u_1 \left(000 + \frac{\overline{111}}{2}\right) - u_1 \left(000 + \frac{\overline{111}}{2}\right) - u_1 \left(000 + \frac{\overline{111}}{2}\right) \right\} \end{aligned} \right\}. $$

These expressions correspond to the tetrahedral difference scheme defined geometrically by the blue lines that emerge from site A in Fig. 1. Consider now the strain components $\delta\epsilon_{11}^{(B1)}$ and $\delta\epsilon_{12}^{(B1)}$ that sit on a site that belongs to strain subgrid $Y_B$, say site (100). According to the definition of the operator $\underline{D}^{B1}$ given in (2.5d), they are obtained by applying the negative of the central symmetry $C_{(100)}$ around site (100) to the two expressions given above, in which the site (000) has been replaced by the site (100). This leads to
$$
\begin{align*}
\delta\epsilon_{11}^{(B1)}(100) &= -\frac{1}{2} C_{(100)} \left\{ u_1 \left(100 + \frac{111}{2}\right) + u_1 \left(100 + \frac{111}{2}\right) - u_1 \left(100 + \frac{111}{2}\right) - u_1 \left(100 + \frac{111}{2}\right) \right\}, \\
&= -\frac{1}{2} \left\{ u_1 \left(100 + \frac{111}{2}\right) + u_1 \left(100 + \frac{111}{2}\right) - u_1 \left(100 + \frac{111}{2}\right) - u_1 \left(100 + \frac{111}{2}\right) \right\}, \\
&= \frac{1}{2} \left\{ u_1 \left(100 + \frac{111}{2}\right) + u_1 \left(100 + \frac{111}{2}\right) - u_1 \left(100 + \frac{111}{2}\right) - u_1 \left(100 + \frac{111}{2}\right) \right\},
\end{align*}
$$

and

$$
\begin{align*}
\delta\epsilon_{12}^{(B1)}(100) &= -\frac{1}{4} C_{(100)} \left\{ u_2 \left(100 + \frac{111}{2}\right) + u_2 \left(100 + \frac{111}{2}\right) - u_2 \left(100 + \frac{111}{2}\right) - u_2 \left(100 + \frac{111}{2}\right) \right\} \\
&\quad -\frac{1}{4} C_{(100)} \left\{ u_1 \left(100 + \frac{111}{2}\right) + u_1 \left(100 + \frac{111}{2}\right) - u_1 \left(100 + \frac{111}{2}\right) - u_1 \left(100 + \frac{111}{2}\right) \right\}, \\
&= -\frac{1}{4} \left\{ u_2 \left(100 + \frac{111}{2}\right) + u_2 \left(100 + \frac{111}{2}\right) - u_2 \left(100 + \frac{111}{2}\right) - u_2 \left(100 + \frac{111}{2}\right) \right\} \\
&\quad -\frac{1}{4} \left\{ u_1 \left(100 + \frac{111}{2}\right) + u_1 \left(100 + \frac{111}{2}\right) - u_1 \left(100 + \frac{111}{2}\right) - u_1 \left(100 + \frac{111}{2}\right) \right\}, \\
&= \frac{1}{4} \left\{ u_2 \left(100 + \frac{111}{2}\right) + u_2 \left(100 + \frac{111}{2}\right) - u_2 \left(100 + \frac{111}{2}\right) - u_2 \left(100 + \frac{111}{2}\right) \right\} \\
&\quad +\frac{1}{4} \left\{ u_1 \left(100 + \frac{111}{2}\right) + u_1 \left(100 + \frac{111}{2}\right) - u_1 \left(100 + \frac{111}{2}\right) - u_1 \left(100 + \frac{111}{2}\right) \right\}.
\end{align*}
$$

These expressions correspond to the tetrahedral difference scheme defined geometrically by the red lines that emerge from site B in Fig. 1.

The above methodology generalizes systematically to all strain components, with the result that the strains $\delta\epsilon^{(A1)}$ and $\delta\epsilon^{(B2)}$ are associated to the blue tetrahedra that emerge from sites A and B, whereas the strains $\delta\epsilon^{(A1)}$ and $\delta\epsilon^{(A2)}$ are associated to the red ones. These blue and red tetrahedra are referred to as tetrahedra of type $T_1$ and $T_2$, respectively.

Appendix B. Definition of the discrete divergence of the stress $\boldsymbol{\sigma} = (\boldsymbol{\sigma}^{(T_1)} + \boldsymbol{\sigma}^{(T_2)})/2$

The finite difference operators $D^{(Kk)}$ permit also to define the discrete divergences of the stress fields $\sigma^{(Kk)}$, which, for a given $k$, live on the displacement subgrid $Y_k$. Examination of Fig. 1, shows immediately that they are given by, for $i = 1, 2, 3$,

$$
\begin{align}
\forall \underline{r} \in Y_1, \quad (\operatorname{div} \boldsymbol{\sigma}^{(A1)})_i(\underline{r}) &= -\sum_j C_{\underline{r}} \{D_j[\boldsymbol{\sigma}_{ij}^{(A1)}](\underline{r})\}, \nonumber \\
\forall \underline{r} \in Y_2, \quad (\operatorname{div} \boldsymbol{\sigma}^{(B2)})_i(\underline{r}) &= -\sum_j C_{\underline{r}} \{D_j[\boldsymbol{\sigma}_{ij}^{(B2)}](\underline{r})\}, \tag{B.1} \\
\forall \underline{r} \in Y_2, \quad (\operatorname{div} \boldsymbol{\sigma}^{(A2)})_i(\underline{r}) &= \sum_j \underline{D}_j[\boldsymbol{\sigma}_{ij}^{(A2)}](\underline{r}), \nonumber \\
\forall \underline{r} \in Y_1, \quad (\operatorname{div} \boldsymbol{\sigma}^{(B1)})_i(\underline{r}) &= \sum_j \underline{D}_j[\boldsymbol{\sigma}_{ij}^{(B1)}](\underline{r}), \nonumber
\end{align}
$$

where the operator $C_{\underline{r}}$, already introduced in 2.1, is the central symmetry around site $\underline{r}$. In Fourier space, one gets

$$
\begin{align}
\widehat{\operatorname{div} \boldsymbol{\sigma}^{(A1)}}(\underline{q}) &= -\hat{\boldsymbol{\sigma}}^{(A1)}(\underline{q}) \cdot \underline{\underline{D}}(\underline{q}), \nonumber \\
\widehat{\operatorname{div} \boldsymbol{\sigma}^{(B2)}}(\underline{q}) &= -\hat{\boldsymbol{\sigma}}^{(B2)}(\underline{q}) \cdot \underline{\underline{D}}(\underline{q}), \nonumber \\
\widehat{\operatorname{div} \boldsymbol{\sigma}^{(A2)}}(\underline{q}) &= \hat{\boldsymbol{\sigma}}^{(A2)}(\underline{q}) \cdot \underline{\underline{D}}(\underline{q}), \nonumber \\
\widehat{\operatorname{div} \boldsymbol{\sigma}^{(B1)}}(\underline{q}) &= \hat{\boldsymbol{\sigma}}^{(B1)}(\underline{q}) \cdot \underline{\underline{D}}(\underline{q}), \tag{B.2}
\end{align}
$$

which leads to the following expression for the divergences of the combined stress fields $\sigma^{T_1}$ and $\sigma^{T_2}$ defined in (2.26):

$$
\widehat{\operatorname{div} \sigma^{(T_1)}}(\underline{q}) = -\hat{\sigma}^{(T_1)}(\underline{q}) \cdot \underline{\underline{D}}(\underline{q}), \qquad (B.3)
$$

$$
\widehat{\operatorname{div} \sigma^{(T_2)}}(\underline{q}) = \hat{\sigma}^{(T_2)}(\underline{q}) \cdot \underline{\underline{D}}(\underline{q}). \qquad (B.4)
$$

We define the divergence of the average stress $\boldsymbol{\sigma} = (\boldsymbol{\sigma}^{(T_1)} + \boldsymbol{\sigma}^{(T_2)})/2$ through

$$
\widehat{\operatorname{div} \boldsymbol{\sigma}}(\underline{q}) = \frac{1}{2} \left( \widehat{\operatorname{div} \boldsymbol{\sigma}^{(T_1)}}(\underline{q}) + \widehat{\operatorname{div} \boldsymbol{\sigma}^{(T_2)}}(\underline{q}) \right), \qquad (B.5)
$$
where it must be kept in mind that each of the two divergence operators that appear on the r.h.s has a mathematical expression that is specific to the field to which it applies, as seen in (B.3) and (B.4), which in turn means that this equation should be interpreted as a definition of the divergence operator when it applies to the average stress field $\sigma$. Finally, one gets

$$ \widehat{\mathrm{div}} \sigma(\underline{q}) = \frac{1}{2} \left( -\hat{\sigma}^{(T_1)}(\underline{q}) \cdot \overline{D}(\underline{q}) + \hat{\sigma}^{(T_2)}(\underline{q}) \cdot \underline{D}(\underline{q}) \right). \quad (B.6) $$

Therefore, the derivative of the elastic energy, given in (2.25), may be written as

$$ \widehat{\mathrm{div}} \sigma(\underline{q}) = - \frac{\partial E}{\partial \hat{u}(\underline{q})}, \quad (B.7) $$

which, in real space, leads to

$$ \mathrm{div} \sigma(\underline{r}) = -N \frac{\partial E}{\partial \hat{u}(\underline{r})}. \quad (B.8) $$

## Appendix C. Derivative of discrete elastic energy $E^{(k)}$ with respect to Fourier component $\underline{u}^{(k)}(\underline{q})$ of the displacement field.

In this section, we demonstrate the relationship given in Eq. (2.23). This is done simply as follows. Using the links between the Fourier components of the fluctuating strains $\hat{u}_s^{(Kk)}(\underline{q})$ and the displacement fields $\hat{u}^{(k)}(\underline{q})$ given in Eq. (2.12), we have

$$ \frac{\partial E^{(k)}}{\partial \hat{u}_i^{(k)}(\underline{q})} = \sum_{K=A,B} \frac{\partial E^{(k)}}{\partial \hat{\epsilon}_{ij}^{(Kk)}(\underline{q})} D_j^{(Kk)}(\underline{q}). \quad (C.1) $$

Using the inverse Fourier transform (2.10b) applied to $\hat{\epsilon}^{(Kk)}(\underline{q})$ and the chain rule, the first term in the right-hand side may be written as

$$ \frac{\partial E^{(k)}}{\partial \hat{\epsilon}_{ij}^{(Kk)}(\underline{q})} = \sum_{\underline{r} \in K} \frac{\partial E^{(k)}}{\partial \hat{\epsilon}_{ij}^{(Kk)}(\underline{r})} \exp(i\underline{q} \cdot \underline{r}). \quad (C.2) $$

On the other hand, the stress $\sigma_{ij}^{(Kk)}(\underline{r})$ is, by definition, related to the total energy $E^{(k)}$ given in Eq. (2.16) by

$$ \sigma_{ij}^{(Kk)}(\underline{r}) = N \frac{\partial E^{(k)}}{\partial \hat{\epsilon}_{ij}^{(Kk)}(\underline{r})}. \quad (C.3) $$

Altogether, Eqs. (C.1)–(C.3) lead to

$$ \frac{\partial E^{(k)}}{\partial \hat{u}_i^{(k)}(\underline{q})} = \sum_{K=A,B} \hat{\sigma}_{ij}^{(Kk)}(\underline{q}) D_j^{(Kk)}(\underline{q}), \quad (C.4) $$

where $\hat{\sigma}^{(Kk)}(\underline{q})$ is the Fourier transform of the stress $\sigma^{(Kk)}(\underline{r})$.

## Appendix D. About the non-analyticity of the strain Green function within the rotated scheme

Within the rotated scheme, the displacement and, consequently, strain Green functions are undefined at some finite q-points in Fourier space. These singularities generate mathematical artefacts in the real space strain and stress fields that cannot be suppressed by any recipe, as we explain here.

The rotated scheme may be summarized as follows. The discretization scheme involves a SC grid for the strain field and a staggered SC grid for the displacement field. The Fourier transforms of the finite difference operators, that link the fluctuating strain and displacement fields through $\delta\hat{\epsilon}(\underline{q}) = \underline{D}(\underline{q}) \otimes_s \hat{u}(\underline{q})$, are given by,<fnref n="2" />

$$ \begin{align} D_1(\underline{q}) &= \frac{1}{4} \{ \exp(iq_1) - 1 \} \{ \exp(iq_2) + 1 \} \{ \exp(iq_3) + 1 \}, \nonumber \\ D_2(\underline{q}) &= \frac{1}{4} \{ \exp(iq_1) + 1 \} \{ \exp(iq_2) - 1 \} \{ \exp(iq_3) + 1 \}, \nonumber \\ D_3(\underline{q}) &= \frac{1}{4} \{ \exp(iq_1) + 1 \} \{ \exp(iq_2) + 1 \} \{ \exp(iq_3) - 1 \}. \tag{D.1} \end{align} $$

The displacement-based Lippmann–Schwinger equation reads

$$ \forall \underline{q} \neq (000)^*, \quad \hat{G}^0(\underline{q})^{-1} \hat{u}(\underline{q}) = -\hat{\tau}(\underline{q}) \cdot \overline{D}(\underline{q}), \quad (D.2) $$

where $\hat{\tau}(\underline{q}) = \hat{\sigma}(\underline{q}) - \lambda^0 : \hat{\tau}(\underline{q})$ is the polarization tensor and $\hat{\sigma}(\underline{q})$ the Fourier transform of the stress field which, in a small deformation setting, is linked to the strain $\epsilon(\underline{r}) = \bar{\epsilon} + \delta\epsilon(\underline{r})$ by $\sigma(\underline{r}) = \hat{\lambda}(\underline{r}) : (\epsilon(\underline{r}) - \epsilon^0(\underline{r}))$. $\bar{\epsilon}$ is the spatial average of the strain and the field $\epsilon^0(\underline{r})$

<footnote n="2">² If the displacement and strain fields are referred to the same origin in real space, as in Fig. 1 each component $D_i(\underline{q})$ should be multiplied by $\exp(-i(q_1 + q_2 + q_3)/2)$. However, this prefactor has no impact on the computed strain field, since, in the inverse of the displacement Green function $\hat{G}^0(\underline{q})$ given in (D.3) and in the strain Green function $\hat{\tau}^0(\underline{q})$ introduced below in (D.7) and defined in (D.8), the components $D_i(\underline{q})$ appear always by product of the form $D_i(\underline{q})D_j(\underline{q})$.</footnote>
is a local eigenstrain. The inverse $\hat{\mathbf{G}}^0(\underline{q})^{-1}$ of the displacement Green function is given by

$$ \hat{\mathbf{G}}^0(\underline{q})^{-1} = \underline{D}(\underline{q}) \lambda_{\underline{q}}^0 \underline{D}(\underline{q}). \qquad (D.3) $$

As noted in [34], if the numbers of voxels along any two cartesian axis, say i and j, are even, all the components of the q-space differential operator $\underline{D}(\underline{q})$ are equal to zero when components $q_i$ and $q_j$ of the q-vector reach the highest frequency, i.e. when $q_i = q_j = \pi$, which makes the Green function, and therefore the Fourier transforms of displacement, strain and stress fields, undefined at the corresponding q-points. In fact, this non-definitiveness is one of the manifestation of the singular behaviour displayed by the Green function in the neighbourhoods of the three lines that emerge from the point $(\pi\pi\pi)^*$ and that are aligned along the directions $(100)^*$, $(010)^*$ and $(001)^*$, to which we refer below as the lines $L_{100}$, $L_{010}$ and $L_{001}$, respectively.

To illustrate this behaviour, we consider the neighbourhood of the line $L_{100}$. We note $\underline{q} = (q_1, \pi + \epsilon \sin\theta, \pi + \epsilon \cos\theta)^*$ a point within this neighbourhood, with $\epsilon \ge 0$ and $\theta \in [0, 2\pi]$. Note that we introduce continuous notations for $\underline{q}$, as required to analyse analytical properties in q-space. For any fixed $q_1 \ne \pi$ and in the limit $\epsilon \to 0$, the differential operators $D_i(\underline{q})$ behave, to the lowest order in $\epsilon$, as

$$ \begin{aligned} D_1(q_1, \pi + \epsilon \sin\theta, \pi + \epsilon \cos\theta) &\sim \frac{\epsilon^2}{2} \cos\left(\frac{q_1 - \pi}{2}\right) \exp\left(i \frac{q_1 - \pi}{2}\right) \cos\theta \sin\theta, \\ D_2(q_1, \pi + \epsilon \sin\theta, \pi + \epsilon \cos\theta) &\sim \epsilon \sin\left(\frac{q_1 - \pi}{2}\right) \exp\left(i \frac{q_1 - \pi}{2}\right) \cos\theta, \\ D_3(q_1, \pi + \epsilon \sin\theta, \pi + \epsilon \cos\theta) &\sim \epsilon \sin\left(\frac{q_1 - \pi}{2}\right) \exp\left(i \frac{q_1 - \pi}{2}\right) \sin\theta. \end{aligned} \qquad (D.4) $$

This limiting behaviour is reflected in the inverse Green operator $\hat{\mathbf{G}}^0(\underline{q})^{-1}$ given in Eq. (D.3) which, in the limit $0 \le \epsilon \ll |\tan \frac{q_1 - \pi}{2}|$ and for $q_1 \ne \pi$, behaves as

$$ \times \begin{pmatrix} C_{44} \sin^2 \frac{q_1 - \pi}{2} & \frac{\epsilon}{4} (C_{12} + C_{44}) \sin(q_1 - \pi) \cos^2 \theta \sin \theta & \frac{\epsilon}{4} (C_{12} + C_{44}) \sin(q_1 - \pi) \cos \theta \sin^2 \theta \\ \frac{\epsilon}{4} (C_{12} + C_{44}) \sin(q_1 - \pi) \cos^2 \theta \sin \theta & (C_{11} \sin^2 \theta + C_{44} \sin^2 \theta) \sin^2 \frac{q_1 - \pi}{2} & (C_{12} + C_{44}) \sin \theta \cos \theta \sin^2 \frac{q_1 - \pi}{2} \\ \frac{\epsilon}{4} (C_{12} + C_{44}) \sin(q_1 - \pi) \cos \theta \sin^2 \theta & (C_{12} + C_{44}) \sin \theta \cos \theta \sin^2 \frac{q_1 - \pi}{2} & (C_{11} \sin^2 \theta + C_{44} \cos^2 \theta) \sin^2 \frac{q_1 - \pi}{2} \end{pmatrix}, \qquad (D.5) $$

where we have considered a cubic elastic anisotropy. As a result, in the limit $0 \le \epsilon \ll |\tan \frac{q_1 - \pi}{2}|$ and for $q_1 \ne \pi$, the Green function $\hat{\mathbf{G}}^0(\underline{q})$ behaves as

$$ \times \begin{pmatrix} \omega_{11}(\theta) & \epsilon \sin^{-2} \frac{q_1 - \pi}{2} \sin(q_1 - \pi)\omega_{12}(\theta) & \epsilon \sin^{-2} \frac{q_1 - \pi}{2} \sin(q_1 - \pi)\omega_{13}(\theta) \\ \epsilon \sin^{-2} \frac{q_1 - \pi}{2} \sin(q_1 - \pi)\omega_{12}(\theta) & \omega_{22}(\theta) & \omega_{23}(\theta) \\ \epsilon \sin^{-2} \frac{q_1 - \pi}{2} \sin(q_1 - \pi)\omega_{13}(\theta) & \omega_{23}(\theta) & \omega_{33}(\theta) \end{pmatrix}. \qquad (D.6) $$

The explicit $\theta$-dependence of the functions $\omega_{ij}(\theta)$ are irrelevant for the rest of the discussion.

The impact of the limiting behaviour of $\hat{\mathbf{G}}^0(\underline{q})$ on the strain and stress fields is readily seen through the examination of the strain Green function $\hat{\mathbf{R}}^0(\underline{q})$ which, when $\hat{\mathbf{G}}^0(\underline{q})$ is defined, controls the strain-based Lippmann-Schwinger equation equivalent to (D.2), i.e.

$$ \forall \underline{q} \ne (000)^*, \quad \hat{\underline{u}}(\underline{q}) = -\hat{\mathbf{R}}^0(\underline{q}) : \hat{\underline{\tau}}(\underline{q}), \qquad (D.7) $$

with

$$ \hat{\mathbf{R}}^0_{ijk}(\underline{q}) = \frac{1}{4} \left\{ D_i(\underline{q}) \hat{\mathbf{G}}^0_{jk}(\underline{q}) D_l(\underline{q}) + D_j(\underline{q}) \hat{\mathbf{G}}^0_{lk}(\underline{q}) D_l(\underline{q}) + D_l(\underline{q}) \hat{\mathbf{G}}^0_{jl}(\underline{q}) D_k(\underline{q}) + D_j(\underline{q}) \hat{\mathbf{G}}^0_{il}(\underline{q}) D_k(\underline{q}) \right\}. \qquad (D.8) $$

From Eq. (D.4) and (D.6), and using Voigt notations, it is easy to see that, for any angle $\theta$ and for $q_1 \ne \pi$, the Green function component $\hat{\mathbf{R}}^0_{IJ}(q_1, \pi + \epsilon \sin\theta, \pi + \epsilon \cos\theta) \to 0$ when $\epsilon \to 0$ except when the pair of indices $IJ$ (or the pair obtained by permutation of I and J) belongs to the set $S_{100}$ defined by

$$ S_{100} = \{22, 23, 24, 33, 34, 44, 55, 56, 66\}, \qquad (D.9) $$

where the notation $S_{100}$ is a reminder that we are considering the behaviour of the Green function in the neighbourhood of the q-line $L_{100}$. However, when $IJ$ belongs to $S_{100}$ and for $q_1 \ne \pi$, component $\hat{\mathbf{R}}^0_{IJ}(q_1, \pi + \epsilon \sin\theta, \pi + \epsilon \cos\theta)$ stays finite when we approach the line. More precisely, within the finite neighbourhood of $L_{100}$ defined by $0 \le \epsilon \ll |\tan \frac{q_1 - \pi}{2}|$ and $q_1 \ne \pi$, $\hat{\mathbf{R}}^0_{IJ}(q_1, \pi + \epsilon \sin\theta, \pi + \epsilon \cos\theta)$ becomes independent of $q_1$ and of $\epsilon$, but varies continuously with the angle $\theta$ around the $L_{100}$ line, i.e.

$$ \text{if } IJ \in S_{100}, \ 0 \le \epsilon \ll |\tan \frac{q_1 - \pi}{2}| \text{ and } q_1 \ne \pi, \quad \hat{\mathbf{R}}^0_{IJ}(q_1, \pi + \epsilon \sin\theta, \pi + \epsilon \cos\theta) \approx \hat{\mathbf{R}}^0_{IJ}(\theta), \qquad (D.10) $$

where the explicit $\theta$-dependence of the functions $\hat{\mathbf{R}}^0_{IJ}(\theta)$ are irrelevant. The important consequence of this behaviour is that, when $IJ$ belongs to $S_{100}$, there is no analytic continuation for $\hat{\mathbf{R}}^0_{IJ}(\underline{q})$ along the line $L_{100}$, which explained why the non-definitiveness of $\hat{\mathbf{R}}^0_{\underline{q}}(\underline{q})$ along $L_{100}$, and the ensuing mathematical artefacts in real space, cannot be solved by any recipe. Conversely, when $IJ$ does not belong to $S_{100}$, and except on the point $(\pi\pi\pi)^*$, the Green function component $\hat{\mathbf{R}}^0_{IJ}(\underline{q})$ may be analytically continued on the line $L_{100}$ by setting $\hat{\mathbf{R}}^0_{IJ}(\underline{q}) = 0$.
A similar analysis along lines $L_{010}$ and $L_{001}$ leads to the same results as above, replacing $S_{100}$ by $S_{010}$ or $S_{001}$ as appropriate, with $S_{010}$ and $S_{001}$ defined by

$$S_{010} = \{11, 13, 15, 33, 35, 44, 46, 55, 66\}, \qquad (D.11)$$

$$S_{001} = \{11, 12, 16, 22, 26, 44, 45, 55, 66\}.$$

When the couple of indices *IJ* does not belong to any one of the sets $S_{100}$, $S_{010}$ and $S_{001}$, the Green function component $\hat{\Gamma}_{IJ}^0(\underline{q})$ may be analytically continued along the lines $L_{100}$, $L_{010}$ and $L_{001}$ using the recipe $\hat{\Gamma}_{IJ}^0(\underline{q}) = 0$ (except at the point $q = (\pi\pi\pi)^*$, which is discussed in the next paragraph). This is the case for *IJ* (or *JI*) = 14, 25 and 36. However, for any other *IJ*, the corresponding $\hat{\Gamma}_{IJ}^0(\underline{q})$ is non-analytic along at least one of the lines $L_{100}$, $L_{010}$ and $L_{001}$.

The point $q = (\pi\pi\pi)^*$ was excluded from the previous analysis because the mathematical behaviour of the strain Green function in the neighbourhood of this point differs from the one observed above. Indeed, a simple analysis shows that, when we approach $q = (\pi\pi\pi)^*$ along any ray and for any couple *IJ*, $\hat{\Gamma}_{IJ}^0(\underline{q})$ stays finite and equal to a value that depends only on the direction of the ray, i.e.

$$\forall IJ, \quad 0 \le \epsilon \ll 1, \quad \hat{\Gamma}_{IJ}^0(\pi + \epsilon \sin\theta \cos\phi, \pi + \epsilon \sin\theta \sin\phi, \pi + \epsilon \cos\theta) \approx h_{IJ}(\theta, \phi), \qquad (D.12)$$

where the angles $\theta$ and $\phi$ define the direction of the ray in a spherical coordinate system around $\underline{q} = (\pi\pi\pi)^*$. Therefore, for any *IJ*, $\hat{\Gamma}_{IJ}^0(\underline{q})$ cannot be analytically continued to $q = (\pi\pi\pi)^*$.

These non-analytical behaviours have important consequences in real space which are best highlighted by analysing the inverse Fourier transform of the Green function $\hat{\Gamma}_{\underline{q}}^0(\underline{q})$, defined by

$$\hat{\Gamma}_{\underline{q}}^0(\underline{r}) = \frac{1}{N} \sum_{\underline{q}} \hat{\Gamma}_{\underline{q}}^0(\underline{q}) \exp(i\underline{q}\cdot\underline{r}). \qquad (D.13)$$

Consider for instance the Green function component $\hat{\Gamma}_{24}^0(\underline{q})$ which, according to the previous analysis, is only singular along the line $L_{100}$. When we approach $L_{100}$ along a ray perpendicular to the line, $\hat{\Gamma}_{24}^0(\underline{q})$ becomes constant and equal to a value that depends only on the direction of the ray defined by the angle $\theta$ (see Eq. (D.10)). Because of this “one-dimensional” non-analytical behaviour, we make the hypothesis that, in (D.13), the contribution of the neighbourhood of the line $L_{100}$ to the real space Green function $\Gamma_{24}^0(\underline{r})$ will dominate and lead to a very slow decrease of $\Gamma_{24}^0(\underline{r})$ with $\|\underline{r}\|$. This hypothesis will be supported later, when we will observe that the numerical computation of $\Gamma_{24}^0(\underline{r})$ through the exact inverse Fourier transform (D.13) will be in agreement with the result of the ongoing approximate analysis. For now, we proceed by extending (D.10) to the entire q-space, which leads to the approximation

$$\Gamma_{24}^0(l_1, l_2, l_3) \approx \frac{1}{N} \exp\{i\pi(l_2 + l_3)\} \sum_{q_1, Q_2, Q_3} \hat{\Gamma}_{24}^0(\frac{Q_2}{\|\underline{Q}\|}, \frac{Q_3}{\|\underline{Q}\|}) \exp\{i(q_1 l_1 + Q_2 l_2 + Q_3 l_3)\}, \qquad (D.14)$$

where the integers $l_1, l_2$ and $l_3$ are the coordinates of the vector $\underline{r}$ and where the vector $\underline{q} = (q_1, q_2, q_3)^*$ has been decomposed into $\underline{q} = (q_1, \pi, \pi)^* + \underline{Q}$ with $\underline{Q} = (0, q_2 - \pi, q_3 - \pi)^*$. Since the vector $\underline{Q}$ is perpendicular to the line $L_{100}$ and univocally defines the angle $\theta$, the $\theta$-dependent function $\hat{\Gamma}_{24}^0(\underline{\theta})$ defined in Eq. (D.10) has been equivalently replaced by the function $\hat{\Gamma}_{24}^0(Q_2/\|\underline{Q}\|, Q_3/\|\underline{Q}\|)$, where $\|\underline{Q}\|$ is the norm of $\underline{Q}$. Now, we proceed by considering $\underline{r}$ vectors that lie within the plane (100) that contains the origin of real space, i.e. vectors of type $\underline{r} = (0, l_2, l_3)$. Since the q-space function $\hat{\Gamma}_{24}^0(Q_2/\|\underline{Q}\|, Q_3/\|\underline{Q}\|)$ does not depend on $q_1$, we get

$$\Gamma_{24}^0(0, l_2, l_3) \approx \frac{1}{N_2 N_3} \exp\{i\pi(l_2 + l_3)\} \sum_{Q_2, Q_3} \hat{\Gamma}_{24}^0(\frac{Q_2}{\|\underline{Q}\|}, \frac{Q_3}{\|\underline{Q}\|}) \exp\{i(Q_2 l_2 + Q_3 l_3)\}. \qquad (D.15)$$

We consider now the continuous limit of the remaining discrete sums and also extend to infinity the domain of integration, which leads to

$$\hat{\Gamma}_{24}^0(0, l_2, l_3) \approx \frac{1}{(2\pi)^2} \exp\{i\pi(l_2 + l_3)\} \int_{Q_2, Q_3} dQ_2 dQ_3 \hat{\Gamma}_{24}^0(\frac{Q_2}{\|\underline{Q}\|}, \frac{Q_3}{\|\underline{Q}\|}) \exp\{i(Q_2 l_2 + Q_3 l_3)\}. \qquad (D.16)$$

We expect these approximations to be valid if $\|\underline{r}\|$ is sufficiently smaller than the system size $L$, in order to minimize the impact of the periodic boundary conditions which limit the continuous setting of the q-space, and sufficiently larger than $d$, in order to minimize the effects of the finite grid spacing which limits the extension of the q-space domain. We immediately see from Eq. (D.16) that we have the following property:

$$\hat{\Gamma}_{24}^0(0, \lambda l_2, \lambda l_3) \approx \frac{1}{\lambda^2} \exp\{i\pi\lambda(l_2 + l_3)\} \hat{\Gamma}_{24}^0(0, l_2, l_3), \qquad (D.17)$$

where, $l_2$ and $l_3$ being integers, the exponential term is equal to $\pm 1$. We conclude that, in the plane (100) which contains the origin of real space, the Green function $\hat{\Gamma}_{24}^0(\underline{r})$ oscillates and decreases as $1/\|\underline{r}\|^2$.

This slow, oscillating decay will be observed for any $\hat{\Gamma}_{IJ}^0(\underline{r})$ component if *IJ* belongs to one (and only one) of the sets $S_{100}$, $S_{010}$ or $S_{001}$ (for $\underline{r}$ vectors in the (100), (010) or (001) plane which contains the origin of real space, as appropriate).

If *IJ* does not belong to any one of the sets $S_{100}$, $S_{010}$ and $S_{001}$, we have seen above that $\hat{\Gamma}_{IJ}^0(\underline{r})$ may be analytically continued along the lines $L_{100}$, $L_{010}$ and $L_{001}$, except at the point $\underline{q} = (\pi\pi\pi)^*$. A simple analysis in the same spirit as above, with Eq. (D.10) replaced by Eq. (D.12), shows that this remaining singularity generates in real space an oscillatory decrease in $1/\|\underline{r}\|^3$.

These oscillatory decays are very different from the behaviour we would observe if the q-space Green function $\hat{\Gamma}_{\underline{q}}^0(\underline{q})$ was regular for any finite $\underline{q}$. In fact, the only singularity that should be observed is that associated with translational invariance, which therefore lies at the point $\underline{q} = (000)^*$ and which, when it is alone, generates a monotonic decrease in $1/\|\underline{r}\|^3$ for all the Green function component $\hat{\Gamma}_{IJ}^0(\underline{r})$.
<figure><img src="image_19.png" /><figcaption>Fig. D.18. Isosurfaces of the modulus of the q-space Green functions $\hat{\Gamma}_{24}^0(\underline{q})$ obtained with the rotated scheme (left panel) and with the tetrahedral stencil (right panel) (the isosurfaces have been defined by $\hat{\Gamma}_{24}^0(\underline{q}) = 0.27 \pm 0.2$, in units of $\mu^{-1}$).</figcaption></figure>

We now present numerical results that will confirm the approximate mathematical analysis that has been outlined above. Green functions obtained with the rotated scheme will be compared to those obtained with the tetrahedral stencil. For the latter, we will consider the real part of half the sum of the Green functions $\hat{\Gamma}_{\approx}^{0,d}(\underline{q})$ and $\hat{\Gamma}_{\approx}^{0,nd}(\underline{q})$ defined in Eq. (2.57). Referring to the mathematical nature of strain Green functions, which is to generate a strain as a response to a given polarization field, and taking that the final strain for the tetrahedral stencil is the arithmetic average defined in Eq. (2.49b), the collapse of the two Lippmann-Schwinger equations given in Eq. (2.60), in which we identify the two polarization fields with a single given field, indeed leads to the real part of half the sum $(\hat{\Gamma}_{\approx}^{0,d}(\underline{q}) + \hat{\Gamma}_{\approx}^{0,nd}(\underline{q}))$, quantity which, for ease of notation, we will also unambiguously designate by $\hat{\Gamma}_{\approx}^0(\underline{q})$. The elastic tensor $\lambda_{\approx}^0$ is isotropic, with shear modulus and Poisson coefficient given by $\mu = 132.30$ GPa and $\nu = 0.26$, and the periodic medium is discretized with $L = 128$ cells along each cartesian coordinates.

First, we show in Fig. D.18 isosurfaces of the modulus of $\hat{\Gamma}_{24}^0(\underline{q})$ obtained with the rotated and tetrahedral stencils (left and right panels, respectively). For the rotated scheme, we clearly observe the singular behaviour along the line $L_{100}$ discussed above: sufficiently close to the line, $\hat{\Gamma}_{24}^0(q_1, q_2, q_3)$ does not depend on $q_1$ or on the distance perpendicular to the line. For the tetrahedral stencil, we only observe, as expected, the isolated singularity associated with the point $\underline{q} = (000)^*$, reproduced by periodicity at the height corners of the q-space domain defined in Eq. (2.9).

Next, we discuss the real space Green functions $\Gamma_{24}^0(\underline{r})$, obtained numerically by inverse Fourier transform, as expressed in Eq. (D.13). In Fig. D.19, the Green functions $\Gamma_{24}^0(\underline{r})$, for $\underline{r} = (0, l_2, l_3)$, have been scaled by $\|\underline{r}\|^2$ in the left panels and by $\|\underline{r}\|^3$ in the right ones. The results obtained with the rotated scheme and those obtained with the tetrahedral stencil are shown in the top and bottom rows, respectively. In line with the analysis presented above, we observe that, with the rotated scheme and if $\|\underline{r}\|$ is sufficiently smaller than the system size $L$, $\Gamma_{24}^0(\underline{r})$ exhibits scaling by $\|\underline{r}\|^2$ and not by $\|\underline{r}\|^3$. Conversely, with the tetrahedral stencil and again if $\|\underline{r}\|$ is sufficiently smaller than the system size, we clearly observe scaling by $\|\underline{r}\|^3$ and no scaling by $\|\underline{r}\|^2$.

Note that the $\|\underline{r}\|^2$-scaling within the rotated scheme observed here numerically justifies a *posteriori* hypothesis we made above to write Eq. (D.14). Indeed, it confirms the dominant role of the $L_{100}$ singularity in the behaviour of $\Gamma_{24}^0(\underline{r})$, since in the absence of this one-dimensional singularity, the slow $1/\|\underline{r}\|^2$ decay would be replaced by a faster $1/\|\underline{r}\|^3$ oscillatory decay associated to the singularity at $\underline{q} = (\pi\pi\pi)^*$ supplemented by the $1/\|\underline{r}\|^3$ monotonic decay associated to the $\underline{q} = (000)^*$ singularity.

In brief, within the rotated scheme, a Green function component $\Gamma_{IJ}^0(\underline{r})$ that is singular along one of the lines $L_{100}$, $L_{010}$ or $L_{001}$ generates, in real space, oscillations that decrease as $1/\|\underline{r}\|^2$, which, in a 3D space, represents a very slow decrease. This is certainly at the source of the long-ranged stress oscillations observed with the rotated scheme. As shown here, the mathematical origin of these artefacts is the behaviour of the Green function in the neighbourhood of the lines $L_{100}$, $L_{010}$ and $L_{001}$, which excludes any analytic continuation of the Green function along those lines. In numerical applications, finite boxes are used, which discretizes the Fourier space on a set of discrete q-points, as in (2.9). If the box dimensions are odd, none of these q-points sit on the lines $L_{100}$, $L_{010}$ and $L_{001}$. However, the underlying singular behaviour in the continuous q-space is still there and manifests itself through the discrete q-points that are close enough to the lines. Therefore, there is no recipe to correct for the real space artefacts, whether the linear dimensions of the simulation box are odd or even.

In summary, within the rotated scheme, the strain Green function components $\Gamma_{IJ}^0(\underline{r})$ exhibit unphysical oscillations. For the few that do not belong to any of the sets $S_{100}$, $S_{010}$ and $S_{001}$, these oscillations decrease as $1/\|\underline{r}\|^3$. However, for the many that belong to one and only one of the sets $S_{100}$, $S_{010}$ or $S_{001}$, these oscillatory artefacts decay as $1/\|\underline{r}\|^2$, which is a very slow decay.

## Data availability

Raw data generated by the numerical experiments presented in the article can be provided by the author on request.
<figure><img src="image_20.png" /><figcaption>Fig. D.19. Green functions $\Gamma_{24}^0(r)$ for $r = (0, l_2, l_3)$ within the rotated scheme (top row) and the tetrahedral stencil (bottom row): scaling by $(\|r\|/L)^2$ (left panels) and by $(\|r\|/L)^3$ (right panels) (units: $\mu^{-1}$); integers $l_2$ and $l_3$ run horizontally and vertically, respectively, from $-L/2$ to $L/2-1$, with $L=128$.</figcaption></figure>

References

[1] H. Moulinec, P. Suquet, A fast numerical method for computing the linear and nonlinear mechanical properties of composites, C. R. Acad. Sci. Sér. II 318 (11) (1994) 1417-1423.

[2] H. Moulinec, P. Suquet, A numerical method for computing the overall response of nonlinear composites with complex microstructure, Comput. Methods Appl. Mech. Engrg. 157 (1-2) (1998) 69-94.

[3] M. Schneider, Review of nonlinear FFT-based computational homogenization methods, Acta Mech. 232 (6) (2021) 2051-2100.

[4] S. Lucarini, M.V. Upadhyay, J. Segurado, FFT based approaches in micromechanics: fundamentals, methods and applications, Modelling Simul. Mater. Sci. Eng. 30 (2) (2022) 023002.

[5] R.A. Lebensohn, A.D. Rollett, Spectral methods for full-field micromechanical modelling of polycrystalline materials, Comput. Mater. Sci. 173 (2020) 109336.

[6] A.G. Khachatryan, Theory of Structural Transformations in Solids, John Wiley and Sons, Inc., New York, 1983.

[7] A. Khachatryan, S. Semenovskaya, T. Tsakalakos, Elastic strain energy of inhomogeneous solids, Phys. Rev. B 52 (22) (1995) 15909.

[8] L.-Q. Chen, Phase-field models for microstructure evolution, Annu. Rev. Mater. Res. 32 (1) (2002) 113-140.

[9] D.J. Eyre, G.W. Milton, A fast numerical scheme for computing the response of composites using grid refinement, Eur. Phys. J.-Appl. Phys. 6 (1) (1999) 41-47.

[10] S. Hu, L. Chen, A phase-field model for evolving microstructures with strong elastic inhomogeneity, Acta Mater. 49 (11) (2001) 1879-1890.

[11] J. Michel, H. Moulinec, P. Suquet, A computational scheme for linear and non-linear composites with arbitrary phase contrast, Internat. J. Numer. Methods Engrg. 52 (1-2) (2001) 139-160.

[12] H. Moulinec, F. Silva, Comparison of three accelerated FFT-based schemes for computing the mechanical response of composite materials, Internat. J. Numer. Methods Engrg. 97 (13) (2014) 960-985.
[13] V. Monchiet, G. Bonnet, A polarization-based FFT iterative scheme for computing the effective properties of elastic composites with arbitrary contrast, Internat. J. Numer. Methods Engrg. 89 (11) (2012) 1419-1436.

[14] M. Schneider, Lippmann-Schwinger solvers for the computational homogenization of materials with pores, Internat. J. Numer. Methods Engrg. 121 (22) (2020) 5017-5041.

[15] V. Vinogradov, G. Milton, An accelerated FFT algorithm for thermoelastic and non-linear composites, Internat. J. Numer. Methods Engrg. 76 (11) (2008) 1678-1695.

[16] J. Michel, H. Moulinec, P. Suquet, A computational method based on augmented Lagrangians and fast Fourier transforms for composites with high contrast, CMES(Comput. Model. Eng. Sci.) 1 (2) (2000) 79-88.

[17] H. Moulinec, P. Suquet, Comparison of FFT-based methods for computing the response of composites with highly contrasted mechanical properties, Physica B 338 (1-4) (2003) 58-60.

[18] H. Moulinec, P. Suquet, G.W. Milton, Convergence of iterative methods based on Neumann series for composite materials: Theory and practice, Internat. J. Numer. Methods Engrg. 114 (2018) 1103-1130.

[19] S. Brisard, L. Dormieux, FFT-based methods for the mechanics of composites: A general variational framework, Comput. Mater. Sci. 49 (3) (2010) 663-671.

[20] S. Brisard, L. Dormieux, Combining Galerkin approximation techniques with the principle of Hashin and Shtrikman to derive a new FFT-based numerical method for the homogenization of composites, Comput. Methods Appl. Mech. Engrg. 217 (2012) 197-212.

[21] D. Rodney, Y. Le Bouar, A. Finel, Phase field methods and dislocations, Acta Mater. 51 (1) (2003) 17-30.

[22] K.S. Eloh, A. Jacques, S. Berbenni, Development of a new consistent discrete green operator for FFT-based methods to solve heterogeneous problems with eigengrains, Int. J. Plast. 116 (2019) 1-23.

[23] J. Zeman, J. Vondřejc, J. Novák, I. Marek, Accelerating a FFT-based solver for numerical homogenization of periodic media by conjugate gradients, J. Comput. Phys. 229 (21) (2010) 8065-8071.

[24] H. Grimm-Strele, M. Kabel, Runtime optimization of a memory efficient CG solver for FFT-based homogenization: implementation details and scaling results for linear elasticity, Comput. Mech. 64 (2019) 1339-1345.

[25] L. Gélébart, R. Mondon-Cancel, Non-linear extension of FFT-based methods accelerated by conjugate gradients to evaluate the mechanical behavior of composite materials, Comput. Mater. Sci. 77 (2013) 430-439.

[26] M. Kabel, T. Böhlke, M. Schneider, Efficient fixed point and Newton-Krylov solvers for FFT-based homogenization of elasticity at large deformations, Comput. Mech. 54 (6) (2014) 1497-1514.

[27] M. Schneider, An FFT-based fast gradient method for elastic and inelastic unit cell homogenization problems, Comput. Methods Appl. Mech. Engrg. 315 (2017) 846-866.

[28] M. Schneider, On the Barzilai-Borwein basic scheme in FFT-based computational homogenization, Internat. J. Numer. Methods Engrg. 118 (8) (2019) 482-494.

[29] M. Schneider, A dynamical view of nonlinear conjugate gradient methods with applications to FFT-based computational micromechanics, Comput. Mech. 66 (2020) 239-257.

[30] D. Wicht, M. Schneider, T. Böhlke, On Quasi-Newton methods in fast Fourier transform-based micromechanics, Internat. J. Numer. Methods Engrg. 121 (8) (2020) 1665-1694.

[31] F.H. Harlow, J.E. Welch, Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface, Phys. Fluids 8 (12) (1965) 2182-2189.

[32] P.-A. Geslin, B. Appolaire, A. Finel, Investigation of coherency loss by prismatic punching with a nonlinear elastic model, Acta Mater. 71 (2014) 80-88.

[33] M. Schneider, F. Ospald, M. Kabel, Computational homogenization of elasticity on a staggered grid, Internat. J. Numer. Methods Engrg. 105 (9) (2016) 693-720.

[34] F. Willot, Fourier-based schemes for computing the mechanical response of composites with accurate local fields, C. R. Méc. 343 (3) (2015) 232-245.

[35] E.H. Saenger, N. Gold, S.A. Shapiro, Modeling the propagation of elastic waves using a modified finite-difference grid, Wave Motion 31 (1) (2000) 77-92.

[36] A. Ruffini, A. Finel, Phase-field model coupling cracks and dislocations at finite strain, Acta Mater. 92 (2015) 197-208.

[37] M. Schneider, D. Merkert, M. Kabel, FFT-based homogenization for microstructures discretized by linear hexahedral elements, Internat. J. Numer. Methods Engrg. 109 (10) (2017) 1461-1489.

[38] M. Leuschner, F. Fritzen, Fourier-accelerated nodal solvers (FANS) for homogenization problems, Comput. Mech. 62 (2018) 359-392.

[39] M. Ladecký, R.J. Leute, A. Falsafi, I. Pultarová, L. Pastewka, T. Junge, J. Zeman, An optimal preconditioned FFT-accelerated finite element solver for homogenization, Appl. Math. Comput. 446 (2023) 127835.

[40] M. Schneider, Voxel-based finite elements with hourglass control in fast Fourier transform-based computational homogenization, Internat. J. Numer. Methods Engrg. 123 (24) (2022) 6286-6313.

[41] F. Willot, Y.-P. Pellegrini, Fast Fourier Transform computations and build-up of plastic deformation in 2D, elastic-perfectly plastic, pixelwise disordered porous media, in: S. Forest, D. Jeulin (Eds.), Procs. Continuum Models and Discrete Systems (CMDS 11, July 30th-August 3rd 2007), Mines ParisTech Les Presses, Paris, 2008, pp. 443-449.

[42] M. Kabel, S. Fliegener, M. Schneider, Mixed boundary conditions for FFT-based homogenization at finite strains, Comput. Mech. 57 (2016) 193-210.

[43] S. Lucarini, J. Segurado, An algorithm for stress and mixed control in Galerkin-based FFT homogenization, Internat. J. Numer. Methods Engrg. 119 (2019) 797-805.

[44] K. Sab, J. Bleyer, S. Brisard, M. Dolbeau, An FFT-based adaptive polarization method for infinitely contrasted media with guaranteed convergence, Comput. Methods Appl. Mech. Engrg. 427 (2024) 117012.

[45] M. Schneider, D. Wicht, T. Böhlke, On polarization-based schemes for the FFT-based computational homogenization of inelastic materials, Comput. Mech. 64 (2019) 1073-1095.

[46] C. Bellis, P. Suquet, Geometric variational principles for computational homogenization, J. Elasticity 137 (2) (2019) 119-149.

[47] M. Frigo, S.G. Johnson, The design and implementation of FFTW3, Proc. IEEE 93 (2) (2005) 216-231, special issue on Program Generation, Optimization, and Platform Adaptation.

[48] Q. Li, P.M. Anderson, A compact solution for the stress field from a cuboidal region with a uniform transformation strain, J. Elasticity 64 (2001) 237-245.

[49] Z. Hashin, S. Shtrikman, A variational approach to the theory of the elastic behaviour of multiphase materials, J. Mech. Phys. Solids 11 (1963) 127-140.

[50] R. Hill, Elastic properties of reinforced solids: Some theoretical principles, J. Mech. Phys. Solids 11 (1963) 357-372.

[51] M. Zecevic, R.A. Lebensohn, L. Capolungo, New large-strain FFT-based formulation and its application to model strain localization in nano-metallic laminates and other strongly anisotropic crystalline materials, Mech. Mater. 166 (2022) 104208.

[52] C. Bellis, R. Ferrier, Numerical homogenization by an adaptive Fourier spectral method on non-uniform grids using optimal transport, Comput. Methods Appl. Mech. Engrg. 419 (2024) 116658.