# Reconstructing displacements from the solution to the periodic Lippmann-Schwinger equation discretized on a uniform grid

Sébastien Brisard

► To cite this version:

Sébastien Brisard. Reconstructing displacements from the solution to the periodic Lippmann-Schwinger equation discretized on a uniform grid. International Journal for Numerical Methods in Engineering, 2017, 109 (4), pp.459-486. 10.1002/nme.5263. hal-01304603

HAL Id: hal-01304603

https://enpc.hal.science/hal-01304603

Submitted on 20 Apr 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Reconstructing displacements from the solution to the periodic Lippmann–Schwinger equation discretized on a uniform grid

S. Brisard¹*

¹Laboratoire Navier, UMR 8205, École des Ponts, IFSTTAR, CNRS, UPE, Champs-sur-Marne, France

## SUMMARY

Uniform grid solvers of the periodic Lippmann–Schwinger equation have been introduced by Moulinec and Suquet for the numerical homogenization of heterogeneous materials. Based on the fast Fourier transform, these methods use the strain as main unknown and usually do not produce displacement fields. While this is generally not perceived as a restriction for homogenization purposes, some tasks might require kinematically admissible displacement fields.

In this paper, we show how the numerical solution to the periodic Lippmann–Schwinger equation can be post-processed to reconstruct a displacement field. Our procedure is general: it applies to any variant of the Moulinec–Suquet solver. The reconstruction is formulated as an auxiliary elastic equilibrium problem of a homogeneous material, which is solved with displacement-based finite elements. Taking advantage of the periodicity, the uniformity of the grid and the homogeneity of the material, the resulting linear system is formulated and solved efficiently in Fourier space. The cost of our procedure is lower than that of one iteration of the Lippmann–Schwinger solver.

An application of this post-processing procedure is proposed, in which the reconstructed displacement field is used to compute a rigorous upper bound on the effective shear modulus of some model microstructure. Copyright © 0000 John Wiley & Sons, Ltd.

Received ...

**KEY WORDS:** Solids; Elasticity; Micromechanics; Integral equations; Variational methods; Finite element methods

This is the pre-peer reviewed version of the following article: *Reconstructing displacements from the solution to the periodic Lippmann–Schwinger equation discretized on a uniform grid*, which has been published in final form at 10.1002/nme.5263. This article may be used for non-commercial purposes in accordance with *Wiley Terms and Conditions for Self-Archiving*. Sections 1 (Introduction) and 5 (Numerical results) have been significantly reworked in the published (final) version of this article. The remainder of this article (including Sections 3 and 4 where the method itself is presented) is essentially unchanged.

## 1. INTRODUCTION

The motivation for the present paper comes from the numerical homogenization of periodic media. The reader is referred to standard textbooks for a detailed account of this theory, from the mechanical point of view [1, 2, 3] and the mathematical point of view [4, 5, 6]. In what follows, only the essential steps of the homogenization process are highlighted.

<footnote>*Correspondence to: S. Brisard, Laboratoire Navier, École des Ponts ParisTech, 6-8 avenue Blaise Pascal, Champs-sur-Marne, F-77455 Marne-la-Valle. Email: sebastien.brisard@ifsttar.fr</footnote>
Let $\Omega = (0, L_1) \times \cdots \times (0, L_d)$ be the unit cell of a heterogeneous, $\Omega$-periodic material, where $L_i$ denotes the length of the unit cell $\Omega$ in the $i$-th direction for ($i = 1, \dots, d$; $d$: dimension of the physical space). Within the framework of linear elasticity, homogenization of periodic media requires the solution to the so-called local (or corrector) problem. This problem amounts to finding the $\Omega$-periodic displacement $\mathbf{u}$ such that

$$ \nabla \cdot (\mathbf{C} : (\mathbf{E} + \nabla^s \mathbf{u})) = \mathbf{0}, \qquad (1) $$

in the whole space $\mathbb{R}^d$. In the above equations, $\mathbf{C}$ denotes the ($\Omega$-periodic) local stiffness and $\mathbf{E}$ denotes the macroscopic (imposed) strain; finally, $\nabla^s \mathbf{u}$ denotes the symmetric part of the gradient of the displacement $\mathbf{u}$. It should be noted that the solution to Eq. (1) is unique up to an additive constant, which need not be specified for the time being. In Sec. 4.4, we will explicitely impose the average displacement to be null.

From the solution to problem (1), the effective stiffness $\mathbf{C}^{\text{eff}}$ is derived as the linear operator which relates the macroscopic stress to the macroscopic strain $\mathbf{E}$

$$ \frac{1}{L_1 \cdots L_d} \int_{\Omega} \mathbf{C}(\mathbf{x}) : (\mathbf{E} + \nabla^s \mathbf{u}(\mathbf{x})) \, dx_1 \cdots dx_d = \mathbf{C}^{\text{eff}} : \mathbf{E}. \qquad (2) $$

Solving problem (1) can be a daunting task if the microstructure is highly heterogeneous (that is, if $\mathbf{C}$ is highly oscillatory within the unit cell $\Omega$). In a seminal paper, Moulinec and Suquet [7] (see also [8]) introduced a new, efficient numerical scheme on a cartesian grid. It relies on the discretization of the so-called Lippmann-Schwinger equation [9, 10, 11] which we now proceed to introduce. It is first observed that problem (1) is equivalent to finding the $\Omega$-periodic displacement $\mathbf{u}$ such that

$$ \nabla \cdot (\mathbf{C}_0 : \nabla^s \mathbf{u} + \tau) = \mathbf{0}, \qquad (3a) $$

$$ \varepsilon = \mathbf{E} + \nabla^s \mathbf{u}, \qquad (3b) $$

$$ \tau = (\mathbf{C} - \mathbf{C}_0) : \varepsilon, \qquad (3c) $$

where $\mathbf{C}_0$ denotes the elastic moduli of some arbitrary, homogeneous, reference material. In the above equations, $\varepsilon$ denotes the total strain, while $\tau$ is the so-called stress polarization. By definition of the Green operator for strains $\Gamma_0$ (see Appendix A), it is readily found that [compare Eqs. (3a) and (62)]

$$ \nabla^s \mathbf{u} = -\Gamma_0 * \tau, \qquad (4) $$

where '*' denotes the standard (periodic) convolution product. Eliminating the stress polarization $\tau$ between Eqs. (3b), (3c) and (4) leads to the following integral equation

$$ \varepsilon + \Gamma_0 * [(\mathbf{C} - \mathbf{C}_0) : \varepsilon] = \mathbf{E}, \qquad (5) $$

which is known as the Lippmann-Schwinger equation. It should be noted that in Eq. (5), the main unknown is no longer the periodic part $\mathbf{u}$ of the displacement, but the total strain $\varepsilon$ defined by Eq. (3b). Solving Eq. (5) numerically proceeds in two steps. First, the continuous equation must be discretized. Then, a suitable strategy must be adopted to find the solution to the resulting linear system.

Regarding the discretization, it is convenient to introduce the concept of discrete Green operator [12], which provides a unified framework to a wide range of discretization strategies: truncation of high frequencies [8], exact evaluation for cell-wise constant stress polarizations [13], low-pass filtering [12], finite elements [14] or finite differences [15, 16]. Although a comparison of these approaches is out of the scope of the present paper, the recent work by Willot [15, 16] deserves a special mention as the most promising strategy.

As for the solution to the resulting discretized Lippmann-Schwinger equation, various iterative schemes have been proposed, namely: the original fixed-point iterations of Moulinec and Suquet [7, 8], the accelerated scheme of Eyre and Milton [17], the augmented Lagrangian solver of Michel et al. [18], Krylov subspace methods [13, 12], or the polarization-based scheme of Monchiet and Bonnet [19] (see also [20] for a comparison of some of these schemes). All these schemes share a common matrix-free approach relying on the fast Fourier transform to efficiently compute matrix-vector products, as was first suggested by Moulinec and Suquet [7, 8].
The combination of a discrete Green operator and an iterative linear solver constitutes what we will call a Uniform Grid Periodic Lippmann-Schwinger solver in the remainder of this paper (in short, UGPLS solver). The exact discretization and solving strategies are not relevant to the present work, as long as it is understood that UGPLS solvers produce an estimate of the strain (or, equivalently, the stress-polarization), not the displacement, on a uniform cartesian grid. While in most use-cases, only the local stresses and strains are required, it is sometimes desirable to produce kinematically admissible displacement fields, as illustrated below with two examples.

As a first example, we consider the limit analysis of periodic, heterogeneous materials, which requires kinematically admissible displacement fields in order to produce rigorous upper bounds on the macroscopic strength criterion. This can be achieved by means of dedicated finite element techniques [21, 22]. Alternatively, the problem can be reformulated as the homogenization of a non-linear, viscous material [2], for which UGPLS solvers can be used iteratively (see [7] for the first application of UGPLS solvers to non-linear materials). The numerical solution is then post-processed, delivering a rigorous upper bound on the collapse load. It should be noted that this approach requires the local strength criterion to be smooth or regularized [23].

As a second example, we mention a posteriori evaluation of the error in constitutive relation as first introduced by Ladevèze and Leguillon [24] within the framework of finite element analysis. From the Prager-Synge theorem, it can be shown that the error in constitutive relation provides an upper bound on the approximation error in the energy norm (see [25] for a review, or the textbook by Ladevèze [26]). Evaluation of the error in constitutive law requires both a kinematically admissible displacement field (which is addressed in the present paper), and a statically admissible stress field (which will be presented in a paper to come).

To sum up the above discussion, UGPLS schemes are attractive techniques to solve the corrector problem in an efficient way, but fail to deliver a kinematically admissible displacement field. In the present paper, we seek to overcome this shortcoming and propose an a posteriori construction of a kinematically admissible displacement field. The proposed method applies to any UGPLS solver (regardless of the discretization strategy and the linear solver). The ideas underlying the construction are described now.

