# Fourier-based computational micromechanics on boundary-conforming transformed grids

François Bignonnet, Martin Ladecký, Ivana Pultarová, Jan Zeman

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

## To cite this version:

François Bignonnet, Martin Ladecký, Ivana Pultarová, Jan Zeman. Fourier-based computational micromechanics on boundary-conforming transformed grids. 2026. hal-05329992

HAL Id: hal-05329992

https://hal.science/hal-05329992v1

Preprint submitted on 2 Feb 2026

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_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Fourier-based computational micromechanics on
boundary-conforming transformed grids

François Bignonnet*¹, Martin Ladecký², Ivana Pultarová³, and Jan Zeman³

¹Nantes Université, École Centrale Nantes, CNRS, GeM, UMR 6183, F-44600 Saint-Nazaire, France

²Department of Microsystems Engineering, University of Freiburg, Georges-Köhler-Allee 103, 79110 Freiburg, Germany

³Faculty of Civil Engineering, Czech Technical University in Prague, Thákurova 7, 166 29 Prague 6, Czech Republic

February 2, 2026

Abstract

Computational homogenization accelerated by Green function preconditioning with Fast Fourier transforms (FFT) is classically performed on a uniform grid, which hinders the discretization accuracy. In this work, we consider a more accurate geometry representation obtained by transforming the uniform computational grid into a boundary-conforming one. The mechanical problem is discretized using the finite element method (FEM) with isoparametric transformation of elements. Boundary adaptation can require large localized geometrical transformations of the grid, which is naturally accounted for in the FEM discretization. Rigorous bounds on the spectrum of eigenvalues of the resulting discrete system with Green preconditioner are provided. For grid transformations with projection of the nearest nodes to boundary, the modified eigenvalues correspond to eigenvectors localized at the material phases boundaries, so that the effective spectrum remains favorable for the preconditioned conjugate gradient solver. Numerical investigations confirm that the accuracy of the homogenized properties and the local fields obtained on boundary-conforming grids are greatly improved over uniform grid ones, at the expense of a moderate increase in computational cost.

**Keywords**— Computational homogenization, Fast Fourier Transform, Green function preconditioning, Finite element discretization, Grid adaptation

# 1 Introduction

## 1.1 State of the art

Computational homogenization provides the macroscopic constitutive response of a heterogeneous material from the resolution of a homogenization boundary value problem defined on a Representative Volume Element (RVE) of its microstructure. Among numerical methods used to solve the homogenization problem, methods based on the Green function and fast Fourier transform (FFT), introduced 30 years ago by Moulinec & Suquet [1, 2], have become popular (see reviews [3, 4]). The efficiency of these Green function based methods results from the combination of a matrix-free implementation made possible by the uniform discretization grid, the cheap evaluation of the Green

<footnote>*email: francois.bignonnet@univ-nantes.fr</footnote>

<footnote>This manuscript has not yet been accepted for publication.</footnote>
operator in the Fourier space using FFT, and the favorable spectrum of the discretized system of equations.

The efficiency of Green function based methods has been improved by the introduction of more efficient iterative solvers than the historical Picard iterations [5, 6, 7, 8, 9]. In the meantime, these methods have been recast in a Galerkin approximation framework, either based on the Hashin-Shtrikman variational formulation [6, 10, 11] or on the classical weak formulation which allows to reuse all the FEM technology at no compromise on efficiency [12, 13, 14, 15].

One of the drawback of the initial method of Moulinec & Suquet [1, 2] is the presence of spurious oscillations, or ringing artifacts, in solution fields. These result from Gibbs phenomena due to the trigonometric collocation method. Oscillations can be suppressed by resorting to alternative discretization strategies leading to energetically consistent or filtered Green operators [6, 10, 16]. Unfortunately, these have a higher computational cost, introduce a dependence of the results on the choice of a reference medium, and loose the rigorous status of projector at the heart of the variational formulation of [12, 14, 15]. Another way to remedy oscillations is the reinterpretation in the Fourier space of discrete differential operators used in finite difference methods [17, 18]. A suppression of oscillations by this strategy is not systematically achieved, but is possible at the expense of the introduction of staggered grid [19], or of a numerical diffusion of the fields due to the inconsistent treatment of first and second order derivatives [20, 21, 22].

More recently, a finite element method (FEM) discretization on a regular grid opened a new way to remove ringing artifacts [23, 24, 25]. This strategy has been extended to voxel meshing by a periodic pattern or higher order finite elements by Leute et al. [26]. Ladecký et al. [27] interpret the use of a FEM based Green operator as a geometrically optimal preconditioner. Besides, they show that bounds on the eigenvalues of the preconditioned system of equations are independent on discretization, which partly explains the efficiency of the resolution by Green preconditioned Krylov solvers.

A yet long-standing drawback of these methods is the need to use a regular discretization grid, with no possibility to refine the grid in regions of high solution gradients. Further, the voxel structure does not allow for a boundary-conforming discretization. It introduces a loss of accuracy and artifacts near boundaries between material phases, and compromises the potential accuracy gains of higher-order FEM discretizations. These limitations have been addressed recently by the enrichment of the discretization space at boundaries using X-FEM [28]. We investigate a different idea here, elaborating on the framework of grid transformation investigated by [29, 30, 31, 32].

Zecevic et al. [29] and Bellis & Ferrier [30] independently suggested to transform the uniform grid into an arbitrary periodic structured grid. For large strain computations, Zecevic et al. [29] distinguish three grid configurations: a reference uniform grid, an initial structured grid and a current structured grid. The motivation is to ensure a better discretization of curved boundaries. The mechanical problem is solved in the current configuration, while the reference grid is used only for the efficient application of the Green operator using FFT. To do so, a linearly viscous reference heterogeneous medium is defined in the current configuration, such that the reference problem becomes homogeneous once transformed to the reference configuration. Investigated applications include a circular inclusion for which an algebraical conformal displacement field is applied to the grid, and a Cu-Nb laminate composite with slightly curved interfaces represented by sinusoids, for which the conformal displacement at every voxel is obtained by linear interpolation of the sinusoidal interface displacements. For the circular inclusion problem, spurious stress concentrations are observed at the inclusion boundary due to the inaccuracy of the discretization of the microstructure on a uniform grid. On a boundary-conforming transformed grid, these artifacts disappear and the numerical solution is very close to the analytical solution.

Bellis & Ferrier [30] show that the homogenization problem defined in the structured physical grid can be transported to a uniform computational grid where Green preconditioning is applied. The transformed problem features transformed constitutive models and loading parameters, but can be handled by Fourier-based methods with minor modifications. Bellis & Ferrier [30] investigated
optimal transport theory to control node density in regions of interest, generating a smooth grid transformation, yet not boundary-conforming. Numerical solutions on adapted grids, as compared to ones on uniform grids with the same number of degrees of freedom, reduce the error on the effective properties. Further, solution fields are qualitatively more satisfying fields, with less oscillations, and with discontinuities being better captured. However, since the grid adaptation modifies the transformed constitutive models, the contrast in transported material properties is usually higher than before transformation, which results in a less favorable eigenvalue spectrum and a greater number of iterations to convergence of the iterative solver. Further, numerical experiments relied on the truncated Fourier series discretization, which prohibits the use of non-smooth transforms generally required to conform to boundaries. Santos-Güemes et. al [31] investigated an analytical transformation of the uniform grid to increase node density around non-singular dislocations. Increasing the number of nodes in the dislocation core region, while keeping a constant total number of grid nodes, reduces fluctuations in the stress fields.

Very recently, Zecevic et al. [32] introduce a method for the generation of boundary-conforming grids. The transformation of the grid is obtained by solving a discrete system of differential equations [33]. Boundary nodes in the uniform grid are first projected on actual boundaries between material phases. The grid is then considered as a spring network, and the stiffness of the springs parametrized as a function of the distance to boundary nodes. The boundary-conforming grid adaptation is achieved by iterating on two steps: mechanical equilibrium of the spring system and projection of boundary nodes on actual boundaries. The method was successfully applied to a variety of microstructure morphologies. As the ratio of stiffness between boundary or interior nodes increases, the density of nodes near boundary increases, and cells are more distorted. As a result, mechanical fields are generally smoother and closer to the reference high-resolution solutions. Numerical experiments indicate that the number of iterations of the FFT-based solver is generally but not systematically increased with the amount of grid distortion, which is quantified by the standard deviation of the Jacobian of the grid transformation.

## 1.2 Contributions

The present work focuses on the following objectives:

1. Introduce a transparent framework for grid adaptation using finite element discretization with isoparametric transformation of elements;

2. Provide bounds on the effective spectrum of eigenvalues of the linear operator obtained from the discretization of the problem on a transformed grid, with preconditioning by the Green function;

3. Investigate the effect of boundary-conforming grid adaptation on the quality of the numerical solution and on the number of iterations of the linear solver.

We show how finite element discretization on regular grids [24, 25, 26, 27] in combination with isoparametric element transformation provides a natural framework for formulating the grid adaptation technique. Spatial convergence estimates of FEM elucidate the impact of boundary-conforming grid adaptation. The rate of convergence of a FEM solution $u_h$ to the exact weak solution $u$ depends on the smoothness of the latter [34, 35]. In the case of a smooth stiffness field, the solution $u$ is in $H^2$ on the whole domain and thus $u_h$ converges to $u$ in the energy norm with the rate $h$. In the case of a discontinuous and piecewise constant stiffness field on subdomains with $C^2$-smooth boundaries, the solution $u$ is not in $H^2$, but the numerical solution $u_h$ obtained from boundary-conforming grid still converges to $u$ with the same rate. Then, first order Lagrange FEM with boundary-conforming elements yields $h^1$ convergence of $u_h$ to $u$ in $H^2$ norm, and $h^2$ convergence of the homogenized stiffness. Although any $C^2$-smooth boundary cannot be approximated by linear elements, this convergence rate is preserved for its piecewise linear approximation (or perturbation boundary) [36]. If the mesh
does not conform to the boundary, as for uniform grids, the rate of convergence of $u_h$ to $u$ can drop to $h^{1/2}$, and that of the homogenized stiffness to $h^1$ [35, 37]. Similarly, an $h^1$ spatial convergence of the homogenized properties has been observed for different types of FFT-based methods on uniform grids, see e.g. [10, 19, 38, 39].

The numerical examples focus on a single grid adaptation strategy: only the nodes closest to the boundaries are projected on the boundaries. This choice departs from transforming all voxels [29, 30], which can generate unfavorable eigenvalue spectrum when used together with Green pre-conditioning. Alternative and promising grid adaptation strategies [40, 32] can be devised within the present framework, but fall outside scope of paper for the sake of conciseness. The resulting grid transformation is localized and non-smooth. This prohibits the use of truncated Fourier series discretizations or alternatives based on centered or staggered finite differences, but it is well-suited for finite element discretization. In particular, elements with several nodes projected on a boundary may undergo a change of topology. For example, conforming to boundary may require that some quadrangular elements have three nodes projected on the same line, and are thus transformed to triangles. Similarly, hexahedral elements with seven nodes projected on the same plane are transformed to tetrahedrons. This feature is shared with the boundary-conforming grid transformation investigated in [32]. The discretization scheme by [18], equivalent to hexahedral finite elements with reduced integration [24], is used in [32] to capture this feature. On the contrary, the types of grid transformation investigated in [30, 31] are not *exactly* boundary-conforming: although large element transformation can arise due to node densification near boundaries, the topology of voxels is not modified.

The article is organized as follows. Section 2 recalls the principles of the FEM discretization of the homogenization problem on uniform grids, and introduces its straightforward extension to transformed, structured grids. Section 3 presents the resolution of the discretized problem using the preconditioned conjugate gradient, with the Green function on a uniform grid as the preconditioner inverse. Bounds on the spectrum of eigenvalues of the preconditioned system obtained with Green preconditioning and grid adaptation are provided to quantify the effect of the localized node relocation. Section 4 details the implementation of the method. Section 5 provides numerical investigations on the accuracy and efficiency of the numerical scheme for both simple examples with analytical solutions and a more complex one with a random microstructure.

**Notations** Scalars are denoted $a$, first order tensors $\mathbf{a}$, second-order tensors $\mathbf{A}$ and fourth order tensors $\mathbf{A}$ (or $\mathbf{\Gamma}$). The double dot product of two second-order tensors $\mathbf{A}$ and $\mathbf{B}$ is $\mathbf{A} : \mathbf{B} = A_{ij}B_{ij}$. In discretized systems of equations, a vector of variables is denoted $\mathbf{a}$ and a matrix $\mathbf{A}$. A tuple of $d$ values is denoted $\mathbf{a}$ or $\mathbf{A}$ and $\mathbf{a} + \mathbf{b}$, $\mathbf{a} - \mathbf{b}$, $\mathbf{ab}$, $\mathbf{a}/\mathbf{b}$ denote element-wise operations.