Let $\epsilon^N$ denote the discrete strain field resulting from a UGPLS analysis, where the supscript 'N' is a discretization parameter which relates to the fineness of the grid (see Sec. 3.1). In a Galerkin setting, $\epsilon^N$ can be considered as a cell-wise constant estimate of the true strain $\epsilon$ [12]. From Eq. (3c), the discrete stress polarization $\tau^N$ is evaluated. Likewise, it is seen as a cell-wise constant estimate of the true stress polarization $\tau$.

If the true stress polarization $\tau$ were known exactly, then the displacement $\mathbf{u}$ would be the periodic solution to Eq. (3a). We therefore seek the approximate displacement field $\mathbf{u}^N$ as the periodic solution to Eq. (3a) where the true stress polarization $\tau$ is replaced with the approximate, cell-wise constant, $\tau^N$

$$ \nabla \cdot (\mathbf{C}_0 : \nabla^s \mathbf{u}^N + \tau^N) = \mathbf{0}. \qquad (6) $$

In other words, reconstructing the displacement field resulting from a UGPLS analysis reduces to solving the problem of the linear elastic equilibrium of a prestressed, homogeneous material with periodic boundary conditions. Of course, the exact solution to this problem is not known; rather, a finite element estimate of this solution is seeked. The remainder of this paper is therefore devoted to the numerical solution of the following problem (formulated on the whole space $\mathbb{R}^d$)

$$ \nabla \cdot (\mathbf{C}_0 : \nabla^s \mathbf{u} + \boldsymbol{\omega}) + \mathbf{b} = \mathbf{0}, \qquad (7) $$

where the unknown displacement $\mathbf{u}$ is $\Omega$-periodic. In the above equation, $\boldsymbol{\omega}$ is the ($\Omega$-periodic) prestress. Problem (7) is slightly more general than problem (6) due to the body forces with $\Omega$-periodic volume density $\mathbf{b}$.

Problem (7) is solved numerically by means of a finite element technique. The displacement is discretized over a $d$-dimensional cartesian grid; interpolation is achieved through P1 shape functions. Periodicity of the discretized displacement is enforced by construction (using periodic shape functions).
Due to the periodic boundary conditions and the homogeneity of the material, the finite element approximation of Eq. (7) results in a *block-circulant* stiffness matrix. This suggests to formulate this approximation in the Fourier space, where all frequencies are uncoupled. More precisely, while the *nodal* displacements are the solution to a ($N_{nodes}d$) × ($N_{nodes}d$) system (where $N_{nodes}$ is the total number of nodes), the *modal* displacements (that is, the discrete Fourier transform of the nodal displacements) are the solution to $N_{nodes}$ systems of size $d \times d$, which can be solved at virtually no cost. Of course, this approach requires the computation of direct and inverse discrete Fourier transforms: using the fast Fourier transform guarantees the overall efficiency of the proposed method.

At this point, the alternative approach of Vondřejc and collaborators (see [27, 28] and references therein) should be mentioned. These authors propose a Galerkin discretization of the corrector problem with trigonometric polynomials: their approach therefore delivers non-local, kinematically admissible, approximations of the displacement field. Local approximations might sometimes be more convenient, and the approach proposed in the present paper should be seen as complementary to that of Vondřejc and collaborators. It should also be noted that, since high frequencies are discarded from the approximation space, the displacement and strain fields delivered by Vondřejc and collaborators might be subject to the Gibbs phenomenon.

The remainder of this paper is organized as follows. Sec. 2 introduces some convenient element-by-element operations on vectors, which allow for more compact expressions; basic results relating to Fourier series and the discrete Fourier transform are also recalled. We then define in Sec. 3 the discretization grid and the discretized displacement and strain. In particular, we relate (without approximation) the Fourier series expansion of the piecewise affine displacement field to the discrete Fourier transform of the nodal displacements. Sec. 4 constitutes the crux of the paper: we show that the global stiffness matrix and the vector of nodal forces have simple expressions in Fourier space. We therefore introduce the *modal* stiffness ($d \times d$) matrices and force ($d \times 1$) vectors. A summary of the method is provided in Sec. 4.5. This paper closes on a 2D (plane strain elasticity) example (see Sec. 5). The reconstructed displacement field is used to produce rigorous upper bounds on the effective elastic moduli.

## 2. PRELIMINARIES

### 2.1. Element-by-element operations on the components of vectors

In the present paper, it will be convenient to consider vectors of $\mathbb{R}^d$ ($\mathbb{N}^d$, $\mathbb{Z}^d$) as tuples of size $d$, and define non-linear operations on the coordinates of the vectors. We first define element-by-element product and division of two vectors $\mathbf{x}$ and $\mathbf{y}$ as follows

$$ \mathbf{x}\mathbf{y} = \sum_{j=1}^{d} x_j y_j \mathbf{e}_j \qquad (8a) $$

$$ \frac{\mathbf{x}}{\mathbf{y}} = \sum_{j=1}^{d} \frac{x_j}{y_j} \mathbf{e}_j. \qquad (8b) $$

Then, $f$ denoting a scalar function of a scalar variable, we define $f(\mathbf{x})$ for any vector $\mathbf{x} \in \mathbb{R}^d$ as the product

$$ f(\mathbf{x}) = \prod_{j=1}^{d} f(x_j). \qquad (9) $$

The above definitions lead to very compact expressions, which can be read naturally, as though they were written in a one-dimensional space. This allows the reader to concentrate on the overall structure of the equations, rather than their (dimension dependent) details. As a first example,

$$ \mathrm{sinc}\left(\frac{\pi}{N}(\mathbf{k} + \mathbf{p}N)\right) = \prod_{i=1}^{d} \mathrm{sinc}\left(\frac{\pi}{N_j}(k_j + p_j N_j)\right), \qquad (10) $$
where the sinc function (cardinal sine) is defined as follows

$$ \text{sinc } x = \begin{cases} 1 & \text{if } x = 0, \\ x^{-1} \sin x & \text{otherwise.} \end{cases} \qquad (11) $$

As further examples, the following identities resulting from rules (8) and (9) will be extensively used below

$$ |\mathbf{x}| = |x_1| \cdots |x_d| = |x_1 \cdots x_d| \qquad (12) $$

$$ \exp\left(2i\pi \frac{\mathbf{k}\mathbf{x}}{\mathbf{L}}\right) = \exp\left(2i\pi \sum_{j=1}^{d} \frac{k_j x_j}{L_j} e_j\right) = \prod_{j=1}^{d} \exp\left(2i\pi \frac{k_j x_j}{L_j}\right) = \exp\left(2i\pi \sum_{j=1}^{d} \frac{k_j x_j}{L_j}\right), \qquad (13) $$

in particular, $|\mathbf{N}| = N_1 \cdots N_d$ is the total number of elements in the cartesian grid described in Sec. 3.1.

It should be noted that rules (8) and (9) have a shortcoming: they are not intrinsic operations. In other words, the results of these operations depend on the basis in which they are performed. However, it is observed that in the present setting, only the basis which is attached to the sides of the rectangular box $\Omega$ is relevant, and no change of basis will be performed or required. Therefore, we will continue to use bold letters to denote vectors and tensors, even if they should really be understood as their arrays of components in the global basis.

## 2.2. Fourier series and the discrete Fourier transform

In the present section, we briefly recall the essential formulas regarding Fourier series and the discrete Fourier transform. In particular, we clarify the normalization constants that we adopt.

**Fourier series** Under some regularity assumptions which are not recalled here, any $\Omega$-periodic function $f$ can be expressed as the sum of a Fourier series

$$ f = \sum_{\mathbf{k} \in \mathbb{Z}^d} \tilde{f}_{\mathbf{k}} \Psi_{\mathbf{k}}, \qquad (14) $$

where the basis functions $\Psi_{\mathbf{k}}$ are defined as follows for $\mathbf{k} \in \mathbb{Z}^d$

$$ \Psi_{\mathbf{k}}(\mathbf{x}) = \exp\left(2i\pi \frac{\mathbf{k}\mathbf{x}}{\mathbf{L}}\right), \qquad (15) $$

and the Fourier coefficients $\tilde{f}_{\mathbf{k}}$ are related to $f$

$$ \tilde{f}_{\mathbf{k}} = \frac{1}{|\mathbf{L}|} \int_{\Omega} f(\mathbf{x}) \exp\left(-2i\pi \frac{\mathbf{k}\mathbf{x}}{\mathbf{L}}\right) d\mathbf{x}_1 \cdots d\mathbf{x}_d. \qquad (16) $$

Parseval's theorem then relates the $L^2$ Hermitian product of two $\Omega$-periodic functions $f$ and $g$ to the $\ell^2$ Hermitian product of their Fourier coefficients

$$ \frac{1}{|\mathbf{L}|} \int_{\Omega} \bar{f}(\mathbf{x}) g(\mathbf{x}) d\mathbf{x}_1 \cdots d\mathbf{x}_d = \sum_{\mathbf{k} \in \mathbb{Z}^d} \bar{f}_{\mathbf{k}} \tilde{g}_{\mathbf{k}}, \qquad (17) $$

where overlined quantities denote complex conjugates.

**The discrete Fourier transform** We now consider a sequence of numbers $(y_n)_{n \in \mathbb{Z}^d}$, which we assume $\mathbf{N}$-periodic ($\mathbf{N} \in \mathbb{N}^d$), that is $y_{\mathbf{n}+\mathbf{p}\mathbf{N}} = y_{\mathbf{n}}$ for all $\mathbf{n}, \mathbf{p} \in \mathbb{Z}^d$. Clearly, this infinite sequence is fully defined by the finite sequence $(y_n)_{n \in \mathcal{N}}$, where $\mathcal{N}$ denotes the following set of multi-indices

$$ \mathcal{N} = \{\mathbf{n} \in \mathbb{N}^d, 0 \le n_i < N_i\}. \qquad (18) $$
The discrete Fourier transform of $(y_n)_{n \in \mathbb{Z}^d}$ is the $\mathbf{N}$-periodic sequence $(\hat{y}_k)_{k \in \mathbb{Z}^d}$

$$ \hat{y}_k = \sum_{n \in \mathcal{N}} \exp\left(-2i\pi \frac{kn}{N}\right) y_n, \qquad (19) $$

and we have the inversion formula