## 2 Finite element discretization of homogenization problems on structured grids

A periodic unit-cell $\Omega$ of the material occupies the $d$-dimensional real space $\mathbb{R}^d$ ($d = 2, 3$). The unit-cell is denoted by $\Omega = [0, L_1] \times \cdots \times [0, L_d]$ where $L_i > 0$ are the side lengths of $\Omega$. Let $\mathcal{V} = \{\mathbf{u} \in H^1(\Omega) \mid \Omega \text{ - periodic}, \bar{\mathbf{u}} = 0\}$ stand for the set of periodic, first-order vector fields with null average over the RVE $\Omega$ belonging to the first order Sobolev space $H^1(\Omega)$. The volume average of a field $\mathbf{f}$ over $\Omega$ is denoted:

$$ \bar{\mathbf{f}} = \frac{1}{|\Omega|} \int_{\Omega} \mathbf{f} \, d\mathbf{x}. \qquad (1) $$

The set $L^2_{\text{sym}}(\Omega)$ denotes the space of symmetric second-order tensor fields, which are $\Omega$-periodic and square integrable over $\Omega$.
## 2.1 Homogenization problem

The homogenization problem for small strain linear elasticity consists in finding a microscopic fluctuation displacement field $u \in \mathcal{V}$ and a Cauchy stress field $\sigma \in L_{\text{sym}}^2(\Omega)$ such that

$$ 
\begin{align} 
\operatorname{div} \sigma &= 0 && \text{in } \Omega, \nonumber \\ 
\sigma &= \mathbb{C} : \epsilon && \text{in } \Omega, \nonumber \\ 
\epsilon &= \mathbf{E} + \nabla^{\text{sym}} u && \text{in } \Omega, \tag{2} 
\end{align} 
$$

where $\epsilon \in L_{\text{sym}}^2(\Omega)$ is the linearized strain, and $\mathbf{E}$ the macroscopic strain. The stiffness tensor field $\mathbb{C}$ is heterogeneous. The homogenization problem (2) can be reformulated as find $u \in \mathcal{V}$ such that

$$ \mathbf{K}u = \mathbf{f}_E \quad \text{where} \quad \mathbf{K} = -\operatorname{div} \mathbb{C} : \nabla^{\text{sym}} \text{ and } \mathbf{f}_E = \operatorname{div} \mathbb{C} : \mathbf{E}. \tag{3} $$

The weak formulation of problem (3) is

$$ \forall v \in \mathcal{V}, \quad \int_{\Omega} \nabla^{\text{sym}} v : \mathbb{C} : (\mathbf{E} + \nabla^{\text{sym}} u) \,dx = 0. \tag{4} $$

Problem (3) will be compared to a reference problem find $u \in \mathcal{V}$ such that

$$ \mathbf{K}_0 u = \mathbf{f} \quad \text{where} \quad \mathbf{K}_0 = -\operatorname{div} \mathbb{C}_0 : \nabla^{\text{sym}}. \tag{5} $$

In the above reference problem, $\mathbb{C}_0$ is a uniform reference stiffness and $\mathbf{f}$ an arbitrary body force field satisfying $\tilde{\mathbf{f}} = 0$. The solution to the reference problem (5) is formally expressed as a convolution with the Green function $G_0$

$$ u = G_0 * \mathbf{f} \quad \text{where} \quad G_0 = \mathbf{K}_0^{-1} \tag{6} $$

After suitable discretization, the Green function can be represented by a matrix and will be used as a preconditioner to solve the homogenization problem (2).

## 2.2 Finite element discretization on uniform grids

The finite element discretization is first presented on uniform grids, i.e. without grid adaptation. The presentation is adapted from [24, 26, 27], with a focus on details suitable for an operator based implementation and Green preconditioning. The modifications required to account for grid adaptation are presented separately in section 2.3.

To avoid proliferation of notations, the presentation is restricted to grids made up of periodic patterns comprising either:

* a single first order tensor product Lagrange finite element $Q_1$ [41], i.e. a quadrangular ($d=2$) or hexahedral ($d=3$) element,

* several first order simplicial Lagrange finite elements $P_1$ [41], i.e. triangular ($d=2$) or tetrahedral ($d=3$) elements.

For a generalization to periodic patterns made up of several elements of order $k$, the reader is referred to [27].

**Periodic pattern** The unit cell $\Omega$ can be decomposed in periodic reproductions $\mathcal{P}_p$ of a reference pattern $\mathcal{P}$, identified to $\mathcal{P}_{(0,...,0)}$. The Cartesian index $p \in \mathcal{I}(\mathbb{N}) = \{0...N_1-1\} \times \dots \times \{0...N_d-1\}$ denotes the pattern position in the grid. The grid consists in $N_p = N_1 \times \dots \times N_d$ patterns.

The reference pattern features $N_e$ reference elements $\mathcal{E}^e$, where $e \in \{1, \dots, N_e\}$ is the element index in the pattern, see Fig. 1. Nodes of the reference pattern are separated in two types. There
<figure><img src="image_3.png" /><figcaption>Figure 1: Finite element discretization: examples of pattern $\mathcal{P}_p$ made up of one quadrangular $Q_1$ element with (a) $N_q = 4$ quadrature points (b) $N_q = 1$ quadrature point or (c) two triangular $P_1$ elements and $N_q = 2$ quadrature points.</figcaption></figure>

is a fundamental node at position $\mathbf{y}^n = \mathbf{0} \in \mathcal{P}$. There also are periodic reproductions of the fundamental node $\mathbf{y}^n$ by translations $\mathbf{sh}$ where shifts $\mathbf{s} \in \mathcal{S} = \{0, 1\}^d$ and $\mathbf{h}$ is the stride of the pattern (equivalent to pixel/voxel edge lengths). There are $N_q$ quadrature points at positions $\mathbf{y}^q \in \mathcal{P}$ where $q \in \{1, \dots, N_q\}$ is the quadrature point index in the pattern.

The periodic reproduction in pattern $\mathcal{P}_p$ of the reference node are $\mathbf{x}_p^n = \mathbf{y}^n + \mathbf{ph}$ and similarly for quadrature points. There are $N_p$ nodes, $N_p N_q$ elements and $N_p N_q$ quadrature points in the grid. Nodal values of the displacement are denoted $\mathbf{u} = \mathbf{u}(\mathbf{x}_p^n)$ and can be assembled in a vector of degrees of freedom $\mathbf{u} = \{\mathbf{u}_p\}_{p \in \mathcal{I}(\mathcal{N})}$ containing $dN_p$ scalars. Similarly, stress at quadrature points are denoted $\boldsymbol{\sigma}_p^q = \boldsymbol{\sigma}(\mathbf{x}_p^q)$ and assembled in a vector $\boldsymbol{\sigma} = \{\boldsymbol{\sigma}_p^q\}_{p \in \mathcal{I}(\mathcal{N}), q \in \{1, \dots, N_q\}}$ containing $d_* N_p N_q$ scalars where $d_* = d(d+1)/2$.

**Periodic shifts** Periodic shifts $\mathbf{s} \in \mathbb{Z}^d$ of Cartesian grid indices are introduced to refer to a pattern $\mathcal{P}_{p+s}$ relatively to a pattern $\mathcal{P}_p$. The sum $\mathbf{p} + \mathbf{s}$ must be understood in a periodic way, i.e. $\mathbf{p} + \mathbf{s}$ is a shorthand notation for $(\mathbf{p} + \mathbf{s})\%\mathcal{N}$, which belongs to $\mathcal{I}(\mathcal{N})$, where % is the element-wise modulo operator. In particular if $\mathbf{s} \in \{-1, 0, 1\}^d \setminus \{0, \dots, 0\}$, the pattern $\mathcal{P}_{p+s}$ or one of its periodic reproduction is adjacent to $\mathcal{P}_p$. These shifts will be used to define discrete differential operators.

**Interpolation** At a point $\mathbf{y}$ of the reference element $\mathcal{E}^e$, displacements are interpolated from nodal values by interpolation functions $\phi_s^e$ such that $\phi_s^e(\mathbf{y}^n + \mathbf{r}\mathbf{h}) = 1$ if $\mathbf{r} = \mathbf{s}$ else 0 for all shifts $\mathbf{r}, \mathbf{s} \in \mathcal{S}$. At a point $\mathbf{x}$ of an element $\mathcal{E}^e$, the displacements are interpolated as:

$$ \mathbf{u}(\mathbf{x}) = \sum_{\mathbf{s} \in \mathcal{S}} \mathbf{u}_{p+s} \phi_s^e(\mathbf{y}) \quad \text{where} \quad \mathbf{y} = \mathbf{x} - \mathbf{p}\mathbf{h} \in \mathcal{E}^e, \qquad (7) $$

where $\mathcal{E}^e$ is a reference element. The gradient of the displacement is defined in the interior of the element $\mathcal{E}^e$ as:

$$ \nabla \mathbf{u}(\mathbf{x}) = \sum_{\mathbf{s} \in \mathcal{S}} \mathbf{u}_{p+s} \otimes \nabla \phi_s^e(\mathbf{y}). \qquad (8) $$

**Discrete gradient** The gradient of displacement at quadrature points $\mathbf{x}_p^q$ in $\mathcal{E}^e$ is then

$$ \nabla \mathbf{u}_p^q = \sum_{\mathbf{s} \in \mathcal{S}} \mathbf{u}_{p+s} \otimes \mathbf{B}_s^q \quad \text{where} \quad \mathbf{B}_s^q = \nabla \phi_s^e(\mathbf{y}^q). \qquad (9) $$
All gradients at quadrature points are gathered in a vector $\nabla u = \{\nabla u_p^q\}_{p \in \mathcal{I}(N), q \in \{1, \dots, N_q\}}$ of size $d^2 N_p N_q$. The discrete gradient operator $\hat{\mathbf{B}}$ of size $d^2 N_p N_q \times dN_p$, obtained by assembly of the $\mathbf{B}_s^q$, is defined such that $\nabla u = \hat{\mathbf{B}}u$. Similarly, the operator $\mathbf{B}$ of size $d_* N_p N_q \times dN_p$ of the symmetric gradient is such that $\nabla^{\text{sym}} u = \mathbf{B}u$.

**Discrete divergence** The divergence is the negative value of the adjoint of the gradient operator, by the application of the Gauss theorem:

$$ \forall \mathbf{u} \in \mathcal{V}, \quad \int_{\Omega} \nabla \mathbf{u} : \sigma \, d\mathbf{x} = - \int_{\Omega} \mathbf{u} \operatorname{div} \sigma \, d\mathbf{x}. \qquad (10) $$

The discrete approximation by Gauss quadrature of the l.h.s. of (10) is for any test function $\mathbf{u}$ interpolated from test nodal values $\mathbf{u}$:

$$ \int_{\Omega} \nabla \mathbf{u} : \sigma \, d\mathbf{x} \approx \sum_{\mathbf{p} \in \mathcal{I}(N)} \sum_{q=1}^{N_q} w_q \nabla \mathbf{u}_p^q : \sigma_p^q, \qquad (11) $$

where $w_q$ are quadrature weights and $\nabla \mathbf{u}_p^q$ is defined in (9). Then $\nabla \mathbf{u}_p^q : \sigma_p^q = \mathbf{u}_{p+s} \cdot \sigma_p^q \cdot \mathbf{B}_s^q$ by symmetry of $\sigma_p^q$. Reorganization of the sum on $\mathbf{p}$ by the shifted Cartesian index $\mathbf{p}' = \mathbf{p} + \mathbf{s}$ yields

$$ \sum_{\mathbf{p}' \in \mathcal{I}(N)} \sum_{q=1}^{N_q} w_q \sum_{\mathbf{s} \in \mathcal{S}} \mathbf{u}_{\mathbf{p}'} \cdot \sigma_{\mathbf{p}'-\mathbf{s}}^q \cdot \mathbf{B}_{\mathbf{s}}^q := - \sum_{\mathbf{p}' \in \mathcal{I}(N)} \mathbf{u}_{\mathbf{p}'} \cdot \operatorname{div} \sigma_{\mathbf{p}'}, \qquad (12) $$

where one defines, by identification with the r.h.s. of (10):

$$ \operatorname{div} \sigma_{\mathbf{p}} := \sum_{q=1}^{N_q} w_q \operatorname{div} \sigma_{\mathbf{p}}^q \quad \text{where} \quad \operatorname{div} \sigma_{\mathbf{p}}^q := - \sum_{\mathbf{s} \in \mathcal{S}} \sigma_{\mathbf{p}-\mathbf{s}}^q \cdot \mathbf{B}_{\mathbf{s}}^q. \qquad (13) $$

The contribution of neighbouring quadrature points $q$ to the divergence $\operatorname{div} \sigma_{\mathbf{p}}$ at the fundamental node of pattern $\mathbf{p}$ is $\operatorname{div} \sigma_{\mathbf{p}}^q$. It is the counterpart of $\nabla \mathbf{u}_{\mathbf{p}}^q$ defined in (9). Finally, gathering in a vector $\operatorname{div} \sigma = \{\operatorname{div} \sigma_{\mathbf{p}}\}_{\mathbf{p} \in \mathcal{I}(N)}$ of size $dN_p$, the discrete divergence operator (for symmetric second order tensors) stemming from the finite element discretization is $-\mathbf{B}^\mathrm{T}\mathbf{W}$ such that $\operatorname{div} \sigma = -\mathbf{B}^\mathrm{T}\mathbf{W}\sigma$, where $\mathbf{W}$ is a diagonal $d_* N_p N_q \times d_* N_p N_q$ matrix built from the quadrature weights $w_q$.

**Discrete problem** The weak formulation (4) of the homogenization problem (2), combined with the finite element discretization and Gauss quadrature approximation (11) leads to the discrete version of (3):

$$ \mathbf{K} \mathbf{u} = \mathbf{f}_E \quad \text{where} \quad \mathbf{K} = \mathbf{B}^\mathrm{T} \mathbf{W} \mathbf{C} \mathbf{B} \quad \text{and} \quad \mathbf{f}_E = - \mathbf{B}^\mathrm{T} \mathbf{W} \mathbf{C} \mathbf{E}, \qquad (14) $$

where $\mathbf{C}$ is the $d_* N_p N_q \times d_* N_p N_q$ block-diagonal matrix containing the stiffness at quadrature points, $\mathbf{K}$ is the $dN_p \times dN_p$ stiffness matrix, $\mathbf{E}$ is the $d_* N_p N_q$ vector built from the uniform the macroscopic strain $\mathbf{E}$ and $\mathbf{f}_E$ the vector of generalized forces.

## 2.3 Finite element discretization on structured grids

**Grid transformation** The geometrical transformation from the uniform grid to the structured grid is performed by the mapping $\mathbf{x} = \mathbf{\Phi}(\tilde{\mathbf{x}})$ from the computational domain $\tilde{\Omega}$ to the physical domain $\Omega$, see Fig. 2. The mapping displacement is $\tilde{\mathbf{u}}_{\mathbf{\Phi}}(\tilde{\mathbf{x}}) = \mathbf{\Phi}(\tilde{\mathbf{x}}) - \tilde{\mathbf{x}}$. This displacement belongs to $H^1(\Omega)$ and is $\Omega$-periodic. The fluctuation displacement $\mathbf{u}(\mathbf{x})$ in the physical domain is mapped from $\tilde{\mathbf{u}}(\tilde{\mathbf{x}}) = \mathbf{u}(\mathbf{\Phi}(\tilde{\mathbf{x}}))$ in $\tilde{\Omega}$. The displacement $\tilde{\mathbf{u}} \in \mathcal{V}$, defined as $\mathcal{V}$ with the modification $\int_{\tilde{\Omega}} \tilde{\mathbf{u}} \, d\tilde{\mathbf{x}} = 0$.
<figure><img src="image_4.png" /><figcaption>Figure 2: Transformation of a reference uniform grid to a boundary-conforming grid by a mapping $\boldsymbol{x} = \boldsymbol{\Phi}(\tilde{\boldsymbol{x}})$.</figcaption></figure>

The gradient of the displacement $\nabla \boldsymbol{u} = \partial \boldsymbol{u}/\partial \boldsymbol{x}$ in the physical domain $\Omega$ is related to the gradient $\tilde{\nabla} \tilde{\boldsymbol{u}} = \partial \tilde{\boldsymbol{u}}/\partial \tilde{\boldsymbol{x}}$ in the computational domain $\tilde{\Omega}$ by chain rule

$$ \tilde{\nabla} \tilde{\boldsymbol{u}} = \frac{\partial \tilde{\boldsymbol{u}}}{\partial \tilde{\boldsymbol{x}}} = \nabla \boldsymbol{u} \cdot \boldsymbol{F} \quad \text{where} \quad \boldsymbol{F} = \tilde{\nabla} \boldsymbol{\Phi} \quad ; \quad J = \det(\boldsymbol{F}). \qquad (15) $$

The weak formulation (4) of the homogenization problem in the physical domain $\Omega$ is pulled back to the reference domain $\tilde{\Omega}$:

$$ \forall \tilde{\boldsymbol{v}} \in \tilde{\mathcal{V}}, \int_{\tilde{\Omega}} \mathrm{sym}(\tilde{\nabla} \tilde{\boldsymbol{v}} \cdot \boldsymbol{F}^{-1}) : \mathbb{C} : (\boldsymbol{E} + \mathrm{sym}(\tilde{\nabla} \tilde{\boldsymbol{u}} \cdot \boldsymbol{F}^{-1})) J \, d\tilde{\boldsymbol{x}} = 0. \qquad (16) $$

**Isoparametric transformation** We restrict the presentation to an isoparametric transformation of the uniform grid [41]. Each node $\tilde{x}_p^n \in \tilde{\Omega}$ of the uniform grid is transported at its physical position $x_p^n = \tilde{x}_p^n + \tilde{u}_\Phi(\tilde{x}_p^n)$. The periodic nodal mapping displacements are stored in a vector $\tilde{\boldsymbol{u}}_\Phi$. A point $\tilde{x}$ in the element $\tilde{\mathcal{E}}_p^e$ is transported at $\boldsymbol{x}$ in the element $\mathcal{E}_p^e$ by interpolation as in (7)

$$ \boldsymbol{x} = \tilde{\boldsymbol{x}} + \tilde{u}_{\Phi}(\tilde{\boldsymbol{x}}) = \tilde{\boldsymbol{x}} + \sum_{\boldsymbol{s} \in \mathcal{S}} (\tilde{u}_{\Phi})_{\boldsymbol{p}+\boldsymbol{s}} \phi_{\boldsymbol{s}}^{\boldsymbol{e}}(\boldsymbol{y}) \quad \text{where} \quad \boldsymbol{y} = \tilde{\boldsymbol{x}} - \boldsymbol{p}\boldsymbol{h} \in \tilde{\mathcal{E}}^{\boldsymbol{e}}. \qquad (17) $$

The transformation gradient at quadrature point is obtained from $\tilde{\nabla} \tilde{\boldsymbol{u}}_{\Phi} = \tilde{\boldsymbol{B}} \tilde{\boldsymbol{u}}_{\Phi}$ following the definition of the discrete gradient operator (9) on the uniform grid as

$$ \boldsymbol{F}_{\boldsymbol{p}}^{\boldsymbol{q}} = \boldsymbol{I} + (\tilde{\nabla} \tilde{\boldsymbol{u}}_{\Phi})_{\boldsymbol{p}}^{\boldsymbol{q}}. \qquad (18) $$

The fluctuation displacement $\boldsymbol{u}$ is interpolated from nodal values $\tilde{\boldsymbol{u}}$ in the computational grid by

$$ \boldsymbol{u}(\boldsymbol{x}) = \tilde{\boldsymbol{u}}(\tilde{\boldsymbol{x}}) = \sum_{\boldsymbol{s} \in \mathcal{S}} \tilde{\boldsymbol{u}}_{\boldsymbol{p}+\boldsymbol{s}} \phi_{\boldsymbol{s}}^{\boldsymbol{e}}(\boldsymbol{y}) \quad \text{where} \quad \boldsymbol{y} = \tilde{\boldsymbol{x}} - \boldsymbol{p}\boldsymbol{h} \in \tilde{\mathcal{E}}^{\boldsymbol{e}}. \qquad (19) $$
Eventually, at quadrature points $x_p^q = \tilde{x}_p^q + \tilde{u}_\Phi(\tilde{x}_p^q)$, the gradient of the displacement follows from (15) as

$$ (\nabla \mathbf{u})_p^q = \tilde{\nabla} \tilde{u}_p^q \cdot (\mathbf{F}_p^q)^{-1}, \qquad (20) $$

where $\tilde{\nabla} \tilde{u} = \tilde{\mathbf{B}} \tilde{u}$ from (9).

**Transformed discrete problem** The Gauss quadrature approximation (11) of the pulled-back weak formulation (16) leads to the discrete version of (3) for a transformed grid

$$ \tilde{\mathbf{K}}\tilde{\mathbf{u}} = \tilde{\mathbf{f}}_E \quad \text{where} \quad \tilde{\mathbf{K}} = \tilde{\mathbf{B}}^T \tilde{\mathbf{W}} \tilde{\mathbf{C}} \tilde{\mathbf{B}} \quad \text{and} \quad \tilde{\mathbf{f}}_E = -\tilde{\mathbf{B}}^T \tilde{\mathbf{W}} \tilde{\mathbf{C}} \tilde{\mathbf{E}}, \qquad (21) $$

where $\tilde{\mathbf{W}}$ is the $d^2 N_p N_q \times d^2 N_p N_q$ diagonal matrix built from the same quadrature weights $w_q$ as for a uniform grid. Note that a quadrature approximation which is exact on the uniform grid (for example $N_q = 8$ on a trilinear hexahedral element) can result in an inexact integration on the transformed grid.

In (21), $\tilde{\mathbf{C}}$ is the $d^2 N_p N_q \times d^2 N_p N_q$ block diagonal matrix built from the pulled-back stiffness tensor $\tilde{\mathbf{C}}$ defined as [30]:

$$ \tilde{\mathbf{C}}_{ijkl} = J F_{jp}^{-1} C_{pikq} F_{ql}^{-T}, \qquad (22) $$

where $\mathbf{F}$ and $J = \det(\mathbf{F})$ are obtained at each Gauss integration point from (18). The fourth order tensor $\tilde{\mathbf{C}}$ is positive definite, possesses the major symmetry but not the minor symmetries, and is expressed in the Voigt notation as a $d^2 \times d^2$ matrix.

The r.h.s. is formally constructed in (21) from the pulled-back stiffness and the pulled-back macroscopic strain $\tilde{\mathbf{E}}$ defined at quadrature points as the second order tensor

$$ \tilde{\mathbf{E}}_p^q = \mathbf{E} \cdot \mathbf{F}_p^q. \qquad (23) $$

The strain $\tilde{\mathbf{E}}$ is neither symmetric nor uniform. In practice, the explicit computation of $\tilde{\mathbf{C}}$ and $\tilde{\mathbf{E}}$ can be avoided by using an operator based implementation as detailed in section 4.

## 2.4 Differential operators and Green function in the Fourier space

**Discrete Fourier Transform (DFT)** The DFT of the array $\mathbf{u}$ at Fourier index $k \in \mathcal{I}(\mathbf{N})$ is defined by:

$$ \hat{\mathbf{u}}_k = \frac{1}{N_p} \sum_{\mathbf{p} \in \mathcal{I}(\mathbf{N})} \mathbf{u}_p \exp(-i2\pi\mathbf{k} \cdot \mathbf{p}/N). \qquad (24) $$

In practice, symmetry properties of the real-to-complex DFT allow restricting the set of Fourier indices to $\tilde{\mathcal{I}}(\mathbf{N}) = \{0, \dots, N_1 \div 2\} \times \{0, \dots, N_2 - 1\} \times \{0, \dots, N_3 - 1\}$ where $\div$ denotes integer division<fnref n="1" />.

**DFT of discrete gradient** From (9), the DFT of the array of discrete gradient $\{\nabla \mathbf{u}_p^q\}_{\mathbf{p} \in \mathcal{I}(\mathbf{N})}$ is consequently

$$ \widehat{\nabla} \mathbf{u}_k^q = \frac{1}{N_p} \sum_{\mathbf{p} \in \mathcal{I}(\mathbf{N})} \sum_{\mathbf{s} \in \mathcal{S}} \mathbf{u}_{\mathbf{p}+\mathbf{s}} \otimes \mathbf{B}_s^q \exp(-i2\pi\mathbf{k} \cdot \mathbf{p}/N). \qquad (25) $$

Application of the DFT shift theorem on Cartesian index $\mathbf{p}$ for shift $\mathbf{s}$ leads to the expression of the discrete gradient in the Fourier domain:

$$ \widehat{\nabla} \mathbf{u}_k^q = \hat{\mathbf{u}}_k \otimes i \xi_k^q \quad \text{and} \quad i \xi_k^q = \sum_{\mathbf{s} \in \mathcal{S}} \mathbf{B}_s^q \exp (i2\pi \mathbf{k} \cdot \mathbf{s}/N). \qquad (26) $$

By direct analogy with the continuous Fourier transform, the interpretation of $\xi_k^q$ is the discrete wave vector associated to the discrete gradient $\nabla \mathbf{u}_p^q$ defined in (9).