$$ y_n = \frac{1}{|\mathbf{N}|} \sum_{\mathbf{k} \in \mathcal{N}} \exp\left(2i\pi \frac{\mathbf{kn}}{N}\right) \hat{y}_k. \qquad (20) $$

The Plancherel theorem relates the $l^2$ Hermitian product of two $\mathbf{N}$-periodic sequences $(y_n)_{n \in \mathbb{Z}^d}$ and $(z_n)_{n \in \mathbb{Z}^d}$ to the $l^2$ Hermitian product of their discrete Fourier transforms

$$ \sum_{n \in \mathcal{N}} \overline{y_n} z_n = \frac{1}{|\mathbf{N}|} \sum_{\mathbf{k} \in \mathcal{N}} \overline{\hat{y}_k} \hat{z}_k. \qquad (21) $$

### 3. P1 DISCRETIZATION ON A PERIODIC, UNIFORM GRID

In the present section, we define the discretized displacement and strain fields. Both are $\Omega$-periodic; as such, they can be expanded into Fourier series. Since the displacement is piecewise affine, we show that the Fourier series expansions take very simple forms. It will be seen that the discrete Fourier transform of the nodal values of the displacement arises naturally in the derivation, without any approximation. The fast Fourier transform can therefore be used to switch efficiently to and from Fourier space. In Sec. 3.1, we define the discretization grid, and show that cell-averages can be computed in Fourier space. Then, we define in Sec. 3.2 the discretized displacement, and show that its Fourier series expansion is completely defined by the discrete Fourier transform of the nodal displacements. Finally, spatial derivation of the discretized displacement leads in Sec. 3.3 to the discretized strain.

#### 3.1. The discretization grid

In the present section, we define the cartesian grid which is to be used to discretize the mechanical fields resulting from problem (7). In the i-th direction ($i = 1, \dots, d$), the side of the unit cell $\Omega$ is divided into $N_i$ equal segments of size $L_i/N_i$ (see Fig. 1). As a result, the vertices of the grid are located at

$$ \mathbf{x}_n^N = \frac{\mathbf{nL}}{N} \quad (\mathbf{n} \in \mathbb{Z}^d), \qquad (22) $$

where $\mathbf{L} = L_i e_i$ and $\mathbf{N} = N_i e_i$. In the above formula, the element-by-element operations introduced in Sec. 2.1 have been used. We further define the grid-cells $\Omega_n^N$ as follows

$$ \Omega_n^N = \{ \mathbf{x} \in \Omega, n_i L_i \le N_i x_i < (n_i + 1) L_i, i = 1, \dots, d \}, \qquad (23) $$

so that

$$ \Omega = \bigcup_{n \in \mathcal{N}} \Omega_n^N, \qquad (24) $$

where the set $\mathcal{N}$ previously defined by Eq. (18) should be understood as the set of node multi-indices.

We will need in what follows expressions of the cell-averages of a $\Omega$-periodic function $f$, as well as their discrete Fourier transform. We introduce the following notations for these quantities

$$ \mathrm{mean}_n^N(f) = \frac{|\mathbf{N}|}{|\mathbf{L}|} \int_{\Omega_n^N} f(\mathbf{x}) \, dx_1 \cdots dx_d, \qquad (25a) $$

$$ \widehat{\mathrm{mean}}_k^N(f) = \sum_{n \in \mathcal{N}} \mathrm{mean}_n^N(f) \exp\left(-2i\pi \frac{kn}{N}\right). \qquad (25b) $$
<figure><img src="image_2.png" /><figcaption>Figure 1. The discretization grid in the 3D case. In the present example, $N_1 = 4$, $N_2 = 5$ and $N_3 = 3$.</figcaption></figure>

It should be noted that the cell-averages of a Ω-periodic function are N-periodic; therefore, it is meaningful to compute the discrete Fourier transform of the cell-averages of a Ω-periodic function (see Sec. 2.2). Definitions (25a) and (25b) are applied below to the Fourier basis functions $\Psi_k$ defined in Eq. (15). We first find by straightforward integration

$$ \text{mean}_{n}^{N}(\Psi_{k}) = \text{sinc}\left(\frac{\pi k}{N}\right) \exp\left(i\frac{\pi k}{N}(2n + 1)\right), \qquad (26) $$

where **1** denotes the vector $\mathbf{1} = \mathbf{e}_1 + \cdots + \mathbf{e}_d$, and the sinc function has been defined above [see Eq. (11)]. Taking the discrete Fourier transform of both sides of Eq. (26), we find, for $\mathbf{h} \in \mathcal{N}$

$$ \widehat{\text{mean}}_{\mathbf{h}}^{N}(\Psi_{\mathbf{k}}) = \begin{cases} |\mathbf{N}| \text{sinc}\left(\frac{\pi \mathbf{k}}{\mathbf{N}}\right) \exp\left(i\frac{\pi \mathbf{k}}{\mathbf{N}}\right) & \text{if } \mathbf{k} = \mathbf{h} + \mathbf{p}\mathbf{N} \quad (\mathbf{p} \in \mathbb{Z}^d), \\ 0 & \text{otherwise.} \end{cases} \qquad (27) $$

The above formulas can be used to compute the cell-averages of any Ω-periodic function $f$ with Fourier coefficients $\tilde{f}_k$. In particular, the following expression is found for the discrete Fourier transform of the cell-averages of $f$

$$ \widehat{\text{mean}}_{\mathbf{k}}^{N}(f) = |\mathbf{N}| \sum_{\mathbf{p} \in \mathbb{Z}^d} \text{sinc}\left(\frac{\pi}{\mathbf{N}}(\mathbf{k} + \mathbf{p}\mathbf{N})\right) \exp\left(i\frac{\pi}{\mathbf{N}}(\mathbf{k} + \mathbf{p}\mathbf{N})\right) \tilde{f}_{\mathbf{k}+\mathbf{p}\mathbf{N}}, \qquad (28a) $$

with

$$ f = \sum_{\mathbf{k} \in \mathbb{Z}^d} \tilde{f}_{\mathbf{k}} \Psi_{\mathbf{k}}. \qquad (28b) $$

## 3.2. The discretized displacement and its Fourier series expansion

Let $\mathbf{u}^N$ denote the P1 finite element approximation of the solution to problem (7), discretized over the cartesian grid $(\Omega_n^N)_{n \in \mathcal{N}}$ defined in Sec. 3.1. It is expressed as follows

$$ \mathbf{u}^{N}(\mathbf{x}) = \sum_{\mathbf{n} \in \mathcal{N}} \Phi^{N}(\mathbf{x} - \mathbf{x}_{\mathbf{n}}^{N}) \mathbf{u}_{\mathbf{n}}^{N}, \qquad (29) $$
<figure><img src="image_3.png" /><figcaption>Figure 2. The periodic shape function $\Phi^N$ in the 2D case, with $N_1 = 3$, $N_2 = 4$.</figcaption></figure>

where $\Phi^N$ denotes the $\Omega$-periodic P1 shape function (to be defined below), and $\mathbf{u}_n^N$ are the nodal values of the discretized displacement

$$ \mathbf{u}_n^N = \mathbf{u}^N(\mathbf{x}_n^N) \quad (\mathbf{n} \in \mathbb{Z}^d), \qquad (30) $$

which are $\mathbf{N}$-periodic owing to the $\Omega$-periodicity of the discretized displacement. As already argued in Sec. 2.2, the discretized displacement is therefore fully defined by its values at nodes $\mathbf{n} \in \mathcal{N}$.

The $\Omega$-periodic shape function $\Phi^N$ is first defined for $\mathbf{x} \in \mathbb{R}^d$ such that $|x_j| \le L_j/2$ for all $j = 1, \dots, d$

$$ \Phi^N(\mathbf{x}) = \begin{cases} \prod_{j=1}^{d} \left(1 - \frac{N_j |x_j|}{L_j}\right) & \text{if } |x_j| \le \frac{L_j}{N_j} \text{ for all } j = 1, \dots, d, \\ 0 & \text{otherwise,} \end{cases} \qquad (31) $$

it is then extended to $\mathbb{R}^d$ by $\Omega$-periodicity (see Fig. 2). Thus defined, $\Phi^N$ is the periodic counterpart of the usual P1 shape function. Its Fourier series expansion reads (see Appendix B)

$$ \Phi^N(\mathbf{x}) = \frac{1}{|\mathbf{N}|} \sum_{\mathbf{k} \in \mathbb{Z}^d} \text{sinc}^2 \left(\frac{\pi \mathbf{k}}{\mathbf{N}}\right) \Psi_k(\mathbf{x}), \qquad (32) $$

and, upon substitution in Eq. (29)

$$ \mathbf{u}^N(\mathbf{x}) = \frac{1}{|\mathbf{N}|} \sum_{\mathbf{k} \in \mathbb{Z}^d} \text{sinc}^2\left(\frac{\pi \mathbf{k}}{\mathbf{N}}\right) \Psi_k(\mathbf{x}) \sum_{\mathbf{n} \in \mathcal{N}} \exp\left(-2i\pi \frac{\mathbf{k} \cdot \mathbf{n}}{\mathbf{N}}\right) \mathbf{u}_n^N, \qquad (33) $$

where Eq. (22) has been used.

In the above expression, the inner sum is recognized as the discrete Fourier transform of the nodal displacements $(\mathbf{u}_n^N)_{n \in \mathcal{N}}$, so that we finally find the compact expression of the Fourier series expansion of the discretized displacement $\mathbf{u}^N$

$$ \mathbf{u}^N(\mathbf{x}) = \sum_{\mathbf{k} \in \mathbb{Z}^d} \Psi_k(\mathbf{x}) \tilde{\mathbf{u}}_k^N, \qquad (34) $$

where the Fourier coefficients $\tilde{\mathbf{u}}_k^N$ are given by the following formula

$$ \tilde{\mathbf{u}}_k^N = \frac{1}{|\mathbf{N}|} \text{sinc}^2\left(\frac{\pi \mathbf{k}}{\mathbf{N}}\right) \hat{\mathbf{u}}_k^N. \qquad (35) $$
Evaluation of the potential energy of the external loads (see Sec. 4.3) further requires the DFT of the cell-averages of the displacement, which can be expressed as follows (see Appendix C.1)

$$ \widehat{\text{mean}}_{\mathbf{k}}^{\mathbf{N}}(\mathbf{u}^{\mathbf{N}}) = \cos\left(\frac{\pi\mathbf{k}}{N}\right) \exp\left(i\frac{\pi\mathbf{k}}{N}\right) \hat{\mathbf{u}}_{\mathbf{k}}^{\mathbf{N}}. \quad (36) $$

### 3.3. The discretized strain and its Fourier series expansion

Let $\boldsymbol{\varepsilon}^{\mathbf{N}} = \nabla^{\mathbf{s}}\mathbf{u}^{\mathbf{N}}$ be the strain associated to the discretized displacement $\mathbf{u}^{\mathbf{N}}$. Its Fourier series expansion is readily obtained from the differentiation of Eq. (34)

$$ \boldsymbol{\varepsilon}(\mathbf{x}) = \sum_{\mathbf{k} \in \mathbb{Z}^d} \Psi_{\mathbf{k}}(\mathbf{x}) \tilde{\boldsymbol{\varepsilon}}_{\mathbf{k}}^{\mathbf{N}}, \quad (37) $$

where the following expression is found for the Fourier coefficients $\tilde{\boldsymbol{\varepsilon}}_{\mathbf{k}}^{\mathbf{N}}$

$$ \tilde{\boldsymbol{\varepsilon}}_{\mathbf{k}}^{\mathbf{N}} = 2i\pi \tilde{\mathbf{u}}_{\mathbf{k}}^{\mathbf{N}} \otimes \frac{\mathbf{k}}{L} = \frac{2i\pi}{|N|} \mathrm{sinc}^2\left(\frac{\pi\mathbf{k}}{N}\right) \hat{\mathbf{u}}_{\mathbf{k}}^{\mathbf{N}} \otimes \frac{\mathbf{k}}{L}. \quad (38) $$

Combining Eqs. (28), (37) and (38) then leads to the following compact expression of the discrete Fourier transform of the cell averages of the strain

$$ \widehat{\text{mean}}_{\mathbf{k}}^{\mathbf{N}}(\boldsymbol{\varepsilon}^{\mathbf{N}}) = \hat{\mathbf{u}}_{\mathbf{n}}^{\mathbf{N}} \otimes \hat{\mathbf{B}}_{\mathbf{k}}^{\mathbf{N}}, \quad (39) $$

where the modal strain-displacement vector $\hat{\mathbf{B}}_{\mathbf{k}}^{\mathbf{N}}$ has been introduced

$$ \hat{\mathbf{B}}_{\mathbf{k}}^{\mathbf{N}} = 2i\pi \sum_{\mathbf{p} \in \mathbb{Z}^d} \mathrm{sinc}^3\left(\frac{\pi}{N}(\mathbf{k} + \mathbf{p}N)\right) \exp\left(i\frac{\pi}{N}(\mathbf{k} + \mathbf{p}N)\right) \frac{\mathbf{k} + \mathbf{p}N}{L}. \quad (40) $$

A lengthy but straightforward (see Appendix C.2) calculation shows that

$$ \hat{\mathbf{B}}_{\mathbf{k}}^{\mathbf{N}} \cdot \mathbf{e}_1 = \frac{2iN_1}{L_1} \sin \frac{\pi k_1}{N_1} \cos \frac{\pi k_2}{N_2} \cdots \cos \frac{\pi k_d}{N_d} \exp i\pi \left( \frac{k_1}{N_1} + \cdots + \frac{k_d}{N_d} \right), \quad (41) $$

and the other components of the modal strain-displacement vector are deduced by circular permutations of the indices.

## 4. FINITE ELEMENT APPROXIMATION OF THE PROBLEM

The approximate nodal displacements for the solution to problem (7) are found by minimization of the total potential energy $\Pi = U - V$. The strain energy $U = U_e + U_p$ is the sum of an elastic contribution $U_e$ and a term arising from the prestress $U_p$, and we have

$$ U_e = \int_{\Omega} \frac{1}{2} \boldsymbol{\varepsilon}(\mathbf{x}) : \mathbf{C} : \boldsymbol{\varepsilon}(\mathbf{x}) \, dx_1 \cdots dx_d, \quad (42a) $$

$$ U_p = \int_{\Omega} \boldsymbol{\sigma}(\mathbf{x}) : \boldsymbol{\varepsilon}(\mathbf{x}) \, dx_1 \cdots dx_d, \quad (42b) $$

$$ V = \int_{\Omega} \mathbf{b}(\mathbf{x}) \cdot \mathbf{u}^{\mathbf{N}}(\mathbf{x}) \, dx_1 \cdots dx_d, \quad (42c) $$

where $V$ denotes the potential of the external loads. The three above contributions are evaluated in Secs. 4.1, 4.2 and 4.3, respectively, by substitution of the discretized displacement $\mathbf{u}^{\mathbf{N}}$ (defined in Sec. 3.2) or strain $\boldsymbol{\varepsilon}^{\mathbf{N}}$ (derived in Sec. 3.3) into Eqs. (42a), (42b) and (42c). The total potential energy is then computed and optimized in Sec. 4.4, leading to a series of low-dimensional linear systems in the Fourier space. A summary of the method is provided in Sec. 4.5.
### 4.1. Evaluation of the strain energy (elastic contribution)

Owing to the homogeneity of the material, the elastic part $U_e$ of the strain energy [see Eq. (42a)] can readily be computed in Fourier space by means of Parseval's identity (17)

$$ 
\begin{aligned}
U_e &= \int_{\Omega} \frac{1}{2} \varepsilon^{\mathrm{N}}(\mathbf{x}) : \mathbf{C}_0 : \varepsilon^{\mathrm{N}}(\mathbf{x}) \, dx_1 \cdots dx_d = \frac{|\mathbf{L}|}{2} \sum_{\mathbf{k} \in \mathbb{Z}^d} \tilde{\varepsilon}_{\mathbf{k}}^{\mathrm{N}} : \mathbf{C}_0 : \tilde{\varepsilon}_{\mathbf{k}}^{\mathrm{N}} \\
&= 2\pi^2 \frac{|\mathbf{L}|}{|\mathbf{N}|^2} \sum_{\mathbf{k} \in \mathbb{Z}^d} \mathrm{sinc}^4 \left(\frac{\pi \mathbf{k}}{\mathbf{N}}\right) \overline{\hat{\mathbf{u}}_{\mathbf{k}}^{\mathrm{N}}} \cdot \mathbf{A}_0(\mathbf{k}/\mathbf{L}) \cdot \hat{\mathbf{u}}_{\mathbf{k}}^{\mathrm{N}}
\end{aligned} 
\quad (43) $$

where the Fourier coefficients $\tilde{\varepsilon}_{\mathbf{k}}^{\mathrm{N}}$ of the strain are given by Eq. (38) and $\mathbf{A}_0$ denotes the acoustic tensor, $\mathbf{A}_0(\mathbf{q}) = \mathbf{q} \cdot \mathbf{C}_0 \cdot \mathbf{q}$. Replacing in the above series $\mathbf{k} \in \mathbb{Z}^d$ with $\mathbf{k} + \mathbf{p}\mathbf{N}$ (with $\mathbf{k} \in \mathcal{N}$ and $\mathbf{p} \in \mathbb{Z}^d$), and taking advantage of the $\mathbf{N}$-periodicity of the discrete Fourier transform $\hat{\mathbf{u}}_{\mathbf{k}}^{\mathrm{N}}$, we find

$$ U_e = 2\pi^2 \frac{|\mathbf{L}|}{|\mathbf{N}|^2} \sum_{\mathbf{k} \in \mathcal{N}} \sum_{\mathbf{p} \in \mathbb{Z}^d} \mathrm{sinc}^4 \left(\frac{\pi}{\mathbf{N}} (\mathbf{k} + \mathbf{p}\mathbf{N})\right) \overline{\hat{\mathbf{u}}_{\mathbf{k}}^{\mathrm{N}}} \cdot \mathbf{A}_0 \left(\frac{\mathbf{k} + \mathbf{p}\mathbf{N}}{\mathbf{L}}\right) \cdot \hat{\mathbf{u}}_{\mathbf{k}}^{\mathrm{N}}. \quad (44) $$

We then introduce the so-called modal stiffness matrices $\hat{\mathbf{K}}_{\mathbf{k}}^{\mathrm{N}}$ defined as follows

$$ \hat{\mathbf{K}}_{\mathbf{k}}^{\mathrm{N}} = 4\pi^2 \sum_{\mathbf{p} \in \mathbb{Z}^d} \mathrm{sinc}^4 \left(\frac{\pi}{\mathbf{N}} (\mathbf{k} + \mathbf{p}\mathbf{N})\right) \mathbf{A}_0 \left(\frac{\mathbf{k} + \mathbf{p}\mathbf{N}}{\mathbf{L}}\right), \quad (45) $$

and expression (44) of the strain energy reduces to

$$ U_e = \frac{|\mathbf{L}|}{2|\mathbf{N}|^2} \sum_{\mathbf{k} \in \mathcal{N}} \overline{\hat{\mathbf{u}}_{\mathbf{k}}^{\mathrm{N}}} \cdot \hat{\mathbf{K}}_{\mathbf{k}}^{\mathrm{N}} \cdot \hat{\mathbf{u}}_{\mathbf{k}}^{\mathrm{N}}. \quad (46) $$

As expected, it is observed that for $\mathbf{k} = \mathbf{0}$, $\hat{\mathbf{K}}_0^N = \mathbf{0}$. Indeed, the strain energy associated with constant displacements (associated with the $\mathbf{k} = \mathbf{0}$ mode) is zero. For isotropic materials with shear modulus $\mu_0$ and Poisson ratio $\nu_0$, a lengthy but straightforward calculation (see Appendix C.3) leads to the following closed-form expression of the modal stiffness matrices