<footnote n="1">¹Depending on the programming language, the convention for memory order of multi-dimensional arrays can change, and division by 2 can be on the first or last dimension.</footnote>
**DFT of discrete divergence** Similarly, the DFT of the array div σ defined in (13) is

$$ \widehat{\operatorname{div} \sigma_k} = \sum_{q=1}^{N_q} w_q \widehat{\operatorname{div} \sigma_k^q} \quad \text{where} \quad \widehat{\operatorname{div} \sigma_k^q} = \hat{\sigma}_k^q \cdot \iota \bar{\xi}_k^q, \qquad (27) $$

where $\bar{\xi}_k^q$ is the complex conjugate of the discrete wave vector $\xi_k^q$ defined in (26).

**Discrete Green function** In the reference problem (5), the heterogeneous material properties $\mathbb{C}$ are replaced by a homogeneous one $\mathbb{C}_0$. The discretization leads to the reference stiffness matrix

$$ \mathbf{K}_0 = \mathbf{B}^T \mathbf{W} \mathbf{C}_0 \mathbf{B}. \qquad (28) $$

By construction, $\mathbf{K}_0$ is block-circulant, and its DFT is derived from the discrete wave vectors defined in (26). The discrete Green operator $\mathbf{G}_0$ is the (pseudo-)inverse of $\mathbf{K}_0$, see (6). The DFT of the discrete Green function can then be computed as follows.

**DFT of discrete Green function** Let the nodal force vector $\mathbf{f}_0 = \mathbf{K}_0 \mathbf{u}$. From the expression of the DFT of the discrete differential operators (26) and (27), the DFT of $\mathbf{K}_0$ operates diagonally at Fourier index $k$ as

$$ \hat{\mathbf{f}}_{0k} = \hat{\mathbf{K}}_{0k} \mathbf{u}_k = - \sum_{q=1}^{N_q} w_q (\mathbb{C}_0 : (\hat{\mathbf{u}}_k \otimes \iota \xi_k^q)) \cdot \iota \bar{\xi}_k^q = \hat{\mathbf{K}}_{0k} \cdot \hat{\mathbf{u}}_k, \qquad (29) $$

where $\hat{\mathbf{K}}_{0k}$ is the acoustic tensor defined by

$$ \hat{\mathbf{K}}_{0k} = \sum_{q=1}^{N_q} w_q (\bar{\xi}_k^q \cdot \mathbb{C}_0 \cdot \xi_k^q). \qquad (30) $$

When $\mathbb{C}_0$ is isotropic with Lamé parameters $(\lambda_0, \mu_0)$

$$ \bar{\xi}_k^q \cdot \mathbb{C}_0 \cdot \xi_k^q = \lambda_0 \bar{\xi}_k^q \otimes \xi_k^q + \mu_0 (\bar{\xi}_k^q \cdot \xi_k^q \mathbf{I} + \xi_k^q \otimes \bar{\xi}_k^q). \qquad (31) $$

The DFT of the discrete Green function is then the real and symmetric second order tensor

$$ \hat{\mathbf{G}}_{0k} = \hat{\mathbf{K}}_{0k}^{-1} \quad \text{if } k \neq (0, \dots, 0) \text{ else } 0, \qquad (32) $$

such that $\hat{\mathbf{u}}_k = \hat{\mathbf{G}}_{0k} \hat{\mathbf{f}}_{0k}$.

When $N_q = 1$ (e.g., $Q_1$ type finite elements with reduced integration [18, 24]), the inverse of the acoustic tensor $\hat{\mathbf{K}}_{0k}$ is known analytically from the discrete wave vector $\xi_k = \xi_k^1$ as

$$ \hat{\mathbf{G}}_{0k} = \frac{1}{\mu_0 \|\xi_k\|^2} \left( \mathbf{I} - \frac{\lambda_0 + \mu_0 \bar{\xi}_k \otimes \xi_k}{\lambda_0 + 2\mu_0 \|\xi_k\|^2} \right) \quad \text{if } \xi_k \neq 0 \text{ else } 0, \qquad (33) $$

and $\bar{\xi}_k \otimes \xi_k$ is real and symmetric. Otherwise, the inversion of a system of $d$ linear equations must be carried out numerically at each Fourier index $k \in \hat{\mathcal{I}}(\mathbb{N})$.
# 3 Fast iterative solvers with Green preconditioning

## 3.1 Preconditioned conjugate gradient method

The homogenization problem discretized on either the uniform or transformed grid leads to the stiffness matrices $K$ (14) or $\tilde{K}$ (21) which are symmetric and positive semi-definite. Modes with eigenvalue zero correspond to rigid body translations, and are not present in the right hand side. The problems (14) or (21) can thus be solved using conjugate gradient (CG) [42, 43, 44], see e.g. [5, 45, 19, 8], or fast gradient solvers [46, 9].

To simplify the discussion of this subsection, $K$ denotes indifferently the stiffness matrix on the uniform or transformed grid in what follows. The eigenvalue spectrum of $K$ is spread, non clustered, and the condition number increases quadratically with grid refinement [41], so that the preconditioned conjugate gradient (PCG) is recommended over CG [47]. The preconditioner should be a matrix $M$ sufficiently close to $K$ to improve convergence, yet easily invertible in the row space of $K$. Green preconditioning refers to the choice $M = K_0$ defined by (28), whose pseudo-inverse in the row space of $K$ is the discrete Green function $G_0$.

The convergence of PCG is governed by the eigenvalue spectrum of the preconditioned system together with the distribution of the modal coordinates of the initial residual. Let $\lambda_i$ and $z_i$ be the $i^{th}$ eigenvalue and eigenvector of the generalized eigenvalue problem for matrix $K$ with preconditioner $K_0$ (restricted to the subspace orthogonal to the shared kernel of $K$ and $K_0$), i.e. $Kz_i = \lambda_i K_0 z_i$. Eigenvectors are normalized such that $z_i^T K_0 z_j = \delta_{ij}$. The modal weight of the initial residual $r^{(0)} = f_E - Ku^{(0)}$ is defined as

$$w_i = \frac{(r^{(0)} \cdot z_i)^2}{\lambda_i} \qquad (34)$$

The initial residual is equal to the r.h.s. $f_E$ for the initial guess $u^{(0)} = 0$. At the $k$-th iteration of the PCG, the error on the solution $u$ of the iterate $u^{(k)}$ is

$$\|u^{(k)} - u\|_{K} = \|r^{(k)}\|_{K^{-1}} = \sum_{i=1}^{dN_p} w_i [P^{(k)}(\lambda_i)]^2. \qquad (35)$$

where $P^{(k)}(\lambda)$ is the polynomial of degree $k$ satisfying $P^{(k)}(0) = 1$ which has the least value of the above error function [48, 49, 50]. Then convergence of PCG depends on the:

* condition number $\kappa = \lambda_{max}/\lambda_{min}$;

* clustering of eigenvalues with non negligible modal weights of the first residual.

Throughout the paper, the iterative solution is computed with the PCG and the first residual $r^{(0)}$ is set equal to the r.h.s vector. Iterations are run until the (too) stringent stopping criterion $\|r^{(k)}\|_{G_0} < 10^{-10} \|r^{(0)}\|_{G_0}$ is met. The preconditioner norm of the residual $\|r^{(k)}\|_{G_0}$ is required by the PCG algorithm to update the search direction, thus comes at no overhead. It is used as a proxy for the inverse operator norm $\|r^{(k)}\|_{K^{-1}}$ appearing in (35).

## 3.2 Green preconditioning on uniform or structured computational grids

**General results** FFT-based methods draw a large part of their efficiency from Green (or Laplacian) preconditioning [27]. The principle of Green preconditioning is to use the solution to the problem (5) defined with a uniform reference medium to precondition the homogenization problem (2). The eigenvalues of the Green preconditioned operator have been analyzed by [51, 52, 53, 54] for scalar elliptic problems of the type $\text{div } c\nabla u$ where $c$ is a heterogeneous field of scalar material properties. For a FEM discretization and preconditioning by the inverse of $\text{div } c_0\nabla$, with $c_0$ a field of reference material properties, the preconditioned discrete operator eigenvalues are contained in intervals built
from the inf and sup of $c(x)/c_0(x)$ in the supports of the conforming finite element basis functions. The approach was extended to elliptic problems with vector primal variable $\mathbf{u}$ and tensor material properties $\mathbb{C}$ by [55, 56]. In particular, for a conformal FEM or finite difference discretization, bounds on the individual eigenvalues of the Green preconditioned operator can be constructed from the inf and sup of the eigenvalues of $C_0^{-1} : \mathbb{C}$ in the support function of the elements.

**Uniform grid** For a discretization on a uniform grid as in (14), i.e. with no grid adaptation, the Green preconditioner for $\mathbf{K}$ is $\mathbf{M} = \mathbf{K}_0$ defined in (28). The pseudo-inverse of $\mathbf{K}_0$ is the discrete Green operator $\mathbf{G}_0$. The later can be applied in the Fourier space from (29), (30), (32) and use of FFT. This is used e.g. in [24, 26, 27], and is equivalent to solving the strain based version $\mathbf{\Gamma}_0 * \mathbb{C} : \boldsymbol{\epsilon} = -\mathbf{\Gamma}_0 * \mathbb{C} : \mathbf{E}$, where $\mathbf{\Gamma}_0 = -\mathbf{B}\mathbf{K}_0\mathbf{B}^\mathrm{T}\mathbf{W}$, provided a Green operator norm of the residual is adopted [24, 25, 27].

**Transformed grid** In turn for the problem (21) discretized on a structured grid and pulled back on a uniform grid, the true Green preconditioner for $\tilde{\mathbf{K}}$ would be $\mathbf{M} = \tilde{\mathbf{K}}_0$, the stiffness matrix obtained from the finite element discretization on the transformed structured grid with uniform stiffness $C_0$. The latter can be constructed as $\tilde{\mathbf{K}}_0 = \tilde{\mathbf{B}}^\mathrm{T}\tilde{\mathbf{W}}\tilde{\mathbf{C}}_0\tilde{\mathbf{B}}$ where $\tilde{\mathbf{C}}_0$ is built from the pulled-back reference stiffness tensor $\tilde{\mathbf{C}}_0$ defined by replacing $\mathbb{C}$ by $C_0$ in (22). In that case, the bounds on the spectrum of eigenvalues of the preconditioned problem would be determined from local eigenvalues of $C_0^{-1} : \mathbb{C}$.

Unfortunately, the inverse of $\tilde{\mathbf{K}}_0$ is not known in general. Instead, the reference matrix $\mathbf{K}_0$ defined on the uniform grid is here used as a preconditioner of $\tilde{\mathbf{K}}$. The bounds on the spectrum of eigenvalues of the preconditioned problem are then determined from local eigenvalues of $C_0^{-1} : \tilde{\mathbb{C}}$, which are usually less favorable than those of $C_0^{-1} : \mathbb{C}$.

## 3.3 Bounds on the eigenvalue spectrum of the preconditioned operator

**Construction of bounds** Since FEM uses base functions with a local support, upper and lower bounds on the individual eigenvalues of the Green preconditioned discrete problem can be constructed a priori [55, 56]. The bounds on the eigenvalues of $\mathbf{G}_0\tilde{\mathbf{K}}$ (or of the generalized eigenvalue problem for $\tilde{\mathbf{K}}$ and $\mathbf{K}_0$) can be computed as follows:

1. For the reference node in pattern p of the structure, compute the largest $((\boldsymbol{\lambda}^+)_p)$ and lowest $((\boldsymbol{\lambda}^-)_p)$ of all eigenvalues of $C_0^{-1} : \tilde{\mathbb{C}}_{p-s}^q$ among all quadrature points in the support of that node, and store them in two vectors $\boldsymbol{\lambda}^\pm$ of size $N_p$.

2. Sort $\boldsymbol{\lambda}^\pm$ by increasing values, with repetition of each value d times.

The above bounds hold for elements with exact integration, in the shared row space of $\tilde{\mathbf{K}}$ and $\mathbf{K}_0$. There are $d$ zero eigenvalues of $\tilde{\mathbf{K}}$ and $\mathbf{K}_0$, and thus $\mathbf{G}_0\tilde{\mathbf{K}}$, corresponding to rigid body translations (rigid body rotations are excluded due to periodicity). For underintegrated elements, additional zero-energy hourglass modes arise [18, 57].

In the first step of the above procedure, for small strain linear elasticity, $\tilde{\mathbb{C}}$ and $C_0^{-1}$ can be expressed as $d^2 \times d^2$ matrices (Voigt notations with no minor symmetries). Clearly, the matrix for $C_0^{-1}$ has $d(d-1)/2$ zero eigenvalues, corresponding eigenvectors obtained from anti-symmetric second order tensors. So does the matrix of $C_0^{-1} : \tilde{\mathbb{C}}$: these zero eigenvalues must be discarded for the construction of the bounds.