$$ \hat{\mathbf{K}}_{\mathbf{k}}^{\mathrm{N}} = \mu_0 \mathrm{tr}(\hat{\mathbf{H}}_{\mathbf{k}}^{\mathrm{N}}) \mathbf{I} + \frac{\mu_0}{1 - 2\nu_0} \hat{\mathbf{H}}_{\mathbf{k}}^{\mathrm{N}}, \quad (47) $$

and the diagonal and off-diagonal coefficients of $\hat{\mathbf{H}}_{\mathbf{k}}^{\mathrm{N}}$ are fully defined (up to a circular permutation of the indices) by the following identities

$$ \mathbf{e}_1 \cdot \hat{\mathbf{H}}_{\mathbf{k}}^{\mathrm{N}} \cdot \mathbf{e}_1 = \frac{2}{3^{d-1}} \left(\frac{N_1}{L_1}\right)^2 \left(1 - \cos \frac{2\pi k_1}{N_1}\right) \prod_{j=2}^{d} \left(2 + \cos \frac{2\pi k_j}{N_j}\right), \quad (48a) $$

$$ \mathbf{e}_1 \cdot \hat{\mathbf{H}}_{\mathbf{k}}^{\mathrm{N}} \cdot \mathbf{e}_2 = \frac{1}{3^{d-2}} \frac{N_1 N_2}{L_1 L_2} \sin \frac{2\pi k_1}{N_1} \sin \frac{2\pi k_2}{N_2} \prod_{j=3}^{d} \left(2 + \cos \frac{2\pi k_j}{N_j}\right). \quad (48b) $$

### 4.2. Evaluation of the strain energy (contribution of the prestress)

We now turn to the second term $U_p$ of the strain energy [see Eq. (42b)]. We assume that the prescribed prestress $\boldsymbol{\omega}$ is cell-wise constant, and denote $\boldsymbol{\omega}_n^N$ its value in cell $\Omega_n^N$. Then, from the Plancherel theorem (21)

$$ U_p = \frac{|\mathbf{L}|}{|\mathbf{N}|} \sum_{\mathbf{n} \in \mathcal{N}} \boldsymbol{\omega}_n^N : \mathrm{mean}_n^N(\boldsymbol{\varepsilon}^N) = \frac{|\mathbf{L}|}{|\mathbf{N}|^2} \sum_{\mathbf{k} \in \mathcal{N}} \overline{\hat{\boldsymbol{\omega}}_{\mathbf{k}}^N} : \overline{\mathrm{mean}_k^N(\boldsymbol{\varepsilon}^N)}, \quad (49) $$

and, upon substitution of Eq. (39)

$$ U_p = \frac{|\mathbf{L}|}{|\mathbf{N}|^2} \sum_{\mathbf{k} \in \mathcal{N}} \overline{\hat{\mathbf{u}}_{\mathbf{k}}^{\mathrm{N}}} \cdot \hat{\boldsymbol{\omega}}_{\mathbf{k}}^{\mathrm{N}} \cdot \overline{\hat{\mathbf{B}}_{\mathbf{k}}^{\mathrm{N}}}. \quad (50) $$
### 4.3. Potential energy of the external loads

The volume density of body forces **b** is assumed cell-wise constant in Eq. (42c): $\mathbf{b}_n^N$ denotes its value in cell $\Omega_n^N$. Then, using again the Plancherel theorem (21)

$$V = \frac{|\mathbf{L}|}{|\mathbf{N}|} \sum_{n \in \mathcal{N}} \mathbf{b}_n^N \cdot \text{mean}_n^N(\mathbf{u}^N) = \frac{|\mathbf{L}|}{|\mathbf{N}|^2} \sum_{n \in \mathcal{N}} \hat{\mathbf{b}}_n^N \cdot \overline{\text{mean}_k^N(\mathbf{u}^N)}, \quad (51)$$

which, upon substitution of Eq. (36), leads to the following expression

$$V = \frac{|\mathbf{L}|}{|\mathbf{N}|^2} \sum_{n \in \mathcal{N}} \hat{\mathbf{F}}_k^N \cdot \overline{\hat{\mathbf{u}}}_k^N, \quad (52)$$

where the modal forces $\hat{\mathbf{F}}_k^N$ have been introduced

$$\hat{\mathbf{F}}_k^N = \cos\left(\frac{\pi\mathbf{k}}{N}\right) \exp\left(-\frac{i\pi\mathbf{k}}{N}\right) \hat{\mathbf{b}}_k^N. \quad (53)$$

### 4.4. Optimization of the total potential energy

Gathering Eqs. (46), (50) and (52), leads to the following expression of the total potential energy $\Pi$

$$\Pi = U - V = \frac{|\mathbf{L}|}{|\mathbf{N}|^2} \sum_{k \in \mathcal{N}} \left( \frac{1}{2} \hat{\mathbf{u}}_k^N \cdot \hat{\mathbf{K}}_k^N \cdot \hat{\mathbf{u}}_k^N + \overline{\hat{\mathbf{u}}}_k^N \cdot \overline{\hat{\mathbf{\omega}}}_k^N \cdot \overline{\hat{\mathbf{B}}}_k^N - \hat{\mathbf{F}}_k^N \cdot \overline{\hat{\mathbf{u}}}_k^N \right), \quad (54)$$

the optimization of which (with respect to the modal displacements $\hat{\mathbf{u}}_k^N$) results in a set of $d \times d$ linear systems of equations (one such system for each value of $k \in \mathcal{N}$)

$$\hat{\mathbf{K}}_k^N \cdot \hat{\mathbf{u}}_k^N = \hat{\mathbf{F}}_k^N - \overline{\hat{\mathbf{\omega}}}_k^N \cdot \overline{\hat{\mathbf{B}}}_k^N. \quad (55)$$

For $k = 0$, the linear system (55) is ill-posed, since $\hat{\mathbf{K}}_0^N = 0$ [see Eq. (45)]. This is consistent with the fact that the solution to Eq. (7) is unique up to an additive constant. It is customary to impose the average displacement to be zero. Then, Eq. (55) must be replaced for $k = 0$ with

$$\hat{\mathbf{u}}_0^N = 0. \quad (56)$$

### 4.5. Summary of the method

The procedure proposed in this paper for the reconstruction of the displacement field results from the above developments (see Secs. 4.1, 4.2 and 4.3). It proceeds in four steps.

**Algorithm 1 (Reconstruction of the displacement)**
1. Compute the discrete Fourier transform $(\hat{\mathbf{\omega}}_n^N)_{n \in \mathcal{N}}$ of the cell values $(\mathbf{\omega}_n^N)_{n \in \mathcal{N}}$ of the prestress.

2. Compute the discrete Fourier transform $(\hat{\mathbf{b}}_k^N)_{k \in \mathcal{N}}$ of the cell values $(\mathbf{b}_n^N)_{n \in \mathcal{N}}$ of the body forces.

3. For each discrete frequency $k \in \mathcal{N}$

(a) Use Eq. (41) to compute the modal strain-displacement vector $\hat{\mathbf{B}}_k^N$.

(b) Use Eqs. (47), (48a) and (48b) to compute the modal stiffness matrix $\hat{\mathbf{K}}_k^N$.

(c) Use Eq. (53) to compute the modal force $\hat{\mathbf{F}}_k^N$.

(d) Solve the $d \times d$ linear system (55) [or use Eq. (56)] to find the modal displacement $\hat{\mathbf{u}}_k^N$.

4. Compute the inverse Fourier transform of the modal displacements $(\hat{\mathbf{u}}_k^N)_{k \in \mathcal{N}}$ to obtain the nodal displacements $(\mathbf{u}_n^N)_{n \in \mathcal{N}}$.

In the absence of body forces, it is observed that the overall cost of Algorithm 1 is dominated by that of $d+s$ discrete Fourier transforms, where $s = d(d+1)/2$. This should be compared to the cost of one iteration of the UGPLS solver, which in general requires $2s$ discrete Fourier transforms. Therefore, the cost of the proposed reconstruction is commensurable with that of one iteration of the UGPLS solver. Of course, the fast Fourier transform is used to compute discrete Fourier transforms; the cost of the present computation therefore grows as $|\mathbf{N}|\log|\mathbf{N}|$, where $|\mathbf{N}|$ is the total number of cells.
<figure><img src="image_4.png" /><figcaption>Figure 3. The model microstructure studied in Sec. 5.</figcaption></figure>

<table><caption>Table I. Summary of the simulations carried out on the microstructure depicted in Fig. 3. The first column is the grid size, the second column is the number of CG iterations. The third column is the direct estimate from Eq. (59) of the 1212 macroscopic modulus and the fourth column is the bound on this modulus, computed by means of the displacement reconstructed from Algorithm 1 and the minimum potential energy principle.</caption><thead><tr><th>N<sub>1</sub> = N<sub>2</sub></th><th>Iterations</th><th>C<sup>eff</sup><sub>1212</sub> from &tau;</th><th>C<sup>eff</sup><sub>1212</sub> from u</th></tr></thead><tbody><tr><td>32</td><td>97</td><td>1.41823</td><td>1.43758</td></tr><tr><td>64</td><td>159</td><td>1.41872</td><td>1.42527</td></tr><tr><td>128</td><td>190</td><td>1.41887</td><td>1.42113</td></tr><tr><td>256</td><td>197</td><td>1.41892</td><td>1.41971</td></tr><tr><td>512</td><td>192</td><td>1.41894</td><td>1.41921</td></tr><tr><td>1024</td><td>185</td><td>1.41894</td><td>1.41904</td></tr><tr><td>2048</td><td>178</td><td>1.41895</td><td>1.41898</td></tr><tr><td>4096</td><td>176</td><td>1.41895</td><td>1.41896</td></tr></tbody></table>

## 5. NUMERICAL RESULTS

This section provides an application of Algorithm 1. The method proposed in this paper applies equally well in two and three dimensions. However, for the sake of simplicity, we selected a plane strain elasticity problem. The geometry of the biphasic microstructure is defined in Fig. 3. This kind of model microstructures has already been used for illustration purposes by many authors [13, 12, 15, 27, 16]. The material constants of the two phases are

$$ \mu_i = 100\mu_m, \qquad v_m = 0.3, \qquad v_i = 0.2, \qquad (57) $$

and the unit cell is subjected to a macroscopic shear strain

$$ \mathbf{E} = E_{12} (\mathbf{e}_1 \otimes \mathbf{e}_2 + \mathbf{e}_2 \otimes \mathbf{e}_1). \qquad (58) $$

The corrector problem (1) is solved by means of a UGPLS solver combining a conjugate gradient (CG) solver [13] and the so-called filtered Green operator [12]. The elastic constants of the reference material are $\mu_0 = 0.9\mu_m$ and $v_0 = 0.3$. The tolerance for the stopping criterion of the CG solver was set to $10^{-10}$, which is extremely low. The computation was carried out on square grids of increasing size ($64 \times 64, 128 \times 128, \dots, 4096 \times 4096$), and was distributed over 8 cores for the largest simulations. It is observed (see Table I, second column) that the total number of iterations is fairly stable regardless of the grid size. Fig. 4 shows maps of the numerical estimates of the stress polarization, $\tau^N$ for $N_1 = N_2 = 32$.

For each grid-size $N$, the piecewise constant estimate of the stress polarization $\tau^N$ provides an estimate of the 1212 effective elastic modulus. Indeed, taking the volume average of Eq. (3c), we have

$$ C_{1212}^{\text{eff}} = \mu_0 + \frac{\bar{\tau}_{12}}{2E_{12}} \approx \mu_0 + \frac{\bar{\tau}_{12}^N}{2E_{12}}, \qquad (59) $$
<figure><img src="image_5.png" /><figcaption>Figure 4. Maps of the normalized components of the stress polarization: $\tau_{11}/\mu_m$ (left), $\tau_{22}/\mu_m$ (middle) and $\tau_{12}/\mu_m$ (right). Owing to the choice of the reference material, the stress polarization is almost null in the “matrix”; therefore, only the $(0, L/2) \times (0, L/2)$ region (“inclusion”) is shown. The same color scale was used for all three maps. The present images correspond to the $32 \times 32$ simulation.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5. Maps of the normalized components of the displacement: $u_1/(E_{12}L)$ (left) and $u_2/(E_{12}L)$ (right). Unlike Fig. 4, the whole unit cell $(0, L) \times (0, L)$ is represented here. The same color scale was used for the two maps (mutually symmetrical about the first bisector); the present images correspond to the $32 \times 32$ simulation. On the left image are also represented the lines along which the displacement profiles are plotted in Fig. 6.</figcaption></figure>

and the resulting estimates are shown in Table I (third column). Then, Algorithm 1 is used to reconstruct the displacement $\mathbf{u}^N$. Fig. 5 shows the resulting maps for $N_1 = N_2 = 32$.

For validation purposes, Fig. 6 compares the $32 \times 32$ UGPLS results with a standard, $256 \times 256$ displacement-based finite element approximation of the displacement, which serves as a reference. A very good agreement is observed, despite the fact that the grid used for the Lippmann-Schwinger based approach is rather coarse. It should be noted that the reference solution was computed with the deal.II finite element library [29, 30].

To close this section, the minimum energy principle is used to produce guaranteed upper bounds on the 1212 elastic modulus with the reconstructed displacements $\mathbf{u}^N$.

$$2C_{1212}^{\text{eff}} E_{12}^2 = \frac{1}{2} \mathbf{E} : \mathbf{C}^{\text{eff}} : \mathbf{E} \le \frac{1}{2} (\mathbf{E} + \nabla^s \mathbf{u}^N) : \mathbf{C} : (\mathbf{E} + \nabla^s \mathbf{u}^N). \quad (60)$$

It is therefore necessary to evaluate the strain energy of the heterogeneous microstructure for the displacement $\mathbf{u}^N$. Unlike the homogeneous case, this is carried out in the real space, using the P1 element stiffness matrix. More precisely,

$$\frac{1}{2}(\mathbf{E} + \nabla^s \mathbf{u}^N) : \mathbf{C} : (\mathbf{E} + \nabla^s \mathbf{u}^N) = \frac{1}{2|\mathbf{L}|} \sum_{\mathbf{n} \in \mathcal{N}} q_{\mathbf{n}}^{\mathrm{T}} K_{\mathbf{n}} q_{\mathbf{n}}, \quad (61)$$

where $K_n$ denotes the element stiffness matrix of element $n$, and $q_n$ is the column vector of the full nodal displacements (including the affine part: $\mathbf{E} \cdot \mathbf{x}_n^N + \mathbf{u}_n^N$). The resulting bounds are tabulated in Table I (fourth column), and plotted in Fig. 7, alongside the estimates of the same quantity derived
<figure><img src="image_7.png" /><figcaption>Figure 6. Profile of the first component of the normalized displacement, $u_1/(E_{12}L)$, along the lines $x_1 = L/4$ (squares) and $x_1 = 3L/4$ (triangles). Solid lines correspond to a reference, displacement-based FEM solution computed on a $256 \times 256$ cartesian grid with P1 elements. Symbols correspond to the displacement reconstructed from the stress polarization $\tau^N$, with $N_1 = N_2 = 32$. Both approaches are in excellent agreement.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7. Comparison of the estimates (squares) of $C_{1212}^{\text{eff}}$ based on Eq. (58) (see also Table I, third column) with the upper bounds (triangles) based on Eq. (60) (see also Table I, fourth column). As expected, both estimates and bounds seem to converge to the same value; however, the estimates are more accurate than the rigorous bounds.</figcaption></figure>

from the average stress polarization [see Eq. (59)]. It is observed that both bounds and estimate seem to converge to the same value; this was an expected result, which provides a quantitative validation of the procedure presented in this paper. It should be noted that, for the same grid size, the estimates are slightly more accurate than the rigorous bounds (the maximum relative distance between estimate and bound is about 1%).

## 6. CONCLUSION

A procedure is presented to reconstruct displacements from the output of a Lippmann-Schwinger solver on a uniform, rectangular grid with periodic boundary conditions. The reconstruction is formulated as a standard problem of elastic equilibrium of a *homogeneous* material (also with periodic boundary conditions).

This auxiliary problem is solved by means of displacement-based P1 finite elements. Owing to the uniformity of the grid and the homogeneity of the material, the global stiffness matrix exhibits a block-circulant structure, which lends itself to an efficient implementation in Fourier space. As a
result, the overall cost of the reconstruction is similar to the cost of one iteration of the Lippmann-Schwinger solver. It is emphasized that the procedure presented here ought to be regarded as a mere post-processing of the output of a Uniform Grid Periodic Lippmann-Schwinger (UGPLS) solver. As such, it applies to any UGPLS solver, regardless of the actual discrete Green operator that is used.

The method applies to both two- and three-dimensional problems. A simple, two-dimensional (plane strain) illustration is provided. The reconstructed displacements are compared to a direct calculation with displacement-based finite elements; both approaches were found to be in excellent agreement. The reconstructed displacement field is then combined with the minimum potential energy principle to produce rigorous upper bounds on the effective elastic moduli.

This work constitutes the first step towards a posteriori error estimates for UGPLS solvers, based on the error in constitutive relation. To achieve this goal, both kinematically admissible displacement field and statically admissible stress field are required. The present paper addresses the former; as for the latter, we are currently investigating how standard techniques commonly used for the aposteriori error analysis of finite element solutions can be adapted to the present case.

## REFERENCES

1. Milton GW. *The Theory of Composites*. Cambridge Monographs on Applied and Computational Mathematics, Cambridge University Press: The Edinburgh Building, Cambridge CB2 8RU, UK, 2002.

2. Dormieux L, Kondo D, Ulm F. *Microporomechanics*. John Wiley & Sons, Ltd/, 2006.

3. Auriault JL, Boutin C, Geindreau C. *Homogenization of Coupled Phenomena in Heterogenous Media*. ISTE, 2009.

4. Bensoussan A, Lions JL, George P (eds.). *Asymptotic Analysis for Periodic Structures, Studies in Mathematics and its Applications*, vol. 5. North-Holland, 1978.

5. Cioranescu D, Donato P. *An Introduction to Homogenization, Oxford Lecture Series in Mathematics and Its Applications*, vol. 17. Oxford University Press, 1999.

6. Jikov VV, Kozlov SM, Oleinik OA. *Homogenization of Differential Operators and Integral Functionals*. Springer, 1994.

7. Moulinec H, Suquet P. A fast numerical method for computing the linear and nonlinear properties of composites. *Comptes-rendus de l'Académie des sciences série II* 1994; **318**(11):1417-1423.

8. Moulinec H, Suquet P. A numerical method for computing the overall response of nonlinear composites with complex microstructure. *Computer Methods in Applied Mechanics and Engineering* 1998; **157**(1-2):69-94, doi: 10.1016/S0045-7825(97)00218-1.

9. Korringa J. *Theory of elastic constants of heterogeneous media*. *Journal of Mathematical Physics* 1973; **14**(4):509-513, doi:10.1063/1.1666346.

10. Zeller R, Dederichs PH. Elastic constants of polycrystals. *Physica Status Solidi (B)* 1973; **55**(2):831-842, doi: 10.1002/pssb.2220550241.

11. Kröner E. On the physics and mathematics of self-stresses. *Topics in Applied Continuum Mechanics*, Zeman JL, Ziegler F (eds.), Springer Verlag Wien, 1974; 22-38.

12. Brisard S, Dormieux L. Combining Galerkin approximation techniques with the principle of Hashin and Shtrikman to derive a new FFT-based numerical method for the homogenization of composites. *Computer Methods in Applied Mechanics and Engineering* 2012; **217-220**(0):197-212, doi:10.1016/j.cma.2012.01.003.

13. Brisard S, Dormieux L. FFT-based methods for the mechanics of composites: A general variational framework. *Computational Materials Science* 2010; **49**(3):663-671, doi:10.1016/j.commatsci.2010.06.009.

14. Yvonnet J. A fast method for solving microstructural problems defined by digital images: a space Lippmann-Schwinger scheme. *International Journal for Numerical Methods in Engineering* 2012; **92**(2):178-205, doi: 10.1002/nme.4334.