**Eigenvalue clusters** In the particular case of piecewise constant material properties, eigenvalues of the FEM-discretized preconditioned operator can be localized. In particular, consider the case where there is a domain $\Omega_\alpha$ in which the stiffness $\mathbb{C}$ is uniform and a scalar multiple of the reference medium; i.e. $\mathbb{C} = \alpha C_0$. Then, each of the $N_\alpha$ basis functions with support in $\Omega_\alpha$ yield $d$ eigenvectors of $\mathbf{G}_0\mathbf{K}$ with eigenvalue $\alpha$ [56]. Hence $dN_\alpha$ eigenvalues cluster exactly at $\alpha$. The modal weight of
the right hand side defined by (34) for the corresponding modes is zero, because div C : E = 0 at these nodes. These modes correspond to type I described in [58] in the case of isotropic conductivity. This situation is extremely favorable for the convergence of the PCG, as only the effective spectrum of eigenvalues with non zero weight is involved in the resolution, see (35). Since eigenvectors are K and K₀-orthogonal, the set of modes with non-zero weight is included in the subspace generated by G₀f, where f are vectors of nodal forces supported in regions where C is not a scalar multiple of C₀, including nodes at interfaces. Such eigenvalue clusters have also been observed for a strain based method with discretization on trigonometric polynomials and numerical integration [8], and in the case of isotropic conductivity [58].

Further, when the stiffness C is bounded by α⁻C₀ ≤ C ≤ α⁺C₀ (in the sense of positive definite quadratic forms) in a domain Ω_α containing N_α basis functions, then the bounds constructed from the procedure in [53, 55, 56] guarantee the existence of a dN_α of eigenvalues within the interval [α⁻, α⁺], which is favorable when α⁻ and α⁺ are close.

**Effect of grid adaptation** Let's illustrate these properties of the spectrum on a simple example. Fig. 3 shows the effect of grid adaptation of the effective spectrum with Green preconditioning. This example features a circular inclusion with radius 3/8 of the periodic unit length. Both phases have uniform isotropic stiffness. The matrix stiffness C_m has Lamé parameters λ_m = 1, μ_m = 0.5. Two cases are considered for the inclusion stiffness C_i:

$$ (a) \quad C_i = 0.01 C_m \text{ (2nd row of Fig. 3)}; $$

$$ (b) \quad \lambda_i = 0.1\lambda_m, \mu_i = 0.01\mu_m \text{ (3rd row of Fig. 3)}. $$

In both cases, the stiffness of the reference medium for the Green preconditioner is chosen as C₀ = √Cᵢ : Cₘ. Thus for a discretization on a uniform grid and plane strain conditions, the eigenvalues of the preconditioned operator are bounded by 0.1 and 10 and the condition number is 100 in the two cases. The main difference is that in case (a), Cᵢ and Cₘ are scalar multiples of C₀, while they are not in case (b).

Three types of grid transformation are investigated, from left to right in Fig. 3:

(left) no transform. Elements in the uniform grid are attributed to the inclusion or the matrix based on the position of their center of mass;

(center) algebraic transform. A square of edge length 2R̃ = 5/8 L in the computational domain ˜Ω = [-L/2, L/2]ᵈ is transformed to a circle of radius R = 3/8 L in the physical domain Ω = [-L/2, L/2]ᵈ. Any point ˜x ∈ ˜Ω is radially transformed to x = ρ ˜x with ρ = 1 if ˜x = 0 else

$$ \rho = \begin{cases} 1 + \tilde{r}(R - R_i)/R_i^2 & \text{if } 0 < \tilde{y} \le \tilde{R} \\ [R + (R_m - R)(\tilde{r} - R_i)/(R_m - R_i)]/\tilde{r} & \text{else,} \end{cases} $$

$$ \text{where } \tilde{r} = \|\tilde{x}\|, \tilde{y} = \max_{i=1...d} |\tilde{x}_i|, R_m = tL/2 \text{ and } R_i = tR \text{ with } t = \tilde{r}/\tilde{y}; $$

(right) boundary node projection transform. Nodes of the uniform grid found at the boundary between elements attributed to matrix and inclusion in the no transform procedure are projected on the actual boundary of the inclusion. Other nodes are not transformed.

Bi-linear elements with 4 quadrature points are used to discretize the problem on a 32×32 grid. In this example with 2048 degrees of freedom, no more than one fourth of the modes (so here, 512) may have a non-zero weight, due to the symmetry of the problem, the transform, the discretization and the load. For the three types of grid transformation investigated, from left to right in Fig. 3:
<figure><img src="image_5.png" /><figcaption>Figure 3: Effect of grid adaptation on the spectrum and effective condition number $\kappa^{\text{eff}}$ of the Green preconditioned operator. Left: no transform, center: algebraic transform from a square to a circle, right: boundary node projection transform. Red and blue curves are lower and upper bounds [55, 56], dots are computed eigenvalues, colored by modal weight (34), gray vertical lines are weights.</figcaption></figure>
(left) *no transform*. The eigenvalue spectrum is bounded by 0.1 and 10. In case (a), 2nd row, eigenvalues cluster exactly at these two values since $C_i = 0.1C_0$ and $C_m = 10C_0$. The corresponding eigenvectors are respectively supported by nodes exactly in the inclusion or in the matrix, and the modal weights (34) are exactly zero (yellow). This situation, where the local stiffness is a scalar multiple of the reference stiffness, is thus similar to the isotropic conductivity situation studied by [58]. Further, numerical investigation of fields confirm that modes with non-zero weight belong to the subspace generated by $G_0f$, where $f$ is supported only by boundary nodes. Since there are 96 boundary nodes and 2 degrees of freedom per node, the maximum number of modes with non-zero weight is, accounting for symmetry, $192/4 = 48$, which is confirmed by Fig. 3a. The effective spectrum of these 48 modes with non-zero modal weight, corresponding to dots colored in green, has a condition number $\kappa^{\text{eff}} = \lambda_{\text{max}}^{\text{eff}}/\lambda_{\text{min}}^{\text{eff}} = 9.4$. This value is roughly 10 times smaller than that of the full spectrum, which is equal to $\kappa = 100$. In case (b), these two clusters diffuse and continuously spread in the intervals formed by the eigenvalues of $C_0^{-1} : C_i$ and $C_0^{-1} : C_m$, and some have non zero modal weight$^2$. This situation thus differs from the isotropic conductivity situation studied by [58]: since $C$ is not a scalar multiple of $C_0$, there is no type I mode. The number of eigenvalues with non-zero weight, 504, nearly attains the maximum limit of 512. The effective condition number raises to 91.9. Accordingly, the number of PCG iterations jumps from 27 in case (a) to 75 in case (b).

(center) *algebraic transform*. All elements are transformed. In both cases (a) and (b), eigenvalues are now continuously distributed, mostly in two diffuse clusters spread around 0.1 and 10. The number of modes with non-zero modal weight nearly attains the maximum limit of 512. The bounds on the largest and smallest non-zero eigenvalues are tight in that case. The effective condition number is larger than 900. The number of PCG iterations increases to 152 in case (a) and 202 in case (b).

(right) *boundary node projection transform*. Most elements are not transformed. In both cases, the spectrum preserves some features of the spectrum of the operator on a uniform grid. In particular for case (a), nodes with support on undeformed elements contribute to eigenvalues exactly 0.1 or 10 with zero modal weight. There are 288 nodes belonging to transformed elements, times 2 degrees of freedom per node, hence, accounting for symmetry, at most 288 × 2/4 = 144 modes with non-zero weight. This limit is attained in the numerical experiment. Further, the bounds on the largest and smallest non-zero eigenvalues are loose in that case, which means that the excessive localized transformation of a few elements is not as problematic as what could be anticipated. As a result, the effective condition number is lower than for the algebraic transform. PCG converges in 70 iterations for case (a) and 95 for case (b).

An important difference between the two types of transforms is observed for the bounds on the smallest and largest non-zero eigenvalues: they provide a good estimate of the effective condition number for the algebraic transform, but overestimate it for the boundary node projection one. To our current knowledge, these bounds cannot be improved. Yet a rationale is proposed to explain the differences obtained for the two types of transforms. Modes have energy contributions from all elements in the support of shape functions of activated nodes. In the *boundary node projection transform*, as one node is projected on the boundary, some of its surrounding elements are stretched and the others shrunk. This has adverse effects on the pulled-back stiffness of elements in the support of its shape function, and thus on the eigenvalue. On the contrary, in the algebraic transform, regions spanning over several elements are either stretched or shrunk, and effects add up. Modes that activate nodes with shape functions entirely supported in these regions may thus generate extreme eigenvalues. The same situation can be encountered for other types of grid transformation where node density is locally modified, such as in [30, 32].

<footnote n="2">²In Fig. 3, dots with non-zero weight are displayed above nodes with zero weight. There is actually an important portion of modes with exactly zero modal weight in all situations investigated in that figure, throughout the spectrum.</footnote>
## 4 Implementation

**Iterative solver** The implementation relies on the standard PCG algorithm [42, 44] with an operator based, matrix free implementation. Main operations are: computation of the right hand side, application of the stiffness operator **K** or **K̃**, and application of the Green preconditioner **G**₀. The first two require application of the discrete differential operators and of the local constitutive model. The last one requires FFT and application of the Green function in the Fourier space. The memory requirement of PCG is four displacement-like vectors, including the solution and the r.h.s.

**Differential operators** The discrete gradient (9) and divergence (13) are short kernel convolutions. These operations are performed in the space domain using stencils built from the $B_g^q$ defined in (9). The quadrature weights are fused in the divergence operation [27, sec. 5.1].

<table><caption>**Algorithm 1** $\tilde{f} \leftarrow \tilde{K}\tilde{u}$ – Application of stiffness operator $\tilde{K}$ (21) on a displacement-like vector $\tilde{u}$.</caption><tr><td>1:</td><td>$\tilde{f} \leftarrow 0$</td><td></td></tr><tr><td>2:</td><td>for all $q \in 1...N_q$ do</td><td>$\triangleright q^{\text{th}}$-quadrature point</td></tr><tr><td>3:</td><td>$F^q \leftarrow I + (\tilde{\nabla}\tilde{u}_{\Phi})^q$</td><td>$\triangleright q^{\text{th}}$-transformation gradient</td></tr><tr><td>4:</td><td>$\epsilon^q \leftarrow \text{sym}((\tilde{\nabla}\tilde{u})^q \cdot (F^q)^{-1})$</td><td>$\triangleright q^{\text{th}}$-transformed gradient</td></tr><tr><td>5:</td><td>$\sigma^q \leftarrow \mathbb{C}^q : \epsilon^q$</td><td>$\triangleright$ constitutive model</td></tr><tr><td>6:</td><td>$\tilde{f} \leftarrow \tilde{f} + w_q \text{div}(\det(F^q)\sigma^q \cdot (F^q)^{-T})$</td><td>$\triangleright$ transformed divergence</td></tr><tr><td>7:</td><td>end for</td><td></td></tr></table>

**Transformed stiffness operator** The transformed stiffness operator $\tilde{K}$ defined in (21) applies to a displacement-like vector $\tilde{u}$ as sketched in Alg. 1. As compared to an implementation on a uniform grid (see e.g. [24, Alg. 2]), the only modification to account for grid transformation is the presence of the transformation gradient **F**, computed in step Alg. 1.3. For reduced integration with a single quadrature point ($N_q = 1$), the inverse of the transformation gradient can be cached to save one gradient computation and inversion at each iteration. For $N_q > 1$, storing $N_q$ copies of $F^{-1}$ is in general prohibitive<fnref n="3" />, so it is recomputed at each iteration. The constitutive model is applied locally at the current quadrature point $q$ of each pattern in step Alg. 1.5. Vectors $\tilde{\nabla}^q \tilde{u}$, $\epsilon$ and $\sigma$ can share the same memory. Hence two buffer $\nabla u$-like vectors are required to apply $\tilde{K}$.

**Right hand side** The r.h.s. is computed as in Alg. 1, except step Alg. 1.4 which is replaced by $\epsilon \leftarrow E$. Once the system is solved, the strain in the physical domain is obtained at quadrature point $q$ as $\epsilon^q \leftarrow \text{sym}((\tilde{\nabla}\tilde{u})^q \cdot (F^q)^{-1}) + E$.