15. Willot F, Abdallah B, Pellegrini YP. Fourier-based schemes with modified green operator for computing the electrical response of heterogeneous media with accurate local fields. *International Journal for Numerical Methods in Engineering* 2014; **98**(7):518-533, doi:10.1002/nme.4641.

16. Willot F. Fourier-based schemes for computing the mechanical response of composites with accurate local fields. *Comptes Rendus Mécanique* 2015; **343**(3):232-245, doi:10.1016/j.crme.2014.12.005.

17. Eyre DJ, Milton GW. A fast numerical scheme for computing the response of composites using grid refinement. *European Physical Journal-Applied Physics* 1999; **6**(1):41-47, doi:10.1051/epjap:1999150. Progress in Electromagnetics Research Symposium (PIERS 98), NANTES, FRANCE, JUL 13-17, 1998.

18. Michel JC, Moulinec H, Suquet P. A computational scheme for linear and non-linear composites with arbitrary phase contrast. *International Journal for Numerical Methods in Engineering* 2001; **52**(1-2):139-160, doi:10.1002/nme.275.

19. Monchiet V, Bonnet G. A polarization-based fft iterative scheme for computing the effective properties of elastic composites with arbitrary contrast. *International Journal for Numerical Methods in Engineering* 2012; **89**(11):1419-1436, doi:10.1002/nme.3295.

20. Moulinec H, Silva F. Comparison of three accelerated fft-based schemes for computing the mechanical response of composite materials. *International Journal for Numerical Methods in Engineering* 2014; **97**(13):960-985, doi: 10.1002/nme.4614.

21. Hassen G, Gueguin M, de Buhan P. A homogenization approach for assessing the yield strength properties of stone column reinforced soils. *European Journal of Mechanics - A/Solids* 2013; **37**:266-280, doi:10.1016/j.euromechsol. 2012.07.003.
22. Gueguin M, Hassen G, de Buhan P. Numerical assessment of the macroscopic strength criterion of reinforced soils using semidefinite programming. *International Journal for Numerical Methods in Engineering* 2014; **99**(7):522–541, doi:10.1002/nme.4716.

23. Barthélémy JF, Dormieux L. A micromechanical approach to the strength criterion of drucker-prager materials reinforced by rigid inclusions. *International Journal for Numerical and Analytical Methods in Geomechanics* 2004; **28**(7-8):565–582, doi:10.1002/nag.368.

24. Ladevèze P, Leguillon D. Error estimate procedure in the finite element method and applications. *SIAM Journal on Numerical Analysis* 1983; **20**(3):485–509, doi:10.1137/0720033.

25. Ainsworth M, Oden J. A posteriori error estimation in finite element analysis. *Computer Methods in Applied Mechanics and Engineering* 1997; **142**(12):1–88, doi:10.1016/S0045-7825(96)01107-3.

26. Ladevèze P, Pelle JP. *Mastering Calculations in Linear and Nonlinear Mechanics*. Mechanical Engineering Series, Springer, 2005.

27. Vondřejc J, Zeman J, Marek I. Accurate upper-lower bounds on homogenized matrix by fft-based galerkin method. *arXiv.org* 2014; abs/1404.3614.

28. Vondřejc J, Zeman J, Marek I. An FFT-based galerkin method for homogenization of periodic media. *Computers & Mathematics with Applications* 2014; **68**(3):156–173, doi:10.1016/j.camwa.2014.05.014.

29. Bangerth W, Hartmann R, Kanschat G. deal.II – a general purpose object oriented finite element library. *ACM Trans. Math. Softw.* 2007; **33**(4):24/1–24/27, doi:10.1145/1268776.1268779.

30. Bangerth W, Heister T, Heltai L, Kanschat G, Kronbichler M, Maier M, Turcksin B, Young TD. The deal.II library, version 8.2. *Archive of Numerical Software* 2015; **3**, doi:10.11588/ans.2015.100.18031.

31. Suquet P. A simplified method for the prediction of homogenized elastic properties of composites with a periodic structure. *Comptes-rendus de l'Académie des sciences série II* 1990; **311**(7):769–774.

32. Nemat-Nasser S, Yu N, Hori M. Bounds and estimates of overall moduli of composites with periodic microstructure. *Mechanics of Materials* 1993; **15**(3):163–181, doi:10.1016/0167-6636(93)90016-K.

33. Aigner M, Ziegler GM. *Proofs from THE BOOK*. 4 edn., Springer Berlin Heidelberg, 2010, doi:10.1007/978-3-642-00856-6.

34. Rice JR. The localization of plastic deformation. *Theoretical and Applied Mechanics*, vol. 1, Koiter WT (ed.), North-Holland Publishing Co.: Delft, 1976; 207–220.

## A. ON THE PERIODIC GREEN OPERATOR FOR STRAINS

The periodic Green operator for strains $\Gamma_0$ of some arbitrary, homogeneous, reference material with elastic moduli $C_0$ was introduced in a periodic setting by e.g. Suquet [31] and Nemat-Nasser et al. [32]. It is defined as the linear operator which maps any $(\Omega$-periodic) prestress $\boldsymbol{\omega}$ to the opposite of the strain $\boldsymbol{\varepsilon}$ induced by this prestress. More precisely, $-\Gamma_0 * \boldsymbol{\omega}$ is by definition geometrically compatible with a periodic displacement $\mathbf{u}$ (that is, $\Gamma_0 * \boldsymbol{\omega} = -\nabla^s \mathbf{u}$) such that

$$ \nabla \cdot (C_0 : \nabla^s \mathbf{u} + \boldsymbol{\omega}) = 0. \qquad (62) $$

The Green operator for strains $\Gamma_0$ is conveniently expressed in the Fourier space

$$ (\Gamma_0 * \boldsymbol{\omega})(\mathbf{x}) = \sum_{\mathbf{k} \in \mathbb{Z}^d} \exp\left(2i\pi \frac{\mathbf{k}\mathbf{x}}{L}\right) \tilde{\Gamma}_0(\mathbf{k}) : \tilde{\boldsymbol{\omega}}_{\mathbf{k}}, \qquad (63) $$

where $\tilde{\boldsymbol{\omega}}_{\mathbf{k}}$ are the Fourier coefficients of the prestress $\boldsymbol{\omega}$ [see Eq. (16)]. For isotropic reference materials, the Fourier coefficients of the Green operator read for $\mathbf{k} \neq \mathbf{0}$

$$ \tilde{\Gamma}_0(\mathbf{k}) = \begin{cases} \mathbf{0} & \text{if } \mathbf{k} = \mathbf{0}, \\ \frac{1}{4\mu_0} (\delta_{ih}n_jn_l + \delta_{il}n_jn_h + \delta_{jh}n_i n_l + \delta_{jl}n_i n_h) - \frac{1}{2\mu_0(1-\nu_0)} n_i n_j n_h n_l & \text{otherwise,} \end{cases} \qquad (64) $$

with $\mathbf{n} = \mathbf{k}/\|\mathbf{k}\|$.

## B. ON THE PERIODIC P1 SHAPE FUNCTION

In the present section we derive the Fourier series expansion (32) of the periodic P1 shape function $\Phi^N$. To do so, it will prove convenient to introduce the so-called one-dimensional periodic hat function $\Lambda^N$. Assuming the (-1/2, 1/2) range has been divided into N equal segments, the periodic
<figure><img src="image_9.png" /><figcaption>Figure 8. The periodic hat function $\Lambda^N(\xi)$ defined in Appendix B. In the above figure, $N = 5$, and $h = 1/N$.</figcaption></figure>

hat function is first defined over this range as follows

$$ \Lambda^N(\xi) = \begin{cases} 1 - N|\xi| & \text{if } 0 \le N|\xi| \le 1, \\ 0 & \text{otherwise.} \end{cases} $$

it is then extended to the whole real line $\mathbb{R}$ by 1-periodicity (see Fig. 8). $\Lambda^N$ thus defined is the periodic counterpart of the standard one-dimensional Lagrangian shape function of degree one. The Fourier coefficients of $\Lambda^N$ are found by means of Eq. (16)

$$ \begin{aligned} \tilde{\Lambda}_k^N &= \int_{-1/2}^{1/2} \Lambda^N(\xi) \exp(-2i\pi k\xi) d\xi = \int_{-1/N}^{1/N} (1 - N|\xi|) \exp(-2i\pi k\xi) d\xi \\ &= 2 \int_{0}^{1/N} (1 - N\xi) \cos(2\pi k\xi) d\xi = \frac{1}{N} \mathrm{sinc}^2 \frac{\pi k}{N}, \end{aligned} \quad (65) $$

and the Fourier series expansion of $\Lambda^N$ therefore reads

$$ \Lambda^N(\xi) = \frac{1}{N} \sum_{k \in \mathbb{Z}} \mathrm{sinc}^2\left(\frac{\pi k}{N}\right) \exp(2i\pi k \xi). \quad (66) $$

Eq. (31) clearly shows that the $d$-dimensional periodic shape function $\Phi^N$ is the tensor product of $d$ one-dimensional periodic hat functions

$$ \Phi^N(\mathbf{x}) = \prod_{j=1}^{d} \Lambda^{N_j}(x_j/L_j), \quad (67) $$

and upon substitution of Eq. (66)

$$ \Phi^N(\mathbf{x}) = \prod_{j=1}^{d} \left[ \frac{1}{N_j} \sum_{k_j \in \mathbb{Z}} \mathrm{sinc}^2 \left( \frac{\pi k_j}{N_j} \right) \exp \left( 2i\pi \frac{k_j x_j}{L_j} \right) \right] = \frac{1}{|\mathbf{N}|} \sum_{\mathbf{k} \in \mathbb{Z}^d} \mathrm{sinc}^2 \left( \frac{\pi \mathbf{k}}{\mathbf{N}} \right) \exp \left( 2i\pi \frac{\mathbf{k} \mathbf{x}}{\mathbf{L}} \right), \quad (68) $$

where the tensorization rules (8), (9), (12) and (13) have been used. Thus Eq. (32) is retrieved.