**Green preconditioner** The FFT operator and its inverse can operate in-place to apply the Green function $\hat{G}_0$ in Fourier space. When $N_q = 1$, the operation $\hat{G}_{0k} \cdot \hat{f}_k$ is directly applied using (33). Otherwise, at each iteration and Fourier index, the matrix $\hat{K}_{0k}$ is assembled on-the-fly and the small $d \times d$ system $\hat{f}_k = \hat{K}_{0k}\hat{u}_k$ solved. Due to the symmetry of $\hat{G}_0$, a viable alternative could be to cache the real and symmetric matrices $\hat{G}_{0k}$ for $k \in \{0...N_1 \div 2, ..., 0...N_d \div 2\}$. The memory overhead would be only $3/4 \times N_p$ floating-point scalars. This optimization is not done in the current implementation, as the on-the-fly procedure did not reveal performance critical.

<footnote n="3">³For the boundary node projection transform, $\tilde{u}_{\Phi}$ and its gradient are sparse, so $F^{-1}$ could be cached (not done in the current implementation).</footnote>
<figure><img src="image_6.png" /><figcaption>Figure 4: Test microstructures (a) diamond inclusions forming a checkerboard along dashed lines ($d = 2$) (b) Hashin's composite inclusion ($d = 2$ or 3) (c) Boolean model of spheres ($d = 3$).</figcaption></figure>

**Memory requirement** To sum up, the memory requirement is 5 displacement-like vectors (including grid nodal displacements $\tilde{u}_{\Phi}$), i.e. $5dN_p$ floating-point scalars. In addition, at least one integer scalar vector is required to encode microstructure phases. As a rule of the thumb, for a careful implementation with 64-bit floating point numbers, simulations on grids of more than $400^3$ nodes can fit on a desktop computer with 16 GB RAM.

**Computer programming** The algorithms have been implemented using the Julia programming language [59] and libraries from the Julia ecosystem [60, 61, 62] with the FFTW backend [63]. All simulations have been run on a 10-year old laptop with 16 GB RAM and a 2.80/3.80 GHz processor with 6 MB cache. The source code for a minimal working example, consisting in a non-optimized yet compact Python implementation for 2D problems, is provided in supplementary materials [64].

# 5 Numerical investigations

Numerical experiments are carried out to investigate the influence of combined Green preconditioning and grid adaptation on the accuracy of homogenized properties, the accuracy of local fields and the number of iterations to convergence of PCG. We first investigate two test microstructures for which exact solutions are available, and then a more complex boolean sphere model.

## 5.1 Diamond inclusions

The first test microstructure is shown in Fig. 4(a). Two diamond inclusions, actually squares rotated by a 45 degree angle w.r.t. the periodic unit cell axes, are embedded in a matrix. The diamonds are placed to form a periodic checkerboard along dashed lines. The boundary between material phases is non-smooth due to right angles at the tips of diamonds.

The homogenization problem is formulated in terms of 2 dimensional thermal conductivity with scalar conductivities (see Appendix A), to benefit from the analytical solution provided in [65]. The inclusion conductivity is set to 0.01, the matrix one to 1, while the choice of the reference conductivity has no influence. The analytical solution of the temperature gradient and flux is not defined at the inclusion boundaries and on the dashed lines in Fig. 4(a). Hence, in practice the diamonds are shifted by half a pixel to the right w.r.t. to the configuration of Fig. 4(a) when $Q_1$ elements are used. This precaution is necessary in order to shift the centers of the pixels of the
<figure><img src="image_7.png" /><figcaption>Figure 5: Diamond inclusions in 2D conduction with macroscopic temperature gradient $E = e_x$.</figcaption></figure>

uniform grid away from these lines where numerical issues can arise while evaluating the analytical solution of [65].

Simulations are carried out with either bilinear rectangular $Q_1$ elements in the reference grid, with either 1 or 4 quadrature points, or linear triangular $P_1$ elements (see Fig. 1). Since inclusions are rotated, rectangular elements on a uniform grid are not boundary-conforming. Similarly, triangular elements following the pattern of Fig. 1(c) on a uniform grid conform to only half of the inclusion edges. The grid adaptation strategy used here corresponds to the *boundary node projection transform* presented in section 3.3. Simulations are carried out on grids ranging from $20^2$ to $2560^2$ nodes to study h-convergence, with unit imposed macroscopic temperature gradient $E = e_x$. Additional simulations with bilinear rectangular $Q_1$ elements on a regular grid are also carried out using the composite pixel technique [38]. The Voigt mixing rule is used, since it has been shown to be the most appropriate when the conductivity of the inclusion is much smaller than that of the matrix [38].

Fig. 5 shows fields of temperature gradient $\epsilon_x$ and error $\delta\epsilon_x = \epsilon_{x\text{num}} - \epsilon_{x\text{ref}}$ on $\epsilon_x$, as compared to the exact solution, for the grid with $40^2$ elements. In the post-processing for $Q_1$ elements, the gradient of $u$ is evaluated at the (transformed position) of the center of each element, for both the simulations with 1 or 4 quadrature points. The overall error on the local fields is quantified in Fig. 6(b) by the root mean square error (RMSE) over all element centers
<figure><img src="image_8.png" /><figcaption>Figure 6: Diamond inclusions in 2D conductivity. (a) relative error on homogenized conductivity (b) root mean square error (RMSE) on temperature gradient field $\epsilon$ for $E = e_x$.</figcaption></figure>

$$ \sqrt{(N_p N_e)^{-1} \sum_{P \in \mathcal{I}(N)} \sum_{e=1 \dots N_e} J_P \delta \epsilon_P^e \cdot \delta \epsilon_P^e} $$

The field on the uniform grid with rectangular $Q_1$ elements and 1 quadrature point presents significant hourglassing, in particular close to diamond tips. Hourglassing is removed by using 4 integration points, which ensures exact integration for the uniform grid, but significant errors occur in elements at the boundary, and to a lesser extent inside the low-conductivity inclusion. Boundary-conforming grid adaptation significantly reduces the error on the fields, as well as hourglassing. The strain field for triangular $P_1$ elements is free of hourglassing, but the error on the $\epsilon_x$ field presents a regular pattern. This pattern is associated to a mismatch between the exact solution, evaluated at element centers, and the derivatives of the numerical solution $\mathbf{u}$ in the $x$ direction, which are more accurate at the middle of horizontal edges on a regular grid. As opposed to simulations on boundary-conforming grids, simulations on uniform grids with composite pixels do not significantly improve the temperature gradient field in this case. For all strategies, largest errors are located at the inclusion tips, where the exact solution presents a singularity.

The error on the homogenized conductivity scales between $h^{1.1}$ and $h^{1.5}$, see Fig. 6(a). The $h^2$ scaling expected for boundary-conforming $P_1$ or $Q_1$ elements is not achieved, because the inclusion boundary is not smooth enough and the exact solution has a singularity at the inclusion tips. Excessive transformation of elements (quadrangles transformed into triangles are no longer shape-regular) could be also blamed; yet results for the two dimensional composite circular inclusion in Appendix B are well behaved despite excessive element transformation. For quadrangular $Q_1$ elements with 1 quadrature point, grid adaptation reduces this error by roughly a factor 4 w.r.t. the uniform grid, and roughly 30-40 (!) for 4 quadrature points. For triangular $P_1$ elements, grid adaptation does not improve the error on $C^{\text{hom}}$ for the diamond inclusion case, while a significant improvement is observed for the composite circular inclusion in Appendix B. The RMSE on the temperature gradient field $\epsilon$ has a sublinear $h$-scaling, and is best for transformed grids. Fig. 7(a) shows that the number of iterations for PCG to converge plateaus for grid sizes above 100$^2$. Grid adaptation doubles to triples the number of iterations. In turn, the use of composite pixels more than doubles the number of iterations, and here does not improve the discretization errors.
<figure><img src="image_9.png" /><figcaption>Figure 7: PCG iteration count (a) Diamond inclusions in 2D conduction. (b) Hashin's composite inclusion in 3D linear elasticity.</figcaption></figure>

## 5.2 Hashin's composite inclusion

The second test microstructure has now smooth boundaries: it is Hashin's composite inclusion [66, 67, 68], see Fig. 4(b). The inclusion core ($i = 1$), shell ($i = 2$) and the matrix ($i = 3$) have uniform isotropic stiffness $C_i$ with Lamé parameters $\lambda_i, \mu_i$ and bulk modulus $\kappa_i = \lambda_i + 2\mu_i/d$. When an isotropic load $E = I$ is applied$^4$, the composite inclusion is neutral provided the matrix bulk modulus $\kappa_3$ is set to the equivalent bulk modulus of the inclusion:

$$ \kappa^{\text{hom}} = \kappa_2 \left( 1 - \frac{d\alpha\phi}{1 + \alpha\phi} \right) \quad \text{where} \quad \phi = (R_1/R_2)^d \quad ; \quad \alpha = \frac{d(\kappa_2 - \kappa_1)}{(d-1)2\mu_2 + d\kappa_1} \quad (36) $$

In that case, $\epsilon = I$ in the matrix, is uniform in the core, and given in each phase $i$ by

$$ [(a_i + b_i/r^d)I - db_i/r^d e_r \otimes e_r] \cdot E \quad \text{where} \quad r = \|x\| \quad \text{and} \quad e_r = x/r, \quad (37) $$

with $a_2 = 1/(1 + \alpha\phi)$, $a_1 = (1 + \alpha)a_2$, $b_2 = \alpha R_1^d a_2$ and $b_1 = 0$. As a result, provided $\kappa_3$ is set to $\kappa^{\text{hom}}$, the solution to Hashin's composite inclusion is exact both on an infinite domain with uniform remote strain, or on a finite domain with periodic boundary conditions.

In the simulations, the inclusion's shell Lamé parameters are set to $\lambda_2 = 1$, $\mu_2 = 0.5$ and the core ones to $\lambda_1 = 0.001$, $\mu_1 = 0.005$. The reference stiffness is $C_0 = \sqrt{C_1 : C_2}$. The matrix bulk modulus $\kappa_3$ is set to $\kappa^{\text{hom}}$ (36). Its shear modulus $\mu_3$ is set such that $C_3$ is a scalar multiple of $C_0$. Figs. 8 and 9 present results for three dimensional linear elasticity. Discretization grids have from $20^3$ to $320^3$ nodes and hexahedral elements with either one, two or eight integration points. In [69, Appendix A], simulations with two superposed integration points for hexahedral elements have been shown to coincide with the finite difference scheme with double tetrahedral stencil introduced by [70]. For these three types of integration, additional simulations on regular grids with the composite voxel technique are provided for comparison, using the Voigt mixing rule. In all cases, the strain and stress are evaluated only at the center of the elements during the post-processing.

<footnote>$^4$Plane strain and plane isotropic load is assumed for the cylindrical inclusion in dimension $d = 2$</footnote>
<figure><img src="image_10.png" /><figcaption>Figure 8: Hashin's composite inclusion in 3D linear elasticity with $\bar{\epsilon} = I$ (a) $\epsilon_{xx}$ (b) norm of absolute error on $\epsilon$. Tri-linear $Q_1$ elements with 1 (top) or 8 (bottom) quadrature points, without (left) or with (right) grid adaptation.</figcaption></figure>

Fig. 8 presents strain fields $\epsilon_{xx}$ on the lowest resolution grid ($20^3$). Without grid adaption, large errors appear (locally exceeding 100% of $\|E\|$), in particular for elements at the boundary which is poorly approximated by cubes. Again, the boundary node projection transform (see section 3.3) is applied. Errors on the strain field are drastically reduced after transformation. Using eight quadrature points, the error drops below $\|\delta\epsilon\| < 10\%\|E\|$ in all elements, while it can exceed 20% in elements at the highly contrasted core-shell boundary while using only 1 quadrature point. However, boundary-conforming grid adaptation failed for hexahedral elements with two integration points, because the boundary-conforming grid adaptation considered here modifies the topology of elements. Indeed, hexahedral elements with three faces and seven nodes belonging to the boundary are transformed into tetrahedrons (see e.g. [32, Figs. 11,13] and Fig. 8), up to the curvature of the boundary. As a result, out of the two tetrahedrons defined in the stencils by [70], one has its four nodes on the boundary. The Jacobian of the transformation for this stencil is zero, so that the pulled-back stiffness defined in (22) is undefined.

Results of Fig. 8 are confirmed upon grid refinement in Fig. 9(b). The RMSE on the strain field, computed as $\sqrt{N_p^{-1} \sum_{p \in \mathcal{I}(N)} J_p \delta\epsilon_p : \delta\epsilon_p}$, scales as $h^{1/2}$ on uniform grids, with or without composite voxels, and $h^{3/2}$ on boundary-conforming transformed grids. It roughly doubles if only one quadrature point is used instead of eight. For uniform grids, using two or eight integration points, the RMSE on the strain field coincides and the strain field is free of hourglassing. Meanwhile, the relative error on the homogenized bulk modulus (Fig. 9(a)) reaches a quadratic $h$-convergence on boundary-conforming grids, while it is only sublinear on uniform grids, with or without composite voxels. Since the inclusion boundaries are smooth enough, the exact weak solution is at least in $H^2$ in each phase. Spatial convergence estimates [35] indeed ensure a $h^2$ rate of convergence of the homogenized stiffness on the boundary-conforming grid, instead of a $h^1$ rate on the uniform grid. Similarly, for the local strain field, the spatial convergence rate is expected to at least $h^1$ on the boundary-conforming grid, instead of $h^{1/2}$ on the uniform grid. Similar results are obtained in plane strain elasticity or two dimensional conductivity, see Figs. 12, 13, 14, 15 in Appendix B.
<figure><img src="image_11.png" /><figcaption>Figure 9: Hashin's composite inclusion in 3D linear elasticity. (a) relative error on homogenized bulk modulus (b) absolute error on strain field $\epsilon$ for $\bar{\epsilon} = I$.</figcaption></figure>

The number of iterations for PCG convergence stabilizes under grid refinement for uniform grids, with or without composite voxels, (Fig. 7(b)), as expected from the properties of the spectrum of the Green preconditioned operator, see section 3.3. The boundary node projection transform grid adaptation more than doubles the number of iterations; yet no stabilization is observed under grid refinement. The increase of iterations is however largely balanced by the considerable gain in accuracy.

## 5.3 Boolean model of spheres

The last investigated microstructure is a Boolean model of spheres, see Fig. 4. The domain $\Omega$ is partitioned in two subsets: inside (phase 1) or outside (phase 2) the spheres. The periodic unit cell comprises 1000 spheres of radius $R = 0.08L$. The spheres are thrown randomly and may overlap. Both phases 1 and 2 are path connected domains. Simulations are carried on grids with $25^3$ to $300^3$ hexahedral elements. To avoid the presence of inconsistent meshes, the distance between (any periodic reproduction of) two sphere centers cannot lie in the open interval $[2R, 2R + \delta]$ [where $\delta = \sqrt{3}L/50$. This implies that spheres either overlap, either are separated by at least one voxel at low resolution ($\ge 50^3$)].

In the simulations, Lamé parameters of the outside spheres phase are set to $\lambda_2 = 1$, $\mu_2 = 0.5$. The stiffness of the inside spheres phase is $C_1 = 0.01C_2$, and the reference one $C_0 = \sqrt{C_1 : C_2}$. The applied macroscopic strain is $E = e_x \otimes e_x$.

Grid adaptation is done as follows. Each element of the uniform grid is first assigned to the phase inside or outside spheres based on its center position. The signed distance of nodes to the spheres is computed, and the identification number of the sphere corresponding to this smallest distance saved in a nodal vector. Then, nodes at the boundary between the two phases are projected to the nearest sphere, according to the previously saved sphere identification number. No specific treatment of the intersection of two or more spheres is done, so that the elements are not conforming in these regions, see Fig. 4. After that procedure, few elements at complex sphere intersection have extremely low
<figure><img src="image_12.png" /><figcaption>Figure 10: Strain field $\epsilon_{xx}$ for $E = e_x \otimes e_x$ for the 3D boolean model of spheres (xy plane cut). Top row: grid $50^3$, middle row: grid $250^3$, bottom row: grid $250^3$ (zoomed in).</figcaption></figure>

Jacobian. A simple mitigation rule is used: if an element has a Jacobian below the tolerance 0.01 in any quadrature point, all corresponding nodal displacements are reduced by 10%. This typically concerns a few tens of elements out of several millions, and has a minor impact on the homogenized stiffness.

The effect of grid adaptation on the quality of the strain field is illustrated in Fig. 10 for two grid sizes, $50^3$ and $250^3$, with one or eight quadrature points. At both resolutions, a close inspection reveals that reduced integration produces hourglassing, which is stronger on the uniform grid than on the transformed grid. Regardless of resolution and number of quadrature points, fields on uniform grids oscillate in elements near boundaries. This is cured by the boundary-conforming grid adaptation.

Fig. 11(a) shows that the number of iterations to PCG convergence is stable with grid size for uniform grids and increasing for transformed ones. Grid transformation can increase the number of iterations by up to a factor 4 or 8 for 8 or 1 quadrature points respectively. In turn, Fig. 11(b) shows that the homogenized stiffness $C^{\text{hom}}$ has a qualitatively better $h$-convergence on transformed grids (near quadratic) than on uniform ones (linear). Using 1 quadrature point, the same error is attained for a transformed grid with $h/R \approx 0.31$ ($40^3$ nodes) or a uniform grid with $h/R = 0.05$
<figure><img src="image_13.png" /><figcaption>Figure 11: Boolean model of spheres in 3D linear elasticity. (a) PCG iteration count (b) homogenized stiffness.</figcaption></figure>

(250³ nodes).

## 5.4 Discussion

Numerical experiments in this section focus on a simple boundary-conforming grid adaptation strategy, the *boundary node projection transform*. This low-cost strategy proved sufficient to reach $h^2$ spatial convergence of the homogenized properties and *h* convergence of the RMSE on the strain field for problems with sufficiently smooth solution. This remarkably departs from the *h* spatial convergence of the homogenized properties observed for different types of FFT-based methods on uniform grids, see e.g. [10, 19, 38, 39]. In turn, this simple grid adaptation preserves favorable features of the spectrum of eigenvalues and induces an acceptable increase in the number of iterations of PCG.

Clearly, the present framework naturally applies to more involved types of grid transformation. In future works, the present grid adaptation strategy can be improved in the light of recent developments. The grid adaptation technique newly introduced in [32] is promising as it provides more control on the transformation of elements, including the possibility to increase the node density near the boundary. Since the obtained transformation is also boundary-conforming, $h^2$ spatial convergence of the homogenized properties is expected. The spring stiffness parameters controlling the grid transformation may provide better compromises on the spatial convergence and number of iterations of the solver than the simple *boundary node projection transform*.

In turn, the grid adaptation strategy based on optimal transport theory introduced in [30] allows to control node density in regions where the solution exhibits high strain gradients. However, numerical experiments in [30] indicate that while the spatial convergence is improved with this strategy, the spatial convergence of the homogenized properties remains of order *h*. Indeed, the transformed grid is not boundary-conforming. Thus, a two-step grid adaptation may prove more beneficial to reach $h^2$ convergence: first transform using optimal transport theory, then conform to boundary. This requires also to shift from a discretization on a trigonometric basis to an FE one.

On the down side, local node densification allowed by transforms considered in [32, 30] can induce
modes with extreme eigenvalues in the regions with highest or lowest mesh density, and thus increase
the number of iterations of the solver. Further, the computation of these transforms requires solving
non-linear problems on the grid, which introduces a non-negligible overhead.

# 6 Conclusion

## Main findings

1. "FFT-based methods" with FEM discretization on structured grid can be efficiently formulated with isoparametric grid transformation.

2. The Green preconditioning is applied in the reference, or pulled-back, grid.

3. The grid adaptation strategy has a critical impact on the spectrum of the Green preconditioned operator and on convergence of PCG. Bounds on the individual eigenvalues of the spectrum have been provided. Grid transformations for which most elements are kept untransformed can preserve favorable features of the spectrum such as clustering. Such grid transformation can however increase PCG iterations by a factor two to eight.

4. Boundary-conforming (or nearly) grids drastically improve the *h*-convergence of the homogenized property $C^{\text{hom}}$ w.r.t. uniform grids with or without composite voxels. For problems with sufficiently smooth boundary, a quadratic *h*-convergence of $C^{\text{hom}}$ is achieved for bi- or tri-linear $Q_1$ elements.

5. Boundary-conforming transforms significantly reduce errors on local strain fields. Hourglassing of bi or tri-linear elements with reduced integration can be diminished.

## Perspectives

1. Boundary-conforming grids benefit from *h*-convergence bounds, hence the use of higher order (quadratic) elements could provide interesting computational cost vs. accuracy compromises. Extension to isogeometric analysis is a promising prospect [71].

2. Efficient grid adaptation remains a bottleneck. The focus has here been on moving only nodes near boundaries, which is simple yet not optimal. Alternative strategies borrowing from e.g. elliptic grid generation techniques [40, 32], or optimal transport theory [30] with subsequent projection of nodes on boundaries, remain to be investigated. The ideal transform should seek optimisation of the accuracy while maintaining a favorable eigenvalue spectrum of the Green preconditioned operator.

# Conflicts of interest

The authors declare no competing financial interest.

# Dedication

The manuscript was written through contributions of all authors. All authors have given approval
to the final version of the manuscript.
Acknowledgments

FB is supported by Nantes Université for a semester as a visiting professor at Czech Technical Uni-
versity in Prague. FB acknowledges funding by l'Agence Nationale de la Recherche (ANR), project
ANR-23-CE08-0036 (MOCAMOR). ML acknowledges funding by the European Commission (Marie
Skłodowska-Curie Fellowship 101106585 — microFFTTO), and Cluster of Excellence livMatS. IP
and JZ acknowledge funding by the European Union under the project ROBOPROX (reg. no.
CZ.02.01.01/00/22_008/0004590).

Open Access

For the purpose of Open Access, a CC-BY public copyright licence has been applied by the authors
to the present document and will be applied to all subsequent versions up to the Author Accepted
Manuscript arising from this submission.

Supplementary data

Supporting information for this article, including example scripts, is available in [64].

A Conductivity

For thermal conductivity, the homogenization problem (2) is modified as follows. The scalar $u$
denotes the periodic fluctuation temperature and the vector $\boldsymbol{\sigma}$ the heat flux. The temperature
gradient $\boldsymbol{\epsilon}$ is split in $\boldsymbol{\epsilon} = \boldsymbol{E} + \nabla u$ where $\boldsymbol{E}$ is the macroscopic temperature gradient. The local
conductivity is the second order tensor $C$, and $\boldsymbol{\sigma} = -\boldsymbol{C} \cdot \boldsymbol{\epsilon}$. For grid transformation, the pulled-back
conductivity is $\tilde{C}_{ij} = J F_{ik}^{-1} C_{k\ell} F_{\ell j}^{-T}$ in place of (22).

For Hashin's composite inclusion problem in section 5.2 and Appendix B, the three phases have isotropic conductivity $C_i = C_i I$. The equivalent conductivity of the inclusion is given by (36) up to the replacement of $\kappa$ by $C$ and modification $\alpha = (C_2 - C_1)/((d-1)C_2 + C_1)$. The temperature gradient is still obtained by (37), where $\boldsymbol{E}$ is now a vector.

B Hashin’s composite inclusion in dimension 2

References

[1] H. Moulinec and Pierre Suquet. A fast numerical method for computing the linear and non-linear properties of composites. *Comptes Rendus de l'Académie des Sciences*, 2(318):1417-1423, 1994.

[2] H. Moulinec and Pierre Suquet. A numerical method for computing the overall response of non-linear composites with complex microstructure. *Computational Methods in Applied Mechanical Engineering*, 157:69-94, 1998.

[3] Matti Schneider. A review of nonlinear fft-based computational homogenization methods. *Acta Mechanica*, 232(6):2051-2010, 2021.

[4] S Lucarini, M V Upadhyay, and J Segurado. Fft based approaches in micromechanics: fundamentals, methods and applications. *Modelling and Simulation in Materials Science and Engineering*, 30(2):023002, dec 2021.
<figure><img src="image_14.png" /><figcaption>Figure 12: Hashin's composite inclusion in 2D conductivity with $E = e_x$, $C_1 = 0.01$ and $C_2 = 1$. Only one quarter or one half of fields is shown by symmetry.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 13: Hashin's composite inclusion in 2D conductivity for $E = e_x$, $C_1 = 0.01$ and $C_2 = 1$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 14: Hashin's composite inclusion in 2D elasticity with $E = I_{2d}$. Only one quarter or one half of fields is shown by symmetry.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 15: Hashin's composite inclusion in 2D linear elasticity for $E = I_{2d}$.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 16: PCG iteration count for Hashin's composite inclusion in 2D for (a) conduction and (b) linear elasticity problems.</figcaption></figure>
[5] Jan Zeman, Jaroslav Vondřejc, Jan Novák, and Ivo Marek. Accelerating a fft-based solver for numerical homogenization of periodic media by conjugate gradients. *Journal of Computational Physics*, 229(21):8065–8071, 2010.

[6] Sébastien Brisard and Luc Dormieux. Fft-based methods for the mechanics of composites: A general variational framework. *Computational Materials Science*, 49:663–671, 2010.

[7] H. Moulinec and F. Silva. Comparison of three accelerated fft-based schemes for computing the mechanical response of composite materials. *International Journal for Numerical Methods in Engineering*, 97(13):960–985, 2014.

[8] N. Mishra, J. Vondřejc, and J. Zeman. A comparative study on low-memory iterative solvers for FFT-based homogenization of periodic media. *Journal of Computational Physics*, 321:151–168, 2016.