### C. ON THE COMPUTATION OF SOME TRIGONOMETRIC SERIES

In this section, it will be convenient to introduce the functions $S_2(z)$, $S_3(z)$ and $S_4(z)$ defined for $z \notin \mathbb{Z}$ as the following series

$$ S_k(z) = \sum_{p \in \mathbb{Z}} (z + p)^{-k}. \quad (69) $$

Starting from the classical identity (see e.g. [33], Chap. 23)

$$ \pi \cot \pi z = \frac{1}{z} + \sum_{p=1}^{+\infty} \left( \frac{1}{z+p} + \frac{1}{z-p} \right), \quad (70) $$
we find by successive differentiation with respect to z

$$S_2(z) = -\pi \frac{d \cot \pi z}{dz} = \frac{\pi^2}{\sin^2 \pi z}, \qquad (71a)$$

$$S_3(z) = -\frac{1}{2} \frac{dS_2}{dz} = \pi^3 \frac{\cos \pi z}{\sin^3 \pi z}, \qquad (71b)$$

$$S_4(z) = -\frac{1}{3} \frac{dS_3}{dz} = \frac{\pi^4}{3} \frac{2 + \cos 2\pi z}{\sin^4 \pi z}. \qquad (71c)$$

Also, the following identity will be frequently invoked for $\mathbf{k} \in \mathcal{N}$ and $\mathbf{p} \in \mathbb{Z}^d$

$$\mathrm{sinc}\left(\frac{\pi(\mathbf{k} + \mathbf{p}\mathbf{N})}{\mathbf{N}}\right) = \pi^{-d}(-1)^{p_1+\cdots+p_d} \sin(\pi\mathbf{z})|\mathbf{z} + \mathbf{p}|^{-1}, \qquad (72)$$

where $\mathbf{z} = \mathbf{k}/\mathbf{N}$.

### C.1. Cell-averages of the displacement

In the present section, we derive the expression (36) of the DFT of the cell-averages of the discretized displacement. Combining Eqs. (28), (34) and (35), and taking advantage of the N-periodicity of the discrete Fourier transform $\hat{\mathbf{u}}_{\mathbf{k}}^{\mathbf{N}}$, we first find

$$\widehat{\mathrm{mean}}_{\mathbf{k}}^{\mathbf{N}}(\mathbf{u}^{\mathbf{N}}) = \sum_{\mathbf{p} \in \mathbb{Z}^d} \mathrm{sinc}^3\left(\frac{\pi}{\mathbf{N}}(\mathbf{k} + \mathbf{p}\mathbf{N})\right) \exp\left(\frac{i\pi}{\mathbf{N}}(\mathbf{k} + \mathbf{p}\mathbf{N})\right) \hat{\mathbf{u}}_{\mathbf{k}}^{\mathbf{N}}, \qquad (73)$$

then, introducing $\mathbf{z} = \mathbf{k}/\mathbf{N}$, and using Eqs. (69), (71b) and (72)

$$\widehat{\mathrm{mean}}_{\mathbf{k}}^{\mathbf{N}}(\mathbf{u}^{\mathbf{N}}) = \pi^{-3} \sin^3(\pi\mathbf{z}) \exp(i\pi\mathbf{z}) S_3(z_1) \cdots S_3(z_d) \hat{\mathbf{u}}_{\mathbf{k}}^{\mathbf{N}} = \cos(\pi\mathbf{z}) \exp(i\pi\mathbf{z}) \hat{\mathbf{u}}_{\mathbf{k}}^{\mathbf{N}}, \qquad (74)$$

and Eq. (36) is proved.

### C.2. Derivation of the modal strain-displacement vector

In order to retrieve the closed-form expression (41) of the modal strain-displacement vector, we start from the definition (40), and introduce $\mathbf{z} = \mathbf{k}/\mathbf{N}$

$$\begin{aligned} \hat{\mathbf{B}}_{\mathbf{k}}^{\mathbf{N}} \cdot \mathbf{e}_1 &= 2i\pi^{1-3d} \sin^3(\pi\mathbf{z}) \exp(i\pi\mathbf{z}) \sum_{\mathbf{p} \in \mathbb{Z}^d} |\mathbf{z} + \mathbf{p}|^{-3} \frac{k_1 + p_1 N_1}{L_1} \\ &= 2i\pi^{1-3d} \frac{N_1}{L_1} \sin^3(\pi\mathbf{z}) \exp(i\pi\mathbf{z}) \sum_{\mathbf{p} \in \mathbb{Z}^d} |\mathbf{z} + \mathbf{p}|^{-3} (z_1 + p_1) \\ &= 2i\pi^{1-3d} \frac{N_1}{L_1} \sin^3(\pi\mathbf{z}) \exp(i\pi\mathbf{z}) S_2(z_1) S_3(z_2) \cdots S_3(z_d) \\ &= 2i\pi^{1-3d} \frac{N_1}{L_1} \sin^3(\pi\mathbf{z}) \exp(i\pi\mathbf{z}) \frac{\pi^2}{\sin^2 \pi z_1} \pi^3 \frac{\cos \pi z_2}{\sin^3 \pi z_2} \cdots \pi^3 \frac{\cos \pi z_d}{\sin^3 \pi z_d} \\ &= 2i \frac{N_1}{L_1} \sin \pi z_1 \cos \pi z_2 \cdots \cos \pi z_d \exp(i\pi\mathbf{z}), \end{aligned} \qquad (75)$$

and the proof of Eq. (41) is complete.

### C.3. Derivation of the modal stiffness matrix

In this section, we derive Eqs. (47), (48a) and (48b) which define the modal stiffness matrix of an isotropic, linearly elastic material with shear modulus $\mu_0$ and Poisson ratio $\nu_0$. The acoustic tensor
$\mathbf{A}_0(\mathbf{q})$ of such materials reads [34]

$$ \mathbf{A}_0(\mathbf{q}) = \mu_0 q^2 \mathbf{I} + \frac{\mu_0}{1 - 2\nu_0} \mathbf{q} \otimes \mathbf{q} = \mu_0 \operatorname{tr}(\mathbf{q} \otimes \mathbf{q}) \mathbf{I} + \frac{\mu_0}{1 - 2\nu_0} \mathbf{q} \otimes \mathbf{q}, \quad (76) $$

which, upon substitution in Eq. (45), leads to Eq. (47), with

$$ \hat{\mathbf{H}}_{\mathbf{k}}^{\mathbf{N}} = 4\pi^2 \sum_{\mathbf{p} \in \mathbb{Z}^d} \mathrm{sinc}^4 \left( \frac{\pi}{N} (\mathbf{k} + \mathbf{p}N) \right) \left( \frac{\mathbf{k} + \mathbf{p}N}{L} \right) \otimes \left( \frac{\mathbf{k} + \mathbf{p}N}{L} \right). \quad (77) $$

The diagonal and off-diagonal coefficients of $\hat{\mathbf{H}}_{\mathbf{k}}^{\mathbf{N}}$ are derived below. Substituting Eq. (72) into Eq. (77)

$$ \hat{\mathbf{H}}_{\mathbf{k}}^{\mathbf{N}} = 4\pi^{2-4d} \sin^4(\pi\mathbf{z}) \sum_{\mathbf{p} \in \mathbb{Z}^d} |\mathbf{z} + \mathbf{p}|^{-4} \left(\frac{\mathbf{k} + \mathbf{p}N}{L}\right) \otimes \left(\frac{\mathbf{k} + \mathbf{p}N}{L}\right), \quad (78) $$

we successively find

$$ \begin{aligned} \mathbf{e}_1 \cdot \hat{\mathbf{H}}_{\mathbf{k}}^{\mathbf{N}} \cdot \mathbf{e}_1 &= 4\pi^{2-4d} \sin^4(\pi\mathbf{z}) \sum_{\mathbf{p} \in \mathbb{Z}^d} |\mathbf{z} + \mathbf{p}|^{-4} \left(\frac{k_1 + p_1 N_1}{L_1}\right)^2 \\ &= 4\pi^{2-4d} \left(\frac{N_1}{L_1}\right)^2 \sin^4(\pi\mathbf{z}) \sum_{\mathbf{p} \in \mathbb{Z}^d} |\mathbf{z} + \mathbf{p}|^{-4} (z_1 + p_1)^2 \\ &= 4\pi^{2-4d} \left(\frac{N_1}{L_1}\right)^2 \sin^4(\pi\mathbf{z}) S_2(z_1) S_4(z_2) \cdots S_4(z_d), \\ &= \frac{2}{3^{d-1}} \left(\frac{N_1}{L_1}\right)^2 (1 - \cos 2\pi z_1) (2 + \cos 2\pi z_2) \cdots (2 + \cos 2\pi z_d), \end{aligned} \quad (79) $$

and

$$ \begin{aligned} \mathbf{e}_1 \cdot \hat{\mathbf{H}}_{\mathbf{k}}^{\mathbf{N}} \cdot \mathbf{e}_2 &= 4\pi^{2-4d} \sin^4(\pi\mathbf{z}) \sum_{\mathbf{p} \in \mathbb{Z}^d} |\mathbf{z} + \mathbf{p}|^{-4} \left(\frac{k_1 + p_1 N_1}{L_1}\right) \left(\frac{k_2 + p_2 N_2}{L_2}\right) \\ &= 4\pi^{2-4d} \frac{N_1 N_2}{L_1 L_2} \sin^4(\pi\mathbf{z}) \sum_{\mathbf{p} \in \mathbb{Z}^d} |\mathbf{z} + \mathbf{p}|^{-4} (z_1 + p_1)(z_2 + p_2) \\ &= 4\pi^{2-4d} \frac{N_1 N_2}{L_1 L_2} \sin^4(\pi\mathbf{z}) S_3(z_1)S_3(z_2) \cdots S_4(z_d), \\ &= \frac{1}{3^{d-2}} \frac{N_1 N_2}{L_1 L_2} \sin 2\pi z_1 \sin 2\pi z_2 \cdots (2 + \cos 2\pi z_d), \end{aligned} \quad (80) $$

and the proof of Eqs. (48a) and (48b) is complete.