[9] Matti Schneider. Lippmann-schwinger solvers for the computational homogenization of materials with pores. *International Journal for Numerical Methods in Engineering*, 121(22):5017–5041, 2020.

[10] Sébastien Brisard and Luc Dormieux. Combining galerkin approximation techniques and the principle of hashin and shtrikman to improve two fft-based numerical methods for the homogenization of composites. *Computer Methods in Applied Mechanics and Engineering*, 217-220:197–212, 2012.

[11] F. Bignonnet and L. Dormieux. Fft-based bounds on the permeability of complex microstructures. *International Journal for Numerical and Analytical Methods in Geomechanics*, 38:1707–1723, 2014.

[12] J. Vondřejc, J. Zeman, and I. Marek. An FFT-based Galerkin method for homogenization of periodic media. *Computers & Mathematics with Applications*, 68(3):157–173, 2014.

[13] J. Vondřejc, J. Zeman, and I. Marek. Guaranteed upper-lower bounds on homogenized properties by FFT-based Galerkin method. *Computer Methods in Applied Mechanics and Engineering*, 297:258–291, 2015.

[14] J. Zeman, T. W. J. de Geus, J. Vondřejc, R. H. J. Peerlings, and M. G. D. Geers. A finite element perspective on nonlinear fft-based micromechanical simulations. *International Journal for Numerical Methods in Engineering*, 111(10):903–926, 2017.

[15] T.W.J. de Geus, J. Vondřejc, J. Zeman, R.H.J. Peerlings, and M.G.D. Geers. Finite strain fft-based non-linear solvers made simple. *Computer Methods in Applied Mechanics and Engineering*, 318:412–430, 2017.

[16] F. Bignonnet. Efficient fft-based upscaling of the permeability of porous media discretized on uniform grids with estimation of rve size. *Computer Methods in Applied Mechanics and Engineering*, 369:113237, 2020.

[17] François Willot, Bassam Abdallah, and Yves-Patrick Pellegrini. 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*, 98(7):518–533, 2014.

[18] François Willot. Fourier-based schemes for computing the mechanical response of composites with accurate local fields. *Comptes Rendus Mécanique*, 343(3):232–245, 2015.
[19] Matti Schneider, Felix Ospald, and Matthias Kabel. Computational homogenization of elasticity on a staggered grid. *International Journal for Numerical Methods in Engineering*, 105(9):693-720, 2016.

[20] Stéphane Berbenni, Vincent Taupin, Komlan Sénam Djaka, and Claude Fressengeas. A numerical spectral approach for solving elasto-static field dislocation and g-disclination mechanics. *International Journal of Solids and Structures*, 51(23):4157-4175, 2014.

[21] Ricardo A. Lebensohn and Alan Needleman. Numerical implementation of non-local polycrystal plasticity using fast fourier transforms. *Journal of the Mechanics and Physics of Solids*, 97:333-351, 2016. SI:Pierre Suquet Symposium.

[22] Y. Chen, L. Gélébart, A. Marano, and J. Marrow. Fft phase-field model combined with cohesive composite voxels for fracture of composite materials with interfaces. *Comput Mech*, 68:433-457, 2021.

[23] 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*, 109(4):459-486, 2017.

[24] Matti Schneider, Dennis Merkert, and Matthias Kabel. Fft-based homogenization for microstructures discretized by linear hexahedral elements. *International Journal for Numerical Methods in Engineering*, 109(10):1461-1489, 2017.

[25] Matthias Leuschner and Felix Fritzen. Fourier-accelerated nodal solvers (fans) for homogenization problems. *Computational Mechanics*, 62(3):359-392, Sep 2018.

[26] R. J. Leute, M. Ladecký, A. Falsafi, I. Jödicke, I. Pultarová, J. Zeman, T. Junge, and L. Pastewka. Elimination of ringing artifacts by finite-element projection in FFT-based homogenization. *Journal of Computational Physics*, 453:110931, 2022.

[27] Martin Ladecký, Richard J. Leute, Ali Falsafi, Ivana Pultarová, Lars Pastewka, Till Junge, and Jan Zeman. An optimal preconditioned fft-accelerated finite element solver for homogenization. *Applied Mathematics and Computation*, 446:127835, 2023.

[28] Flavia Gehrig and Matti Schneider. An x-fft solver for two-dimensional thermal homogenization problems. *International Journal for Numerical Methods in Engineering*, 126(7):e70022, 2025.

[29] Miroslav Zecevic, Ricardo A. Lebensohn, and Laurent Capolungo. New large-strain fft-based formulation and its application to model strain localization in nano-metallic laminates and other strongly anisotropic crystalline materials. *Mechanics of Materials*, 166:104208, 2022.

[30] Cédric Bellis and Renaud Ferrier. Numerical homogenization by an adaptive fourier spectral method on non-uniform grids using optimal transport. *Computer Methods in Applied Mechanics and Engineering*, 419:116658, 2024.

[31] Rodrigo Santos-Güemes, Gonzalo Álvarez, and Javier Segurado. An fft based adaptive grid framework to represent non-singular dislocations. *Mechanics of Materials*, 194:105004, 2024.

[32] Miroslav Zecevic, Ricardo A. Lebensohn, and Laurent Capolungo. Achieving geometric accuracy in fft-based micromechanical models using conformal grid. *Mechanics of Materials*, 212:105512, 2026.

[33] Vladimir D. Liseikin. *Grid Generation Methods*. Springer International Publishing, Cham, 2017.
[34] Susanne C. Brenner and L. Ridgway Scott. *The Mathematical Theory of Finite Element Methods*. Springer New York, New York, NY, 2008.

[35] Ivo Babuška. The finite element method for elliptic equations with discontinuous coefficients. *Computing*, 5(3):207–213, Sep 1970.

[36] Ivo Babuška, Uday Banerjee, and Kenan Kergrene. Strongly stable generalized finite element method: Application to interface problems. *Computer Methods in Applied Mechanics and Engineering*, 327:58–92, 2017. Advances in Computational Mechanics and Scientific Computation—the Cutting Edge.

[37] S. A. Voitsekhovskii and V. L. Makarov. A rate of convergence bound of regular-grid difference schemes for elliptical equations with discontinuous coefficients. *Journal of Soviet Mathematics*, 58(1):17–21, Jan 1992.

[38] Matthias Kabel, Dennis Merkert, and Matti Schneider. Use of composite voxels in fft-based homogenization. *Computer Methods in Applied Mechanics and Engineering*, 294:168–188, 2015.

[39] Jaroslav Vondřejc. Improved guaranteed computable bounds on homogenized properties of periodic media by the fourier–galerkin method with exact integration. *International Journal for Numerical Methods in Engineering*, 107(13):1106–1135, 2016.

[40] Joe F. Thompson, Bharat K. Soni, and Nigel P. Weatherill. *Handbook of Grid Generation*. CRC Press, 1998.

[41] Alexandre Ern and Jean-Luc Guermond. *Theory and Practice of Finite Elements*. Springer New York, NY, 2004.

[42] Yousef Saad. *Iterative Methods for Sparse Linear Systems*. Society for Industrial and Applied Mathematics, second edition, 2003.

[43] Jörg Liesen and Zdenek Strakos. *Krylov Subspace Methods: Principles and Analysis*. Oxford University Press, 10 2012.

[44] Henk A. van der Vorst. *Iterative Krylov Methods for Large Linear Systems*. Cambridge Monographs on Applied and Computational Mathematics. Cambridge University Press, 2003.

[45] Matthias Kabel, Thomas Böhlke, and Matti Schneider. Efficient fixed point and newton–krylov solvers for fft-based homogenization of elasticity at large deformations. *Computational Mechanics*, 56(6):1497–1514, 2014.

[46] Matti Schneider. On the barzilai-borwein basic scheme in fft-based computational homogenization. *International Journal for Numerical Methods in Engineering*, 118(8):482–494, 2019.

[47] Dietrich Braess. *Finite Elements: Theory, Fast Solvers, and Applications in Solid Mechanics*. Cambridge University Press, 2007.

[48] A. Jennings. Influence of the eigenvalue spectrum on the convergence rate of the conjugate gradient method. *IMA Journal of Applied Mathematics*, 20(1):61–72, 08 1977.

[49] Owe Axelsson and Gunhild Lindskog. On the rate of convergence of the preconditioned conjugate gradient method. *Numerische Mathematik*, 48:499–523, 1986.

[50] A. van der Sluis and H.A. van der Vorst. The rate of convergence of conjugate gradients. *Numerische Mathematik*, 48:543–560, 1986.
[51] Bjørn Fredrik Nielsen, Aslak Tveito, and Wolfgang Hackbusch. Preconditioning by inverting the laplacian: an analysis of the eigenvalues. *IMA Journal of Numerical Analysis*, 29(1):24–42, 2008.

[52] Tomáš Gergelits, Kent-André Mardal, Bjørn Fredrik Nielsen, and Zdeněk Strakoš. Laplacian preconditioning of elliptic pdes: Localization of the eigenvalues of the discretized operator. *SIAM Journal on Numerical Analysis*, 57(3):1369–1394, 2019.

[53] Tomáš Gergelits, Bjørn Fredrik Nielsen, and Zdeněk Strakoš. Generalized spectrum of second order differential operators. *SIAM Journal on Numerical Analysis*, 58(4):2193–2211, 2020.

[54] Tomáš Gergelits, Bjørn Fredrik Nielsen, and Zdeněk Strakoš. Numerical approximation of the spectrum of self-adjoint operators in operator preconditioning. *Numerical Algorithms*, 91:301–325, 2022.

[55] Ivana Pultarová and Martin Ladecký. Two-sided guaranteed bounds to individual eigenvalues of preconditioned finite element and finite difference problems. *Numerical Linear Algebra with Applications*, 28(5):e2382, 2021.

[56] Martin Ladecký, Ivana Pultarová, and Jan Zeman. Guaranteed two-sided bounds on all eigenvalues of preconditioned diffusion and elasticity problems solved by the finite element method. *Applications of Mathematics*, 66(1):21–42, 2021.

[57] Matti Schneider. Voxel-based finite elements with hourglass control in fast fourier transform-based computational homogenization. *International Journal for Numerical Methods in Engineering*, 123(24):6286–6313, 2022.

[58] C. Bellis and H. Moulinec. Lippmann-schwinger spectrum, composite materials eigenstates and their role in computational homogenization. *International Journal for Numerical Methods in Engineering*, 126(19):e70130, 2025.

[59] Jeff Bezanson, Alan Edelman, Stefan Karpinski, and Viral B Shah. Julia: A fresh approach to numerical computing. *SIAM Review*, 59(1):65–98, 2017.

[60] Juan Ignacio Polanco. "pencilarrays.jl: Distributed julia arrays using the mpi protocol", 2021.

[61] Juan Ignacio Polanco. Pencilffts.jl: Ffts of mpi-distributed julia arrays, 2021.

[62] Kristoffer Carlsson and Fredrik Ekre. Tensors.jl — tensor computations in julia. *Journal of Open Research Software*, 2019.

[63] Matteo Frigo and Steven G. Johnson. The design and implementation of FFTW3. *Proceedings of the IEEE*, 93(2):216–231, 2005. Special issue on "Program Generation, Optimization, and Platform Adaptation".

[64] François Bignonnet, Martin Ladecký, Ivana Pultarová, and Jan Zeman. Supplementary material to "fourier-based computational micromechanics on boundary-conforming transformed grids", February 2026.

[65] R. V. Craster and Yu. V. Obnosov. Four-phase checkerboard composites. *SIAM Journal on Applied Mathematics*, 61(6):1839–1856, 2001.

[66] Zvi Hashin. The elastic moduli of heterogeneous materials. *Journal of Applied Mechanics*, 29(1):143–150, 03 1962.

[67] Zvi Hashin and B. Walter Rosen. The elastic moduli of fiber-reinforced materials. *Journal of Applied Mechanics*, 31(2):223–232, 06 1964.
[68] Zvi Hashin and B. Walter Rosen. Erratum: “the elastic moduli of fiber-reinforced materials” (journal of applied mechanics, 1964, 31, pp. 223–232). *Journal of Applied Mechanics*, 32(1):219–219, 03 1965.

[69] Flavia Gehrig and Matti Schneider. Element-based internal variable formulations for finite element discretizations in fft-based homogenization methods. *International Journal for Numerical Methods in Engineering*, 126(21):e70170, 2025.

[70] A. Finel. A tetrahedron-based discretization for fft-based computational homogenization with smooth solution fields. *Computer Methods in Applied Mechanics and Engineering*, 436:117703, 2025.

[71] P. Antolin and T. Hirschler. Quadrature-free immersed isogeometric analysis. *Engineering with Computers*, 38(5):4475–4499, Oct 2022.