%~Mouliné par MaN_auto v.0.40.4 (2aa04cce) 2026-03-19 13:41:41
\documentclass[CRMECA,Unicode,biblatex]{cedram}

\addbibresource{CRMECA_Bignonnet_20250849.bib}

\usepackage{subcaption}
\usepackage{upgreek}
\usepackage{algorithm}
\usepackage{algpseudocode}
\usepackage{array}

\DeclareMathOperator{\sym}{sym}

\let\grad\nabla
\DeclareMathOperator{\gradsym}{\grad^{\sym}}

\let\divisionsymbol\div
\let\div\relax
\DeclareMathOperator{\div}{div}

\renewcommand{\hom}{\mathrm{hom}}
\newcommand{\eff}{\mathrm{eff}}
\newcommand{\num}{\mathrm{num}}
\newcommand{\rmref}{\mathrm{ref}}

\newcommand{\uuuu}[1]{\mathbb{#1}}
\newcommand{\vv}[1]{\mathbf{\boldsymbol{#1}}}
\newcommand{\uu}[1]{\mathbf{\boldsymbol{#1}}}
\newcommand{\uv}[1]{\mathbf{\boldsymbol{#1}}}

\newcommand{\ov}[1]{\overline{#1}}

\newcommand\mathbftt[1]{\textnormal{\ttfamily\bfseries #1}}
\newcommand\mathbfsf[1]{\textnormal{\sffamily\bfseries #1}}
\newcommand\vecdof[1]{\mathbfsf{#1}}
\newcommand\matdof[1]{\mathbfsf{#1}}
\newcommand\tuple[1]{\mathbftt{#1}}

\newcommand{\diff}{\mathop{}\!{\operatorfont{d}}}

\makeatletter
\renewcommand\p@subfigure{}
\makeatother

%%%---------------------------------------------------------------------------

%%% DELIMITERS

\usepackage{mathtools}
\DeclarePairedDelimiter{\parens}{\lparen}{\rparen}
\DeclarePairedDelimiter{\braces}{\{}{\}}
\DeclarePairedDelimiter{\abs}{\lvert}{\rvert}
\DeclarePairedDelimiter{\norm}{\lVert}{\rVert}
\DeclarePairedDelimiter{\bracks}{[}{]}

%%%---------------------------------------------------------------------------

%%% VERTICAL BAR

%%% ''Such that'' macro
\newcommand*{\st}[1][]{\;#1\vert\;}

%%%---------------------------------------------------------------------------

\graphicspath{{./figures/}}

\newcommand*{\mk}{\mkern -1mu}
\newcommand*{\Mk}{\mkern -2mu}
\newcommand*{\mK}{\mkern 1mu}
\newcommand*{\MK}{\mkern 2mu}

\hypersetup{urlcolor=purple, linkcolor=blue, citecolor=red}

\newcommand*{\relabel}{\renewcommand{\labelenumi}{(\theenumi)}}
\newcommand*{\romanenumi}{\renewcommand*{\theenumi}{\roman{enumi}}\relabel}
\newcommand*{\Romanenumi}{\renewcommand*{\theenumi}{\Roman{enumi}}\relabel}
\newcommand*{\alphenumi}{\renewcommand*{\theenumi}{\alph{enumi}}\relabel}
\newcommand*{\Alphenumi}{\renewcommand*{\theenumi}{\Alph{enumi}}\relabel}
\let\oldtilde\tilde
\renewcommand*{\tilde}[1]{\mathchoice{\widetilde{#1}}{\widetilde{#1}}{\oldtilde{#1}}{\oldtilde{#1}}}
\let\oldhat\hat
\renewcommand*{\hat}[1]{\mathchoice{\widehat{#1}}{\widehat{#1}}{\oldhat{#1}}{\oldhat{#1}}}
\let\oldforall\forall
\renewcommand*{\forall}{\mathrel{\oldforall}}

\title{Fourier-based computational micromechanics on boundary-conforming transformed grids}
\alttitle{Homogénéisation par FFT sur grilles transformées conformes aux interfaces}

\author{\firstname{François} \lastname{Bignonnet} \CDRorcid{0000-0002-5279-9112}\IsCorresp}
\address{Nantes Université, École Centrale Nantes, CNRS, GeM, UMR 6183, F-44600 Saint-Nazaire, France}
\email{francois.bignonnet@univ-nantes.fr}

\author{\firstname{Martin} \lastname{Ladecký}\CDRorcid{0000-0002-1337-7603}}
\address{Department of Microsystems Engineering, University of Freiburg, Georges-K\"ohler-Allee 103, 79110 Freiburg, Germany}
\email{martin.ladecky@imtek.uni-freiburg.de}

\author{\firstname{Ivana} \lastname{Pultarov{\'a}}\CDRorcid{0000-0003-0460-9270}}
\address{Faculty of Civil Engineering, Czech Technical University in Prague, Thákurova 7, 166~29 Prague 6, Czech Republic}
\email{ivana.pultarova@cvut.cz}

\author{\firstname{Jan} \lastname{Zeman}\CDRorcid{0000-0003-2503-8120}}
\address[3]{Faculty of Civil Engineering, Czech Technical University in Prague, Thákurova 7, 166~29 Prague 6, Czech Republic}
\email{jan.zeman@cvut.cz}

\thanks{Agence Nationale de la Recherche (ANR), project ANR-23-CE08-0036 (MOCAMOR). European Commission (Marie Skłodowska-Curie Fellowship 101106585 -- microFFTTO). Cluster of Excellence livMatS. European Union project ROBOPROX (reg.~no.~CZ.02.01.01/00/22\_008/0004590). The first author thanks Nantes Université for supporting a semester as a visiting professor at Czech Technical University in Prague.}
\CDRGrant[ANR]{ANR-23-CE08-0036}
\CDRGrant[European Commission]{101106585}
\CDRGrant[European Union]{CZ.02.01.01/00/22\_008/0004590}

\keywords{\kwd{Computational homogenization} \kwd{fast Fourier transform} \kwd{Green function preconditioning} \kwd{finite element discretization} \kwd{grid adaptation}}
\altkeywords{\kwd{Homogénéisation numérique} \kwd{transformée de Fourier rapide} \kwd{préconditionnement par opérateur de Green} \kwd{discrétisation par éléments finis} \kwd{adaptation de grille}}

\begin{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.
\end{abstract}

\begin{altabstract}
Les méthodes d'homogénéisation accélérées par un préconditionnement par l'opérateur de Green avec transformée de Fourier rapide (FFT) sont usuellement mises en œuvre sur une grille uniforme, ce qui limite la qualité de la discrétisation. Nous considérons ici une représentation plus précise de la géométrie, obtenue en transformant la grille uniforme en une grille conforme aux interfaces. Le problème mécanique est discrétisé par la méthode des éléments finis (FEM) avec une transformation isoparamétrique des éléments. L'adaptation aux interfaces peut nécessiter localement une importante transformation géométrique de la grille, ce qui est naturellement pris en compte par la discrétisation par éléments finis. Des bornes rigoureuses sont fournies sur le spectre des valeurs propres du système discret préconditionné par l'opérateur de Green. Pour des transformations de grille où seuls les nœuds les plus près des interfaces y sont projetés, les valeurs propres modifiées correspondent à des vecteurs propres localisés aux interfaces entre constituants, de sorte que le spectre effectif reste favorable pour le solveur par gradient conjugué préconditionné. Des investigations numériques confirment que la précision des propriétés homogénéisées et des champs locaux obtenus sur une grille conforme aux interfaces est grandement améliorée par rapport à celle sur une grille uniforme, pour une augmentation modérée du temps de calcul.
\end{altabstract}

\begin{document}
%\input{CR-pagedemetas}
%\end{document}
\maketitle

\section{Introduction}\label{sec:introduction}

\subsection{State of the art}\label{sec:state-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~\cite{Moulinec1994,Moulinec1998}, have become popular (see reviews~\cite{Schneider2021a,Lucarini2022}). 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 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~\cite{Zeman2010,Brisard2010,Moulinec2014,Mishra2016,Schneider2020}. In the meantime, these methods have been recast in a Galerkin approximation framework, either based on the Hashin--Shtrikman variational formulation~\cite{Brisard2010,Brisard2012,Bignonnet2014} or on the classical weak formulation which allows to reuse all the FEM technology at no compromise on efficiency~\cite{Vondrejc2014,Vondrejc2015,Zeman2017,deGeus2017}.

One of the drawbacks of the initial method of Moulinec \& Suquet~\cite{Moulinec1994,Moulinec1998} 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~\cite{Brisard2010,Brisard2012,Bignonnet2020}. 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~\cite{Vondrejc2014,Zeman2017,deGeus2017}. Another way to remedy oscillations is the reinterpretation in the Fourier space of discrete differential operators used in finite difference methods~\cite{Willot2014,Willot2015}. A suppression of oscillations by this strategy is not systematically achieved, but is possible at the expense of the introduction of staggered grid~\cite{Schneider2015}, or of a numerical diffusion of the fields due to the inconsistent treatment of first and second order derivatives~\cite{Berbenni2014,Lebensohn2016,Chen2021}.

More recently, a finite element method (FEM) discretization on a regular grid opened a new way to remove ringing artifacts~\cite{Brisard2017,Schneider2017,Leuschner2018}. This strategy has been extended to voxel meshing by a periodic pattern or higher order finite elements by Leute et al.~\cite{Leute2022}. Ladeck\'y et al.~\cite{Ladecky2023} 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~\cite{Gehrig2025}. We investigate a different idea here, elaborating on the framework of grid transformation investigated by~\cite{Zecevic2022,Bellis2024,SantosGuemes2024,Zecevic2026}.

Zecevic et al.~\cite{Zecevic2022} and Bellis \& Ferrier~\cite{Bellis2024} independently suggested to transform the uniform grid into an arbitrary periodic structured grid. For large strain computations, Zecevic et al.~\cite{Zecevic2022} 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~\cite{Bellis2024} 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~\cite{Bellis2024} 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.~\cite{SantosGuemes2024} 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.~\cite{Zecevic2026} introduced a method for the generation of boundary-conforming grids. The transformation of the grid is obtained by solving a discrete system of differential equations~\cite{Liseikin2017}. 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.

\subsection{Contributions}\label{sec:contributions}

The present work focuses on the following objectives:
\begin{enumerate}
\item introduce a transparent framework for grid adaptation using finite element discretization with isoparametric transformation of elements;
\item 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;
\item 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.
\end{enumerate}
We show how finite element discretization on regular grids~\cite{Schneider2017, Leuschner2018, Leute2022, Ladecky2023} 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~\cite{Brenner2008, Babuska1970}. 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) \cite{Babuska2017}. 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$~\cite{Babuska1970,Voitsekhovskii1992}. 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.~\cite{Brisard2012,Schneider2015,Kabel2015,Vondrejc2016}.

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~\cite{Zecevic2022,Bellis2024}, which can generate unfavorable eigenvalue spectrum when used together with Green preconditioning. Alternative and promising grid adaptation strategies~\cite{Thompson1998,Zecevic2026} 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 shape. 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~\cite{Zecevic2026}. The discretization scheme by~\cite{Willot2015}, equivalent to hexahedral finite elements with reduced integration~\cite{Schneider2017}, is used in~\cite{Zecevic2026} to capture this feature. On the contrary, the types of grid transformation investigated in~\cite{Bellis2024,SantosGuemes2024} are not \emph{exactly} boundary-conforming: although large element transformation can arise due to node densification near boundaries, the shape of voxels is not modified.

The article is organized as follows. Section~\ref{sec:fe-discretization-homogenization-problem} recalls the principles of the FEM discretization of the homogenization problem on uniform grids, and introduces its straightforward extension to transformed, structured grids. Section~\ref{sec:green-preconditioning-eigenvalues} 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~\ref{sec:implementation} details the implementation of the method. Section~\ref{sec:numerical-investigations} 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.

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

\section{Finite element discretization of homogenization problems on structured grids}\label{sec:fe-discretization-homogenization-problem}

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 \dots \times [0, L_d]$ where $L_i > 0$ are the side lengths of $\Omega$. Let $\mathcal{V} = \braces[\big]{\uv{u} \in H^1(\Omega) \st[\big] \text{$\Omega$-periodic}, \ \ov{\uv{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 $\uv{f}$ over $\Omega$ is denoted:
\begin{equation}\label{eq:volume-average}
\ov{\uv{f}} = \dfrac{1}{\abs{\Omega}} \int_\Omega \uv{f} \diff \uv{x}.
\end{equation}
The set $L^2_{\sym}(\Omega)$ denotes the space of symmetric second-order tensor fields, which are $\Omega$-periodic and square integrable over $\Omega$.

\subsection{Homogenization problem}\label{sec:homogenization-problem}

The homogenization problem for small strain linear elasticity consists in finding a microscopic fluctuation displacement field $\uv{u} \in \mathcal{V}$ and a Cauchy stress field $\uu{\sigma} \in L^2_{\sym}(\Omega)$ such that
\begin{equation}\label{eq:homogenization-problem}
\begin{aligned}
	& \div \uu{\sigma} = 0 & \text{in $\Omega$},
\\	& \uu{\sigma} = \uuuu{C} : \uu{\epsilon} & \text{in $\Omega$},
\\	& \uu{\epsilon} = \uu{E} + \gradsym \uv{u} & \text{in $\Omega$},
\end{aligned}
\end{equation}
where $\uu{\epsilon} \in L^2_{\sym}(\Omega)$ is the linearized strain, and $\uu{E}$ the macroscopic strain. The stiffness tensor field~$\uuuu{C}$ is heterogeneous. The homogenization problem~\eqref{eq:homogenization-problem} can be reformulated as: find $\uv{u} \in \mathcal{V}$ such that
\begin{equation}\label{eq:stiffness-operator}
\uu{K} \uv{u} = \uv{f}_{\uu{E}}
\quad \text{where}
\quad \uu{K} = - \div \uuuu{C} : \gradsym
\quad \text{and}
\quad \uv{f}_{\uu{E}} = \div \uuuu{C} : \uu{E}.
\end{equation}
The weak formulation of problem~\eqref{eq:stiffness-operator} is
\begin{equation}\label{eq:homogenization-problem-weak}
\forall \uv{v} \in \mathcal{V}, \quad \int_\Omega \gradsym \uv{v} : \uuuu{C} : (\uu{E} + \gradsym \uv{u}) \diff \uv{x} = 0.
\end{equation}
Problem~\eqref{eq:stiffness-operator} will be compared to a reference problem: find $\uv{u} \in \mathcal{V}$ such that
\begin{equation}\label{eq:reference-stiffness-operator}
\uu{K}_0 \uv{u} = \uv{f}
\quad \text{where}
\quad \uu{K}_0 = - \div \uuuu{C}_0 : \grad^{\sym}.
\end{equation}
In the above reference problem, $\uuuu{C}_0$ is a uniform reference stiffness and $\uv{f}$ an arbitrary body force field satisfying $\ov{\uv{f}} = 0$. The solution to the reference problem~\eqref{eq:reference-stiffness-operator} is formally expressed as a convolution with the Green function $\uu{G}_0$
\begin{equation}\label{eq:green-function}
\uv{u} = \uu{G}_0 \ast \uv{f}
\quad \text{where}
\quad \uu{G}_0 = \uu{K}_0^{-1}.
\end{equation}
After suitable discretization, the Green function can be represented by a matrix and will be used as a preconditioner to solve the homogenization problem~\eqref{eq:homogenization-problem}.

\subsection{Finite element discretization on uniform grids}\label{sec:fe-discretization-uniform}

The finite element discretization is first presented on uniform grids, i.e.\ without grid adaptation. The presentation is adapted from~\cite{Schneider2017,Leute2022,Ladecky2023}, 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~\ref{sec:fe-discretization-structured}.

To avoid proliferation of notations, the presentation is restricted to grids made up of periodic patterns comprising either:
\begin{itemize}
\item a single first order tensor product Lagrange finite element $Q_1$~\cite{Ern2004}, i.e.\ a quadrangular ($d=2$) or hexahedral ($d=3$) element;
\item several first order simplicial Lagrange finite elements $P_1$~\cite{Ern2004}, i.e.\ triangular ($d=2$) or tetrahedral ($d=3$) elements.
\end{itemize}
For a generalization to periodic patterns made up of several elements of order $k$, the reader is referred to~\cite{Ladecky2023}.

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

\begin{figure}[!h]
\begin{subfigure}{.3\textwidth}
\includegraphics[width=\textwidth]{fem/fem_pattern_square_4_quad}
\caption{1 quadrangular $Q_1$ element with 4 quadrature points ($N_e=1$, $N_q=4$)}\label{fig:fem-discretization_a}
\end{subfigure}
\hfill
\begin{subfigure}{.3\textwidth}
\includegraphics[width=\textwidth]{fem/fem_pattern_square_1_quad}
\caption{1 quadrangular $Q_1$ element with 1 quadrature point ($N_e=1$, $N_q=1$)}\label{fig:fem-discretization_b}
\end{subfigure}
\hfill
\begin{subfigure}{.3\textwidth}
\includegraphics[width=\textwidth]{fem/fem_pattern_triangle_2_quad}
\caption{2 triangular $P_1$ elements and 2 quadrature points ($N_e=2$, $N_q=2$)}\label{fig:fem-discretization_c}
\end{subfigure}
\caption{Finite element discretization: examples of pattern $\mathcal{P}_{\tuple{p}}$ made up of one quadrangular $Q_1$ element with \eqref{fig:fem-discretization_a}~$N_q = 4$ quadrature points, \eqref{fig:fem-discretization_b}~$N_q = 1$ quadrature point, or \eqref{fig:fem-discretization_c} two triangular $P_1$ elements and $N_q = 2$ quadrature points.}\label{fig:fem-discretization}
\end{figure}

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 Figure~\ref{fig:fem-discretization}. Nodes of the reference pattern are separated in two types. There is a fundamental node at position $\uv{y}^n = \uv{0} \in \mathcal{P}$. There also are periodic reproductions of the fundamental node $\uv{y}^n$ by translations $\tuple{s} \uv{h}$ where shifts $\tuple{s} \in \tuple{S} = \{0, 1\}^{d}$ and $\uv{h}$ is the stride of the pattern (equivalent to pixel/voxel edge lengths). There are $N_q$ quadrature points at positions $\uv{y}^q \in \mathcal{P}$ where $q \in \{1, \dotsc, N_q\}$ is the quadrature point index in the pattern.

The periodic reproduction in pattern $\mathcal{P}_{\tuple{p}}$ of the reference node are $\uv{x}^n_{\tuple{p}} = \uv{y}^n + \tuple{p} \uv{h}$ and similarly for quadrature points. There are $N_p$ nodes, $N_pN_e$ elements and $N_p N_q$ quadrature points in the grid. Nodal values of the displacement are denoted $\vecdof{u}_{\tuple{p}} = \uv{u}(\uv{x}^n_{\tuple{p}})$ and can be assembled in a vector of degrees of freedom $\vecdof{u} =\{\vecdof{u}_{\tuple{p}}\}_{\tuple{p} \in \mathcal{I}(\tuple{N})}$ containing $d N_p$ scalars. Similarly, stress at quadrature points are denoted $\uu{\upsigma}^q_{\tuple{p}} = \uu{\sigma}(\uv{x}^q_{\tuple{p}})$ and assembled in a vector $\uu{\upsigma} = \{\uu{\upsigma}^q_{\tuple{p}}\}_{\tuple{p} \in \mathcal{I}(\tuple{N}), q \in \{1, \dotsc, N_q\}}$ containing $d_\ast N_p N_q$ scalars where $d_\ast = d(d+1)/2$.

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

\paragraph{Interpolation}
At a point $\uv{y}$ of the reference element $\mathcal{E}^e$, displacements are interpolated from nodal values by interpolation functions $\phi^e_{\tuple{s}}$ such that $\phi^e_{\tuple{s}} (\uv{y}^n + \tuple{r} \uv{h}) = 1$ if $\tuple{r} = \tuple{s}$ else 0 for all shifts $\tuple{r}, \tuple{s} \in \tuple{S}$. At a point $\uv{x}$ of an element $\mathcal{E}^e_{\tuple{p}}$, the displacements are interpolated as:
\begin{equation}\label{eq:fe-interpolation}
\uv{u}(\uv{x}) = \sum_{\tuple{s} \in \tuple{S}} \vecdof{u}_{\tuple{p}+\tuple{s}} \phi^e_{\tuple{s}}(\uv{y})
\quad \text{where}
\quad \uv{y} = \uv{x} - \tuple{p} \uv{h} \in \mathcal{E}^e,
\end{equation}
where $\mathcal{E}^e$ is a reference element. The gradient of the displacement is defined in the interior of the element $\mathcal{E}^e_{\tuple{p}}$ as:
\begin{equation}\label{eq:fe-interpolation-gradient}
\grad \uv{u}(\uv{x}) = \sum_{\tuple{s} \in \tuple{S}} \vecdof{u}_{\tuple{p}+\tuple{s}} \otimes \grad \phi^e_{\tuple{s}}(\uv{y}).
\end{equation}

\paragraph{Discrete gradient}
The gradient of displacement at quadrature points $\uv{x}_{\tuple{p}}^q$ in $\mathcal{E}^e$ is then
\begin{equation}\label{eq:fe-gradient-quadrature}
\grad \vecdof{u}^q_{\tuple{p}} = \sum_{\tuple{s} \in \tuple{S}} \vecdof{u}_{\tuple{p}+\tuple{s}} \otimes \matdof{B}^{q}_{\tuple{s}}
\quad \text{where}
\quad \matdof{B}^{q}_{\tuple{s}} = \grad \phi^e_{\tuple{s}}(\uv{y}^q).
\end{equation}
All gradients at quadrature points are gathered in a vector $\grad \vecdof{u} = \{\grad \vecdof{u}^q_{\tuple{p}}\}_{\tuple{p} \in \mathcal{I}(\tuple{N}), q \in \{1, \dotsc, N_q\}}$ of size $d^2 N_p N_q$. The discrete gradient operator $\tilde{\matdof{B}}$ of size $d^2 N_p N_q \times d N_p$, obtained by assembly of the $\matdof{B}^{q}_{\tuple{s}}$, is defined such that $\grad \vecdof{u} = \tilde{\matdof{B}} \vecdof{u}$. Similarly, the operator $\matdof{B}$ of size $d_\ast N_p N_q \times d N_p $ of the symmetric gradient is such that $\gradsym \vecdof{u} = \matdof{B} \vecdof{u}$.

\paragraph{Discrete divergence}
The divergence is the negative value of the adjoint of the gradient operator, by the application of the Gauss theorem:
\begin{equation}\label{eq:gauss-theorem}
\forall \uv{u} \in \mathcal{V}, \quad \int_\Omega \grad \uv{u} : \uu{\sigma} \diff \uv{x} = - \int_\Omega \uv{u} \div \uu{\sigma} \diff \uv{x}.
\end{equation}
The discrete approximation by Gauss quadrature of the l.h.s.\ of~\eqref{eq:gauss-theorem} is for any test function $\uv{u}$ interpolated from test nodal values $\vecdof{u}$:
\begin{equation}\label{eq:gauss-quadrature}
\int_\Omega \grad \uv{u} :\uu{\sigma} \diff \uv{x} \approx \sum_{\tuple{p} \in \mathcal{I}(\tuple{N})} \sum_{q=1}^{N_q} w_q \grad \vecdof{u}^q_{\tuple{p}}: \uu{\upsigma}_{\tuple{p}}^q,
\end{equation}
where $w_q$ are quadrature weights and $\grad \vecdof{u}^{q}_{\tuple{p}}$ is defined in~\eqref{eq:fe-gradient-quadrature}. Then $\grad \vecdof{u}^q_{\tuple{p}}: \uu{\upsigma}_{\tuple{p}}^q = \vecdof{u}_{\tuple{p}+\tuple{s}} \cdot \uu{\upsigma}_{\tuple{p}}^q \cdot \matdof{B}^{q}_{\tuple{s}}$ by symmetry of $\uu{\upsigma}_{\tuple{p}}^q$. Reorganization of the sum on $\tuple{p}$ by the shifted Cartesian index $\tuple{p}' = \tuple{p} +\tuple{s}$ yields
\begin{equation}\label{eq:divergence-adjoint-of-gradient}
\sum_{\tuple{p}'\in \mathcal{I}(\tuple{N})} \sum_{q=1}^{N_q} w_q \sum_{\tuple{s} \in \tuple{S}} \vecdof{u}_{\tuple{p}'} \cdot \uu{\upsigma}_{\tuple{p}'-\tuple{s}}^q \cdot \matdof{B}^{q}_{\tuple{s}} \coloneqq - \sum_{\tuple{p}' \in \mathcal{I}(\tuple{N})} \vecdof{u}_{\tuple{p}'} \cdot \div \uu{\upsigma}_{\tuple{p}'},
\end{equation}
where one defines, by identification with the r.h.s.\ of~\eqref{eq:gauss-theorem}:
\begin{equation}\label{eq:fe-divergence-nodes}
\div \uu{\upsigma}_{\tuple{p}} \coloneqq \sum_{q=1}^{N_q} w_q \div \uu{\upsigma}^{q}_{\tuple{p}}
\quad \text{where}
\quad \div \uu{\upsigma}^{q}_{\tuple{p}} \coloneqq - \sum_{\tuple{s} \in \tuple{S}} \uu{\upsigma}^q_{\tuple{p}-\tuple{s}} \cdot \matdof{B}^{q}_{\tuple{s}}.
\end{equation}
The contribution of neighboring quadrature points $q$ to the divergence $\div \uu{\upsigma}_{\tuple{p}}$ at the fundamental node of pattern $\tuple{p}$ is $\div \uu{\upsigma}^{q}_{\tuple{p}}$. It is the counterpart of $\grad \vecdof{u}^{q}_{\tuple{p}}$ defined in~\eqref{eq:fe-gradient-quadrature}. Finally, gathering in a vector $\div \uu{\upsigma} = \{\div \uu{\upsigma}_{\tuple{p}}\}_{\tuple{p} \in \mathcal{I}(\tuple{N})}$ of size $d N_p$, the discrete divergence operator (for symmetric second order tensors) stemming from the finite element discretization is $- \matdof{B}^{\mathsf{T}} \matdof{W}$ such that $\div \uu{\upsigma} = - \matdof{B}^{\mathsf{T}} \matdof{W} \uu{\upsigma}$, where $\matdof{W}$ is a diagonal $d_\ast N_p N_q \times d_\ast N_p N_q$ matrix built from the quadrature weights $w_q$.

\paragraph{Discrete problem}
The weak formulation~\eqref{eq:homogenization-problem-weak} of the homogenization problem~\eqref{eq:homogenization-problem}, combined with the finite element discretization and Gauss quadrature approximation~\eqref{eq:gauss-quadrature} leads to the discrete version of~\eqref{eq:stiffness-operator}:
\begin{equation}\label{eq:fe-discrete-problem}
\matdof{K} \vecdof{u} = \vecdof{f}_{\uu{E}}
\quad \text{where}
\quad \matdof{K} = \matdof{B}^{\mathsf{T}}\matdof{W} \matdof{C} \matdof{B} 
\quad \text{and}
\quad \vecdof{f}_{\uu{E}} = - \matdof{B}^{\mathsf{T}}\matdof{W} \matdof{C} \vecdof{E},
\end{equation}
where $\matdof{C}$ is the $d_\ast N_p N_q \times d_\ast N_p N_q$ block-diagonal matrix containing the stiffness at quadrature points, $\matdof{K}$ is the $d N_p \times d N_p$ stiffness matrix, $\vecdof{E}$ is the $d_\ast N_p N_q$ vector built from the uniform the macroscopic strain $\uu{E}$ and $\vecdof{f}_{\uu{E}}$ the vector of generalized forces.

\subsection{Finite element discretization on structured grids}\label{sec:fe-discretization-structured}

\paragraph{Grid transformation}
The geometrical transformation from the uniform grid to the structured grid is performed by the mapping $\uv{x} = \uv{\Phi}(\tilde{\uv{x}})$ from the computational domain $\tilde{\Omega}$ to the physical domain $\Omega$, see Figure~\ref{fig:grid-adaptation}. The mapping displacement is $\widetilde{\uv{u}_{\Phi}}(\tilde{\uv{x}}) = \uv{\Phi}(\tilde{\uv{x}}) - \tilde{\uv{x}}$. This displacement belongs to $H^1(\Omega)$ and is $\Omega$-periodic. The fluctuation displacement $\uv{u}(\uv{x})$ in the physical domain is mapped from $\tilde{\uv{u}}(\tilde{\uv{x}}) = \uv{u}\parens[\big]{\uv{\Phi}(\tilde{\uv{x}})}$ in $\tilde{\Omega}$. The displacement $\tilde{\uv{u}} \in \tilde{\mathcal{V}}$, defined as $\mathcal{V}$ with the modification $\int_{\tilde{\Omega}} \tilde{\uv{u}} \diff \tilde{\uv{x}} = 0$. The gradient of the displacement $\grad \uv{u} = \partial \uv{u} / \partial \uv{x}$ in the physical domain $\Omega$ is related to the gradient $\tilde{\grad} \tilde{\uv{u}}= \partial \tilde{\uv{u}} / \partial \tilde{
\uv{x}}$ in the computational domain $\tilde{\Omega}$ by chain rule
\begin{equation}\label{eq:gradient-transformed}
\tilde{\grad} \tilde{\uv{u}} = \dfrac{\partial \tilde{\uv{u}}}{\partial \tilde{\uv{x}}} = \grad \uv{u} \cdot \uu{F}
\quad \text{where}
\quad \uu{F} = \tilde{\grad} \uv{\Phi}; \quad J = \det(\uu{F}).
\end{equation}
The weak formulation~\eqref{eq:homogenization-problem-weak} of the homogenization problem in the physical domain $\Omega$ is pulled back to the reference domain $\tilde{\Omega}$:
\begin{equation}\label{eq:weak-formulation-transformed}
\forall \tilde{\uv{v}} \in \tilde{\mathcal{V}},
\quad \int_{\tilde{\Omega}} \sym(\tilde{\grad} \tilde{\uv{v}} \cdot \uu{F}^{-1}) : \uuuu{C} : \parens[\big]{\uu{E} +\sym(\tilde{\grad} \tilde{\uv{u}} \cdot \uu{F}^{-1})} J \diff \tilde{\uv{x}}= 0.
\end{equation}

\begin{figure}[!h]
\includegraphics[width=.9\textwidth]{microstructure/grid_adaptation}
\caption{Transformation of a reference uniform grid to a boundary-conforming grid by a mapping $\uv{x} = \uv{\Phi}(\tilde{\uv{x}})$.}\label{fig:grid-adaptation}
\end{figure}

\paragraph{Isoparametric transformation}
We restrict the presentation to an isoparametric transformation of the uniform grid~\cite{Ern2004}. Each node $\tilde{\uv{x}}^n_{\tuple{p}} \in \tilde{\Omega}$ of the uniform grid is transported at its physical position $\uv{x}^n_{\tuple{p}} = \tilde{\uv{x}}^n_{\tuple{p}} + \widetilde{\uv{u}_{\Phi}}(\tilde{\uv{x}}^n_{\tuple{p}})$. The periodic nodal mapping displacements are stored in a vector $\widetilde{\vecdof{u}_{\Phi}}$. A point $\tilde{\uv{x}}$ in the element $\tilde{\mathcal{E}}^e_{\tuple{p}}$ is transported at $\uv{x}$ in the element $\mathcal{E}^e_{\tuple{p}}$ by interpolation as in~\eqref{eq:fe-interpolation}
\begin{equation}\label{eq:fe-interpolation-node-transformation}
\uv{x} = \tilde{\uv{x}} + \widetilde{\uv{u}_{\Phi}}(\tilde{\uv{x}}) = \tilde{\uv{x}} + \sum_{\tuple{s} \in \tuple{S}} \parens*{\widetilde{\vecdof{u}_{\Phi}} }_{\tuple{p}+\tuple{s}} \phi^e_{\tuple{s}}(\uv{y})
\quad \text{where}
\quad \uv{y} = \tilde{\uv{x}} - \tuple{p} \uv{h} \in \tilde{\mathcal{E}}^e.
\end{equation}
The transformation gradient at quadrature point is obtained from $\tilde{\grad} \widetilde{\vecdof{u}_\Phi} = \tilde{\matdof{B}} \widetilde{\vecdof{u}_\Phi}$ following the definition of the discrete gradient operator~\eqref{eq:fe-gradient-quadrature} on the uniform grid as
\begin{equation}\label{eq:fe-transformation-gradient}
\vecdof{F}^q_{\tuple{p}} = \uu{I} + \parens{\tilde{\grad} \widetilde{\vecdof{u}_\Phi}}^q_{\tuple{p}}.
\end{equation}
The fluctuation displacement $\uv{u}$ is interpolated from nodal values $\tilde{\vecdof{u}}$ in the computational grid by
\begin{equation}\label{eq:fe-interpolation-variable-transformation}
\uv{u}(\uv{x}) = \tilde{\uv{u}}(\tilde{\uv{x}}) = \sum_{\tuple{s} \in \tuple{S}} \tilde{\vecdof{u}}_{\tuple{p}+\tuple{s}} \phi^e_{\tuple{s}}(\uv{y})
\quad \text{where}
\quad \uv{y} = \tilde{\uv{x}} - \tuple{p} \uv{h} \in \tilde{\mathcal{E}}^e.
\end{equation}
Eventually, at quadrature points $\uv{x}^q_{\tuple{p}} = \tilde{\uv{x}}^q_{\tuple{p}} + \widetilde{\uv{u}_\Phi}(\tilde{\uv{x}}^q_{\tuple{p}})$, the gradient of the displacement follows from~\eqref{eq:gradient-transformed} as
\begin{equation}\label{eq:fe-discrete-gradient-transformed}
\parens*{\grad \vecdof{u}}^{q}_{\tuple{p}} = \tilde{\grad} \tilde{\vecdof{u}}^{q}_{\tuple{p}} \cdot \parens{\vecdof{F}^q_{\tuple{p}} }^{-1},
\end{equation}
where $\tilde{\grad} \tilde{\vecdof{u}} = \tilde{\matdof{B}} \tilde{\vecdof{u}}$ from~\eqref{eq:fe-gradient-quadrature}.

\paragraph{Transformed discrete problem}
The Gauss quadrature approximation~\eqref{eq:gauss-quadrature} of the pulled-back weak formulation~\eqref{eq:weak-formulation-transformed} leads to the discrete version of~\eqref{eq:stiffness-operator} for a transformed grid
\begin{equation}\label{eq:fe-discrete-problem-transformed}
\tilde{\matdof{K}} \tilde{\vecdof{u}} = \tilde{\vecdof{f}}_{\uu{E}}
\quad \text{where}
\quad \tilde{\matdof{K}} = \tilde{\matdof{B}}^{\mathsf{T}} \tilde{\matdof{W}} \tilde{\matdof{C}} \tilde{\matdof{B}}
\quad \text{and}
\quad \tilde{\vecdof{f}}_{\uu{E}} =- \tilde{\matdof{B}}^{\mathsf{T}} \tilde{\matdof{W}} \tilde{\matdof{C}} \tilde{\vecdof{E}},
\end{equation}
where $\tilde{\matdof{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~\eqref{eq:fe-discrete-problem-transformed}, $\tilde{\matdof{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{\uuuu{C}}$ defined as~\cite{Bellis2024}:
\begin{equation}\label{eq:transformed-stiffness}
\tilde{C}_{\mathit{ijk\ell}} = J F^{-1}_{jp} C_{\mathit{pikq}} F^{-T}_{q \ell},
\end{equation}
where $\uu{F}$ and $J = \det(\uu{F})$ are obtained at each Gauss integration point from~\eqref{eq:fe-transformation-gradient}. The fourth order tensor $\tilde{\uuuu{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~\eqref{eq:fe-discrete-problem-transformed} from the pulled-back stiffness and the pulled-back macroscopic strain $\tilde{\vecdof{E}}$ defined at quadrature points as the second order tensor
\begin{equation}\label{eq:transported-macro-strain}
\tilde{\vecdof{E}}^{q}_{\tuple{p}} = \uu{E} \cdot \vecdof{F}^q_{\tuple{p}}.
\end{equation}
The strain $\tilde{\vecdof{E}}$ is neither symmetric nor uniform. In practice, the explicit computation of $\tilde{\uuuu{C}}$ and $\tilde{\vecdof{E}}$ can be avoided by using an operator based implementation as detailed in Section~\ref{sec:implementation}.

\subsection{Differential operators and Green function in the Fourier space}\label{sec:differential-green-fourier}

\paragraph{Discrete Fourier Transform (DFT)}
The DFT of the array $\vecdof{u}$ at Fourier index $\tuple{k} \in \mathcal{I}(\tuple{N})$ is defined by:
\begin{equation}\label{eq:DFT-primal-fem}
\hat{\vecdof{u}}_{\tuple{k}} = \dfrac{1}{N_p} \sum_{\tuple{p} \in \mathcal{I}(\tuple{N})} \vecdof{u}_{\tuple{p}} \exp \parens{-\imath 2\pi \tuple{k} \cdot \tuple{p} / \tuple{N} }.
\end{equation}
In practice, symmetry properties of the real-to-complex DFT allow restricting the set of Fourier indices to $\hat{\mathcal{I}}(\tuple{N}) = \{0, \dotsc, N_1 \divisionsymbol 2\} \times \{0, \dotsc, N_2 -1\} \bracks[\big]{\times \{0, \dotsc, N_3 -1\}}$ where $\divisionsymbol$ denotes integer division.\footnote{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.}

\paragraph{DFT of discrete gradient}
From~\eqref{eq:fe-gradient-quadrature}, the DFT of the array of discrete gradient $\{\grad \vecdof{u}^q_{\tuple{p}}\}_{\tuple{p} \in \mathcal{I}(\tuple{N})}$ is consequently
\begin{equation}\label{eq:DFT-fe-gradient-primal}
\widehat{\grad \vecdof{u}}^q_{\tuple{k}} = \dfrac{1}{N_p} \sum_{\tuple{p} \in \mathcal{I}(\tuple{N})} \sum_{\tuple{s} \in \tuple{S}} \vecdof{u}_{\tuple{p}+\tuple{s}} \otimes \matdof{B}^{q}_{\tuple{s}} \exp \parens{-\imath 2\pi \tuple{k} \cdot \tuple{p} / \tuple{N} }.
\end{equation}
Application of the DFT shift theorem on Cartesian index $\tuple{p}$ for shift $\tuple{s}$ leads to the expression of the discrete gradient in the Fourier domain:
\begin{equation}\label{eq:DFT-fe-gradient}
\widehat{\grad \vecdof{u}}^q_{\tuple{k}} = \hat{\vecdof{u}}_{\tuple{k}} \otimes \imath \uv{\xi}^{q}_{\tuple{k}}
\quad \text{and}
\quad
\imath \uv{\xi}^{q}_{\tuple{k}} =
\sum_{\tuple{s} \in \tuple{S}} \matdof{B}^{q}_{\tuple{s}}\exp \parens*{\imath 2\pi \tuple{k}\tuple{s} / \tuple{N} }.
\end{equation}
By direct analogy with the continuous Fourier transform, the interpretation of $\uv{\xi}^{q}_{\tuple{k}}$ is the discrete wave vector associated to the discrete gradient $ \grad \vecdof{u}^{q}_{\tuple{p}}$ defined in~\eqref{eq:fe-gradient-quadrature}.

\paragraph{DFT of discrete divergence}
Similarly, the DFT of the array $\div \uu{\upsigma}$ defined in~\eqref{eq:fe-divergence-nodes} is
\begin{equation}\label{eq:DFT-fe-divergence}
\widehat{\div \uu{\upsigma}}_{\tuple{k}} = \sum_{q=1}^{N_q} w_q \widehat{\div \uu{\upsigma}}^{q}_{\tuple{k}}
\quad \text{where}
\quad \widehat{\div \uu{\upsigma}}^{q}_{\tuple{k}} = \hat{\uu{\upsigma}}^q_{\tuple{k}} \cdot \imath \ov{\uv{\xi}}^{q}_{\tuple{k}},
\end{equation}
where $\ov{\uv{\xi}}^{q}_{\tuple{k}}$ is the complex conjugate of the discrete wave vector $\uv{\xi}^{q}_{\tuple{k}}$ defined in~\eqref{eq:DFT-fe-gradient}.

\paragraph{Discrete Green function}
In the reference problem~\eqref{eq:reference-stiffness-operator}, the heterogeneous material properties $\uuuu{C}$ are replaced by a homogeneous one $\uuuu{C}_0$. The discretization leads to the reference stiffness matrix
\begin{equation}\label{eq:reference-stiffness-fe}
\matdof{K}_0 = \matdof{B}^{\mathsf{T}} \matdof{W} \matdof{C}_0 \matdof{B}.
\end{equation}
By construction, $\matdof{K}_0$ is block-circulant, and its DFT is derived from the discrete wave vectors defined in~\eqref{eq:DFT-fe-gradient}. The discrete Green operator $\matdof{G}_0$ is the (pseudo-)inverse of $\matdof{K}_0$, see~\eqref{eq:green-function}. The DFT of the discrete Green function can then be computed as follows.

\paragraph{DFT of discrete Green function}
Let the nodal force vector $\vecdof{f}_0 = \matdof{K}_0 \vecdof{u}$. From the expression of the DFT of the discrete differential operators~\eqref{eq:DFT-fe-gradient} and~\eqref{eq:DFT-fe-divergence}, the DFT of $\matdof{K}_0$ operates diagonally at Fourier index $\tuple{k}$ as
\begin{equation}\label{eq:fourier-green-operation}
\widehat{\vecdof{f}_0}_{\tuple{k}} = \widehat{\matdof{K}_0 \vecdof{u}}_{\tuple{k}} = - \sum_{q=1}^{N_q} w_q \parens*{\uuuu{C}_0 : \parens{\hat{\vecdof{u}}_{\tuple{k}} \otimes \imath \uv{\xi}^{q}_{\tuple{k}}}} \cdot \imath \ov{\uv{\xi}}^{q}_{\tuple{k}} =\widehat{\matdof{K}_0}_{\tuple{k}} \cdot \hat{\vecdof{u}}_{\tuple{k}},
\end{equation}
where $\widehat{\matdof{K}_0}_{\tuple{k}}$ is the acoustic tensor defined by
\begin{equation}\label{eq:acoustic-tensors}
\widehat{\matdof{K}_0}_{\tuple{k}} = \sum_{q=1}^{N_q} w_q \parens{\ov{\uv{\xi}}^{q}_{\tuple{k}} \cdot \uuuu{C}_0 \cdot \uv{\xi}^{q}_{\tuple{k}}}.
\end{equation}
When $\uuuu{C}_0$ is isotropic with Lamé parameters $(\lambda_0, \mu_0)$,
\begin{equation}\label{eq:acoustic-tensors-iso}
\ov{\uv{\xi}}^{q}_{\tuple{k}} \cdot \uuuu{C}_0 \cdot \uv{\xi}^{q}_{\tuple{k}} = \lambda_0 \ov{\uv{\xi}}^{q}_{\tuple{k}} \otimes \uv{\xi}^{q}_{\tuple{k}} + \mu_0 (\ov{\uv{\xi}}^{q}_{\tuple{k}} \cdot \uv{\xi}^{q}_{\tuple{k}} \uu{I} + \uv{\xi}^{q}_{\tuple{k}} \otimes \ov{\uv{\xi}}^{q}_{\tuple{k}}).
\end{equation}
The DFT of the discrete Green function is then the real and symmetric second order tensor
\begin{equation}\label{eq:DFT-green-function-fe}
\widehat{\matdof{G}_0}_{\tuple{k}} = \widehat{\matdof{K}_0}^{-1}_{\tuple{k}}
\quad \text{if $\tuple{k} \neq (0, \dotsc, 0)$ else $0$},
\end{equation}
such that $\hat{\vecdof{u}}_{\tuple{k}} = \widehat{\matdof{G}_0}_{\tuple{k}} \widehat{\vecdof{f}_0}_{\tuple{k}}$.

When $N_q = 1$ (e.g., $Q_1$ type finite elements with reduced integration~\cite{Willot2015,Schneider2017}), the inverse of the acoustic tensor $\widehat{\matdof{K}_0}_{\tuple{k}}$ is known analytically from the discrete wave vector $\uv{\xi}_{\tuple{k}} = \uv{\xi}^{1}_{\tuple{k}}$ as
\begin{equation}\label{eq:DFT-green-function-iso}
\widehat{\matdof{G}_0}_{\tuple{k}} =
\dfrac{1}{\mu_0 \norm{\uv{\xi}_{\tuple{k}}}^2} \parens[\Bigg]{\uu{I} - \dfrac{\lambda_0 + \mu_0}{\lambda_0 + 2\mu_0} \dfrac{\ov{\uv{\xi}_{\tuple{k}}} \otimes \uv{\xi}_{\tuple{k}}}{\norm{\uv{\xi}_{\tuple{k}}}^2} }
\quad \text{if $\uv{\xi}_{\tuple{k}} \neq 0$ else $0$},
\end{equation}
and $\ov{\uv{\xi}_{\tuple{k}}} \otimes \uv{\xi}_{\tuple{k}}$ is real and symmetric. Otherwise, the inversion of a system of $d$ linear equations must be carried out numerically at each Fourier index $\tuple{k} \in \hat{\mathcal{I}}(\tuple{N})$.

\section{Fast iterative solvers with Green preconditioning}\label{sec:green-preconditioning-eigenvalues}

\subsection{Preconditioned conjugate gradient method}\label{sec:solvers}

The homogenization problem discretized on either the uniform or transformed grid leads to the stiffness matrices $\matdof{K}$~\eqref{eq:fe-discrete-problem} or $\tilde{\matdof{K}}$~\eqref{eq:fe-discrete-problem-transformed} 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~\eqref{eq:fe-discrete-problem} or~\eqref{eq:fe-discrete-problem-transformed} can thus be solved using conjugate gradient (CG) \cite{Saad2003, Liesen2012, vanderVorst2003}, see e.g.~\cite{Zeman2010,Kabel2014,Schneider2015,Mishra2016}, or fast gradient solvers~\cite{Schneider2019,Schneider2020}.

To simplify the discussion of this subsection, $\matdof{K}$ denotes indifferently the stiffness matrix on the uniform or transformed grid in what follows. The eigenvalue spectrum of $\matdof{K}$ is spread, non clustered, and the condition number increases quadratically with grid refinement~\cite{Ern2004}, so that the preconditioned conjugate gradient (PCG) is recommended over CG~\cite{Braess2007}. The preconditioner should be a matrix $\matdof{M}$ sufficiently close to $\matdof{K}$ to improve convergence, yet easily invertible in the row space of $\matdof{K}$. Green preconditioning refers to the choice $\matdof{M} = \matdof{K}_0$ defined by~\eqref{eq:reference-stiffness-fe}, whose pseudo-inverse in the row space of $\matdof{K}$ is the discrete Green function $\matdof{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 $\vecdof{z}_i$ be the $i^{\mathrm{th}}$ eigenvalue and eigenvector of the generalized eigenvalue problem for matrix $\matdof{K}$ with preconditioner $\matdof{K}_0$ (restricted to the subspace orthogonal to the shared kernel of $\matdof{K}$ and $\matdof{K}_0$), i.e.\ $\matdof{K}\vecdof{z}_i = \lambda_i\matdof{K}_0\vecdof{z}_i$. Eigenvectors are normalized such that $\vecdof{z}_i^{\mathsf{T}} \matdof{K}_0 \vecdof{z}_j = \delta_{ij}$. The modal weight of the initial residual $\vecdof{r}^{(0)} = \vecdof{f}_{\uu{E}} - \matdof{K} \vecdof{u}^{(0)}$ is defined as
\begin{equation}\label{eq:conjugate-gradient-weight}
w_i = \dfrac{(\vecdof{r}^{(0)} \cdot \vecdof{z}_i)^2}{\lambda_i}.
\end{equation}
The initial residual is equal to the r.h.s.\ $\vecdof{f}_{\uu{E}}$ for the initial guess $\vecdof{u}^{(0)} = 0$. At the $k$-th iteration of the PCG, the error on the solution $\vecdof{u}$ of the iterate $\vecdof{u}^{(k)}$ is
\begin{equation}\label{eq:conjugate-gradient-error}
\norm{\vecdof{u}^{(k)} - \vecdof{u}}_{\matdof{K}} = \norm{\vecdof{r}^{(k)}}_{\matdof{K}^{-1}} =\sum_{i=1}^{d N_p} w_i \, \bracks[\big]{P^{(k)}(\lambda_i) }^2,
\end{equation}
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~\cite{Jennings1977,Axelsson1986,vanderSluis1986}. Then convergence of PCG depends on the:
\begin{itemize}
\item condition number $\kappa = \lambda_{\max} / \lambda_{\min}$;
\item clustering of eigenvalues with non negligible modal weights of the first residual.
\end{itemize}
Throughout the paper, the iterative solution is computed with the PCG and the first residual $\vecdof{r}^{(0)}$ is set equal to the r.h.s.\ vector. Iterations are run until the (too) stringent stopping criterion $\norm{\vecdof{r}^{(k)}}_{\matdof{G}_0} < 10^{-10} \norm{\vecdof{r}^{(0)}}_{\matdof{G}_0}$ is met. The preconditioner norm of the residual $\norm{\vecdof{r}^{(k)}}_{\matdof{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 $\norm{\vecdof{r}^{(k)}}_{\matdof{K}^{-1}}$ appearing in~\eqref{eq:conjugate-gradient-error}.

\subsection{Green preconditioning on uniform or structured computational grids}\label{sec:green-preconditioning}

\paragraph{General results}
FFT-based methods draw a large part of their efficiency from Green (or Laplacian) preconditioning~\cite{Ladecky2023}. The principle of Green preconditioning is to use the solution to the problem~\eqref{eq:reference-stiffness-operator} defined with a uniform reference medium to precondition the homogenization problem~\eqref{eq:homogenization-problem}. The eigenvalues of the Green preconditioned operator have been analyzed by~\cite{Nielson2008, Gergelits2019, Gergelits2020, Gergelits2022} for scalar elliptic problems of the type $\div c \grad u$ where $c$ is a heterogeneous field of scalar material properties. For a FEM discretization and preconditioning by the inverse of $\div c_0 \grad $, 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 $\uv{u}$ and tensor material properties $\uuuu{C}$ by~\cite{Pultarova2021, Ladecky2021}. 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 $\uuuu{C}_0^{-1}:\uuuu{C}$ in the support function of the elements.

\paragraph{Uniform grid}
For a discretization on a uniform grid as in~\eqref{eq:fe-discrete-problem}, i.e.~with no grid adaptation, the Green preconditioner for $\matdof{K}$ is $\matdof{M} = \matdof{K}_0$ defined in~\eqref{eq:reference-stiffness-fe}. The pseudo-inverse of $\matdof{K}_0$ is the discrete Green operator $\matdof{G}_0$. The latter can be applied in the Fourier space from~\eqref{eq:fourier-green-operation}, \eqref{eq:acoustic-tensors}, \eqref{eq:DFT-green-function-fe} and use of FFT. This is used e.g.~in~\cite{Schneider2017, Leute2022, Ladecky2023}, and is equivalent to solving the strain based version $\vv{\Gamma}_0 \ast \uuuu{C} :\uu{\epsilon} = - \vv{\Gamma}_0 \ast \uuuu{C}:\uu{E}$, where $\vv{\Gamma}_0 = - \matdof{B} \matdof{K}_0 \matdof{B}^{\mathsf{T}} \matdof{W}$, provided a Green operator norm of the residual is adopted~\cite{Schneider2017,Leuschner2018,Ladecky2023}.

\paragraph{Transformed grid}
In turn for the problem~\eqref{eq:fe-discrete-problem-transformed} discretized on a structured grid and pulled back on a uniform grid, the true Green preconditioner for $\tilde{\matdof{K}}$ would be $\matdof{M} = \widetilde{\matdof{K}_0}$, the stiffness matrix obtained from the finite element discretization on the transformed structured grid with uniform stiffness $\uuuu{C}_0$. The latter can be constructed as $\widetilde{\matdof{K}_0} = \tilde{\matdof{B}}^{\mathsf{T}} \tilde{\matdof{W}} \widetilde{\matdof{C}_0} \tilde{\matdof{B}}$ where $\widetilde{\matdof{C}_0}$ is built from the pulled-back reference stiffness tensor $\widetilde{\uuuu{C}_0}$ defined by replacing $\uuuu{C}$ by $\uuuu{C}_0$ in~\eqref{eq:transformed-stiffness}. In that case, the bounds on the spectrum of eigenvalues of the preconditioned problem would be determined from local eigenvalues of $\uuuu{C}_0^{-1}:\uuuu{C}$.

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

\subsection{Bounds on the eigenvalue spectrum of the preconditioned operator}\label{sec:eigenvalues-spectrum}

\paragraph{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 \emph{a priori}~\cite{Pultarova2021,Ladecky2021}. The bounds on the eigenvalues of $\matdof{G}_0 \tilde{\matdof{K}}$ (or of the generalized eigenvalue problem for $\tilde{\matdof{K}}$ and $\matdof{K}_0$) can be computed as follows:
\begin{enumerate}
\item for the reference node in pattern $\tuple{p}$ of the structure, compute the largest ($(\uu{\uplambda}^+)_{\tuple{p}}$) and lowest ($(\uu{\uplambda}^-)_{\tuple{p}}$) of all eigenvalues of $\uuuu{C}_0^{-1}:\tilde{\uuuu{C}}^q_{\tuple{p} - \tuple{s}}$ among all quadrature points in the support of that node, and store them in two vectors $\uu{\uplambda}^\pm$ of size $N_p$;
\item sort $\uu{\uplambda}^\pm$ by increasing values, with repetition of each value $d$ times.
\end{enumerate}
The above bounds hold for elements with exact integration, in the shared row space of $\tilde{\matdof{K}}$ and $\matdof{K}_0$. There are $d$ zero eigenvalues of $\tilde{\matdof{K}}$ and $\matdof{K}_0$, and thus $\matdof{G}_0 \tilde{\matdof{K}}$, corresponding to rigid body translations (rigid body rotations are excluded due to periodicity). For underintegrated elements, additional zero-energy hourglass modes arise~\cite{Willot2015,Schneider2022}.

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

\paragraph{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 $\uuuu{C}$ is uniform and a scalar multiple of the reference medium; i.e.\ $\uuuu{C} = \alpha \uuuu{C}_0$. Then, each of the $N_\alpha$ basis functions with support in $\Omega_\alpha$ yields $d$ eigenvectors of $\matdof{G}_0 \matdof{K}$ with eigenvalue $\alpha$~\cite{Ladecky2021}. Hence $dN_\alpha$ eigenvalues cluster exactly at $\alpha$. The modal weight of the right-hand side defined by~\eqref{eq:conjugate-gradient-weight} for the corresponding modes is zero, because $\div \uuuu{C} : \uu{E} = 0$ at these nodes. These modes correspond to type~I described in~\cite{Bellis2025} in the case of isotropic conductivity. This situation is extremely favorable for the convergence of the PCG, as only the \emph{effective} spectrum of eigenvalues with non-zero weight is involved in the resolution, see~\eqref{eq:conjugate-gradient-error}. Since eigenvectors are $\matdof{K}$ and $\matdof{K}_0$-orthogonal, the set of modes with non-zero weight is included in the subspace generated by $\matdof{G}_0 \vecdof{f}$, where $\vecdof{f}$ are vectors of nodal forces supported in regions where $\uuuu{C}$ is not a scalar multiple of $\uuuu{C}_0$, including nodes at interfaces. Such eigenvalue clusters have also been observed for a strain based method with discretization on trigonometric polynomials and numerical integration~\cite{Mishra2016}, and in the case of isotropic conductivity~\cite{Bellis2025}.

Further, when the stiffness $\uuuu{C}$ is bounded by $\alpha^- \uuuu{C}_0 \leqslant \uuuu{C} \leqslant \alpha^+ \uuuu{C}_0$ (in the sense of positive definite quadratic forms) in a domain $\Omega_\alpha$ containing $N_\alpha$ basis functions, then the bounds constructed from the procedure in~\cite{Gergelits2020,Pultarova2021,Ladecky2021} guarantee the existence of a $dN_\alpha$ of eigenvalues within the interval $[\alpha^-, \alpha^+]$, which is favorable when $\alpha^-$ and $\alpha^+$ are close.

\begin{figure}[!h]
\includegraphics[width=0.295\textwidth]{microstructure/microstructure_circle_uniform_32}~
\includegraphics[width=0.295\textwidth]{microstructure/microstructure_square_to_circle_32}~
\includegraphics[width=0.295\textwidth]{microstructure/microstructure_circle_transformed_32}~~~\mbox{}

\vskip .25\baselineskip

\begin{subfigure}{\textwidth}
\includegraphics[width=\textwidth]{spectrum/spectruma_rev.pdf}
\caption{Lamé parameters $\lambda_i = 0.01 \lambda_m$, $\mu_i = 0.01 \mu_m $}\label{fig:spectrum-transformation-circle-elasticity_a}
\end{subfigure}

\vskip .5\baselineskip

\begin{subfigure}{\textwidth}
\includegraphics[width=\textwidth]{spectrum/spectrumb_rev.pdf}
\caption{Lamé parameters $\lambda_i = 0.1 \lambda_m$, $\mu_i = 0.01 \mu_m$}\label{fig:spectrum-transformation-circle-elasticity_b}
\end{subfigure}
\caption{Effect of grid adaptation on the spectrum and effective condition number $\kappa^{\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~\cite{Pultarova2021,Ladecky2021}, dots are computed eigenvalues, colored by modal weight~\eqref{eq:conjugate-gradient-weight}, gray vertical lines are weights.}\label{fig:spectrum-transformation-circle-elasticity}
\end{figure}

\paragraph{Effect of grid adaptation}
Let's illustrate these properties of the spectrum on a simple example. Figure~\ref{fig:spectrum-transformation-circle-elasticity} 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 $\uuuu{C}_m$ has Lamé parameters $\lambda_m=1$, $\mu_m = 0.5$. Two cases are considered for the inclusion stiffness $\uuuu{C}_i$:
\begin{enumerate}\alphenumi
\item \label{effect_case_a} $\uuuu{C}_i = 0.01\uuuu{C}_m$ (2\textsuperscript{nd} row of Figure~\ref{fig:spectrum-transformation-circle-elasticity});
\item \label{effect_case_b} $\lambda_i = 0.1 \lambda_m, \; \mu_i = 0.01 \mu_m $ (3\textsuperscript{rd} row of Figure~\ref{fig:spectrum-transformation-circle-elasticity}).
\end{enumerate}
In both cases, the stiffness of the reference medium for the Green preconditioner is chosen as $\uuuu{C}_0 = \sqrt{\uuuu{C}_i : \uuuu{C}_m}$. 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~\eqref{effect_case_a}, $\uuuu{C}_i$ and $\uuuu{C}_m$ are scalar multiples of~$\uuuu{C}_0$, while they are not in case~\eqref{effect_case_b}.

\goodbreak

Three types of grid transformation are investigated, from left to right in Figure~\ref{fig:spectrum-transformation-circle-elasticity}:
\begin{enumerate}[leftmargin=5em,labelwidth=\widthof{(center)},align=left]
\item[(left)] \emph{No transform}. Elements in the uniform grid are attributed to the inclusion or the matrix based on the position of their center of mass.
\item[(center)] \emph{Algebraic transform}. A square of edge length $2\tilde{R} = 5/8 \, L$ in the computational domain $\tilde{\Omega}=[-L/2, L/2]^d$ is transformed to a circle of radius $R = 3/8 \, L$ in the physical domain $\Omega=[-L/2, L/2]^d$. Any point $\tilde{\uv{x}} \in \tilde{\Omega}$ is radially transformed to $\uv{x} = \rho \tilde{\uv{x}} $ with $\rho = 1$ if $\tilde{\uv{x}} = 0$ else
\[
\rho =
\begin{cases}
	1 + \tilde{r}(R - R_i) / R_i^2 & \text{if $0 <\tilde{y} \leq \tilde{R}$},
\\	\bracks[\big]{R + (R_m - R) (\tilde{r} - R_i) / (R_m - R_i)} \big/ \tilde{r} & \text{else},
\end{cases}
\]
where $\tilde{r} = \norm{\tilde{\uv{x}}}$, $\tilde{y} = \max_{i=1\dots d} \abs{\tilde{x}_i}$, $R_m = t L/2$ and $R_i = t \tilde{R}$ with $t = \tilde{r} / \tilde{y}$.
\item[(right)] \emph{Boundary node projection transform}. Nodes of the uniform grid found at the boundary between elements attributed to matrix and inclusion in the \emph{no transform} procedure are projected on the actual boundary of the inclusion. Other nodes are not transformed.
\end{enumerate}

Bi-linear elements with four quadrature points are used to discretize the problem on a $32{\times}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 Figure~\ref{fig:spectrum-transformation-circle-elasticity}:
\begin{enumerate}[leftmargin=5em,labelwidth=\widthof{(center)},align=left]
\item[(left)] \emph{No transform}. The eigenvalue spectrum is bounded by $0.1$ and $10$. In case~\eqref{effect_case_a}, 2\textsuperscript{nd} row, eigenvalues cluster exactly at these two values since $\uuuu{C}_i = 0.1\uuuu{C}_0$ and $\uuuu{C}_m = 10\uuuu{C}_0$. The corresponding eigenvectors are respectively supported by nodes exactly in the inclusion or in the matrix, and the modal weights~\eqref{eq:conjugate-gradient-weight} 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~\cite{Bellis2025}. Further, numerical investigation of fields confirm that modes with non-zero weight belong to the subspace generated by $\matdof{G}_0 \vecdof{f}$, where $\vecdof{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 Figure~\ref{fig:spectrum-transformation-circle-elasticity}\eqref{fig:spectrum-transformation-circle-elasticity_a}. The effective spectrum of these 48 modes with non-zero modal weight, corresponding to dots colored in green, has a condition number $\kappa^{\eff} = \lambda^{\eff}_{\max} / \lambda^{\eff}_{\min} = 9.4$. This value is roughly 10 times smaller than that of the full spectrum, which is equal to $\kappa=100$. In case~\eqref{effect_case_b}, these two clusters diffuse and continuously spread in the intervals formed by the eigenvalues of $\uuuu{C}_0^{-1}:\uuuu{C}_i$ and $\uuuu{C}_0^{-1}:\uuuu{C}_m$, and some have non-zero modal weight.\footnote{In Figure~\ref{fig:spectrum-transformation-circle-elasticity}, 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.} This situation thus differs from the isotropic conductivity situation studied by~\cite{Bellis2025}: since $\uuuu{C}$ is not a scalar multiple of $\uuuu{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~\eqref{effect_case_a} to 75 in case~\eqref{effect_case_b}.
\item[(center)] \emph{Algebraic transform}. All elements are transformed. In both cases~\eqref{effect_case_a} and~\eqref{effect_case_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~\eqref{effect_case_a} and 202 in case~\eqref{effect_case_b}.
\item[(right)] \emph{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~\eqref{effect_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 \times 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~\eqref{effect_case_a} and 95 for case~\eqref{effect_case_b}.
\end{enumerate}
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 \emph{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~\cite{Bellis2024,Zecevic2026}.

\section{Implementation}\label{sec:implementation}

\paragraph{Iterative solver}
The implementation relies on the standard PCG algorithm~\cite{Saad2003,vanderVorst2003} with an operator based, matrix free implementation. Main operations are: computation of the right-hand side, application of the stiffness operator $\matdof{K}$ or $\tilde{\matdof{K}}$, and application of the Green preconditioner~$\matdof{G}_0$. 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.

\paragraph{Differential operators}
The discrete gradient~\eqref{eq:fe-gradient-quadrature} and divergence~\eqref{eq:fe-divergence-nodes} are short kernel convolutions. These operations are performed in the space domain using stencils built from the $\matdof{B}^{q}_{\tuple{s}}$ defined in~\eqref{eq:fe-gradient-quadrature}. The quadrature weights are fused in the divergence operation~\cite[Section~5.1]{Ladecky2023}.

\begin{algorithm}
\begin{algorithmic}[1]
\State $\tilde{\vecdof{f}} \gets 0$
\ForAll{$q \in 1\dots N_q$} \Comment{$q^{\mathrm{th}}$-quadrature point}
\State \label{alg:F}
$ \vecdof{F}^q \gets \uu{I} + \parens{\tilde{\grad} \widetilde{\vecdof{u}_\Phi}}^q$
\Comment{$q^{\mathrm{th}}$-transformation gradient}
\State \label{alg:epsilon}
$\uu{\upepsilon}^q \gets
\sym\parens[\big]{\parens{\tilde{\grad} \tilde{\vecdof{u}}}^{q} \cdot \parens{\vecdof{F}^q }^{-1} } $
\Comment{$q^{\mathrm{th}}$-transformed gradient}
\State \label{alg:sigma}
$ \uu{\upsigma}^q \gets \uuuu{C}^q : \uu{\upepsilon}^q $
\Comment{constitutive model}
\State \label{alg:div}
$\tilde{\vecdof{f}} \gets \tilde{\vecdof{f}} + w_q \widetilde{\div} \parens*{\det(\vecdof{F}^q) \uu{\upsigma}^q \cdot (\vecdof{F}^q)^{-T}}$
\Comment{transformed divergence}
\EndFor
\end{algorithmic}
\caption{$\tilde{\vecdof{f}} \gets \tilde{\matdof{K}} \tilde{\vecdof{u}}$ --- Application of stiffness operator $\tilde{\matdof{K}}$~\eqref{eq:fe-discrete-problem-transformed} on a displacement-like vector $\tilde{\vecdof{u}}$.}\label{alg:stiffness-operator-transformed}
\end{algorithm}

\paragraph{Transformed stiffness operator}
The transformed stiffness operator $\tilde{\matdof{K}}$ defined in~\eqref{eq:fe-discrete-problem-transformed} applies to a displacement-like vector $\tilde{\vecdof{u}}$ as sketched in Algorithm~\ref{alg:stiffness-operator-transformed}. As compared to an implementation on a uniform grid (see e.g.~\cite[Algorithm~2]{Schneider2017}), the only modification to account for grid transformation is the presence of the transformation gradient $\vecdof{F}$, computed in line~\ref{alg:F}. 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 $\matdof{F}^{-1}$ is in general prohibitive,\footnote{For the \emph{boundary node projection transform}, $\widetilde{\vecdof{u}_{\Phi}}$ and its gradient are sparse, so $\matdof{F}^{-1}$ could be cached (not done in the current implementation).} so it is recomputed at each iteration. The constitutive model is applied locally at the current quadrature point $q$ of each pattern in line~\ref{alg:sigma}. Vectors $\tilde{\grad}^q \tilde{\vecdof{u}}$, $\uu{\upepsilon}$ and $\uu{\upsigma}$ can share the same memory. Hence two buffer $\grad \vecdof{u}$-like vectors are required to apply $\tilde{\matdof{K}}$.

\paragraph{Right-hand side}
The r.h.s.~is computed as in Algorithm~\ref{alg:stiffness-operator-transformed}, except line~\ref{alg:epsilon} which is replaced by $\uu{\upepsilon} \gets \uu{E}$. Once the system is solved, the strain in the physical domain is obtained at quadrature point $q$ as $\uu{\upepsilon}^q \gets \sym\parens[\big]{(\tilde{\grad} \tilde{\vecdof{u}})^q \cdot (\vecdof{F}^q)^{-1}} + \uu{E}$.

\paragraph{Green preconditioner}
The FFT operator and its inverse can operate in-place to apply the Green function $\widehat{\matdof{G}_0}$ in Fourier space. When $N_q = 1$, the operation $\widehat{\matdof{G}_0}_\tuple{k} \cdot \hat{\vecdof{f}}_{\tuple{k}}$ is directly applied using~\eqref{eq:DFT-green-function-iso}. Otherwise, at each iteration and Fourier index, the matrix $\widehat{\matdof{K}_0}_{\tuple{k}}$ is assembled on-the-fly and the small $d\times d$ system $\hat{\vecdof{f}}_{\tuple{k}} = \widehat{\matdof{K}_0}_{\tuple{k}} \hat{\vecdof{u}}_{\tuple{k}}$ solved. Due to the symmetry of $\widehat{\matdof{G}_0}$, a viable alternative could be to cache the real and symmetric matrices $\widehat{\matdof{G}_0}_{\tuple{k}}$ for $\tuple{k} \in \{0 \dots N_1 \divisionsymbol 2, \dotsc, 0 \dots N_d \divisionsymbol 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.

\paragraph{Memory requirement}
To sum up, the memory requirement is five displacement-like vectors (including grid nodal displacements $\widetilde{\vecdof{u}_{\Phi}}$), i.e.\ $5d N_p$ floating-point scalars. In addition, at least one integer scalar vector is required to encode microstructure phases. As a rule of 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.

\paragraph{Computer programming}
The algorithms have been implemented using the Julia programming language~\cite{Julia2017} and libraries from the Julia ecosystem~\cite{PencilArrays2021,PencilFFTs2021,Tensors2019} with the FFTW backend~\cite{FFTW2005}. 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~\cite{Bignonnet_2026_18458599}.

\section{Numerical investigations}\label{sec: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.

\begin{figure}
\begin{subfigure}{0.32\textwidth}
\includegraphics[width=\textwidth]{microstructure/diamonds.pdf}
\caption{Diamond inclusions}\label{fig:test-microstructures_a}
\end{subfigure}
~
\begin{subfigure}{0.32\textwidth}
\includegraphics[width=\textwidth]{microstructure/hashin.pdf}
\caption{Composite inclusion}\label{fig:test-microstructures_b}
\end{subfigure}
~
\begin{subfigure}{0.32\textwidth}
\includegraphics[width=\textwidth]{microstructure/microstructure_boolean_spheres_num_1000_radius_0.08_discretization_100.png}
\caption{Boolean model of spheres}\label{fig:test-microstructures_c}
\end{subfigure}
\caption{Test microstructures: \eqref{fig:test-microstructures_a}~diamond inclusions forming a checkerboard along dashed lines ($d=2$); \eqref{fig:test-microstructures_b}~Hashin's composite inclusion ($d = 2$ or $3$); \eqref{fig:test-microstructures_c}~Boolean model of spheres ($d=3$).}\label{fig:test-microstructures}
\end{figure}

\subsection{Diamond inclusions}\label{sec:diamond}

The first test microstructure is shown in Figure~\ref{fig:test-microstructures}\eqref{fig:test-microstructures_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~\ref{sec:conductivity}), to benefit from the analytical solution provided in~\cite{Craster2001}. 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 Figure~\ref{fig:test-microstructures}\eqref{fig:test-microstructures_a}. Hence, in practice the diamonds are shifted by half a pixel to the right w.r.t.\ to the configuration of Figure~\ref{fig:test-microstructures}\eqref{fig:test-microstructures_a} when $Q_1$ elements are used. This precaution is necessary in order to shift the centers of the pixels of the uniform grid away from these lines where numerical issues can arise while evaluating the analytical solution of~\cite{Craster2001}.

\begin{figure}
\begin{subfigure}{.48\textwidth}
\includegraphics[width=\textwidth]{field/diamond_conduction_2_gradux_rev.pdf}
\caption{Temperature gradient ${\epsilon_x}_{\num}$}
\end{subfigure}
\hfill
\begin{subfigure}{.48\textwidth}
\includegraphics[width=\textwidth]{field/diamond_conduction_2_error_gradux_rev.pdf}
\caption{Error on temperature gradient ${\epsilon_x}_{\num} - {\epsilon_x}_{\rmref}$}
\end{subfigure}
\caption{Diamond inclusions in 2D conduction with macroscopic temperature gradient $\uv{E} = \uv{e}_x$.}\label{fig:diamond-field-conduction}
\end{figure}

Simulations are carried out with either bilinear rectangular $Q_1$ elements in the reference grid, with either one or four quadrature points, or linear triangular $P_1$ elements (see Figure~\ref{fig:fem-discretization}). Since inclusions are rotated, rectangular elements on a uniform grid are not boundary-conforming. Similarly, triangular elements following the pattern of Figure~\ref{fig:fem-discretization}\eqref{fig:fem-discretization_c} on a uniform grid conform to only half of the inclusion edges. The grid adaptation strategy used here corresponds to the \emph{boundary node projection transform} presented in Section~\ref{sec:eigenvalues-spectrum}. Simulations are carried out on grids ranging from $20^2$ to $2560^2$ nodes to study $h$-convergence, with unit imposed macroscopic temperature gradient $\uv{E} = \uv{e}_x$. Additional simulations with bilinear rectangular $Q_1$ elements on a regular grid are also carried out using the composite pixel technique~\cite{Kabel2015}. 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~\cite{Kabel2015}.

Figure~\ref{fig:diamond-field-conduction} shows fields of temperature gradient ${\epsilon}_x$ and error $\delta {\epsilon}_x={\epsilon_x}_{\num} - {\epsilon_x}_{\rmref}$ 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 one or four quadrature points. The overall error on the local fields is quantified in Figure~\ref{fig:diamond-convergence-conduction}\eqref{fig:diamond-convergence-conduction_b} by the root mean square error (RMSE) over all element centers $\sqrt{(N_pN_e)^{-1}\sum_{\tuple{p} \in \mathcal{I}(\tuple{N})} \sum_{e=1\dots N_e} J_\tuple{p} \delta \uv{\epsilon}^e_\tuple{p} \cdot \delta \uv{\epsilon}^e_\tuple{p}}$.

\begin{figure}
\begin{subfigure}{.48\textwidth}
\includegraphics[width=\textwidth]{convergence/error_hom_diamond_Conduction_2_rev.pdf}
\caption{$\abs{C^\hom_{\num}-C^\hom_{\rmref}}/C^\hom_{\rmref}$}\label{fig:diamond-convergence-conduction_a}
\end{subfigure}
\hfill
\begin{subfigure}{.48\textwidth}
\includegraphics[width=\textwidth]{convergence/error_field_diamond_Conduction_2_rev.pdf}
\caption{RMSE on temperature gradient field $\uv{\epsilon}$}\label{fig:diamond-convergence-conduction_b}
\end{subfigure}
\caption{Diamond inclusions in 2D conductivity. \eqref{fig:diamond-convergence-conduction_a}~Relative error on homogenized conductivity. \eqref{fig:diamond-convergence-conduction_a}~Root mean square error (RMSE) on temperature gradient field $\uv{\epsilon}$ for $\uv{E} = \uv{e}_x$.}\label{fig:diamond-convergence-conduction}
\end{figure}

The field on the uniform grid with rectangular $Q_1$ elements and one quadrature point presents significant hourglassing, in particular close to diamond tips. Hourglassing is removed by using four 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 $\vecdof{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 Figure~\ref{fig:diamond-convergence-conduction}\eqref{fig:diamond-convergence-conduction_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 2-dimensional composite circular inclusion in Appendix~\ref{sec:hashin-2d} are well behaved despite excessive element transformation. For quadrangular $Q_1$ elements with one quadrature point, grid adaptation reduces this error by roughly a factor 4 w.r.t.\ the uniform grid, and roughly 30--40 (!) for four quadrature points. For triangular $P_1$ elements, grid adaptation does not improve the error on $\uu{C}^\hom$ for the diamond inclusion case, while a significant improvement is observed for the composite circular inclusion in Appendix~\ref{sec:hashin-2d}. The RMSE on the temperature gradient field $\uv{\epsilon}$ has a sublinear $h$-scaling, and is best for transformed grids. Figure~\ref{fig:PCG-convergence}\eqref{fig:PCG-convergence_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.

\begin{figure}[!h]
\begin{subfigure}{.48\textwidth}
\includegraphics[width=\textwidth]{convergence/convergence_PCG_diamond_Conduction_2_rev.pdf}
\caption{Diamond inclusions}\label{fig:PCG-convergence_a}
\end{subfigure}
\hfill
\begin{subfigure}{.48\textwidth}
\includegraphics[width=\textwidth]{convergence/convergence_PCG_hashin_Elasticity_3_rev.pdf}
\caption{Hashin's composite inclusion}\label{fig:PCG-convergence_b}
\end{subfigure}
\caption{PCG iteration count. \eqref{fig:PCG-convergence_a}~Diamond inclusions in 2D conduction. \eqref{fig:PCG-convergence_b}~Hashin's composite inclusion in 3D linear elasticity.}\label{fig:PCG-convergence}
\end{figure}

\subsection{Hashin's composite inclusion}\label{sec:hashin}

The second test microstructure has now smooth boundaries: it is Hashin's composite inclusion~\cite{Hashin1962,Hashin1964a,Hashin1964b}, see Figure~\ref{fig:test-microstructures}\eqref{fig:test-microstructures_b}. The inclusion core ($i=1$), shell ($i=2$) and the matrix ($i=3$) have uniform isotropic stiffness $\uuuu{C}_i$ with Lamé parameters $\lambda_i, \mu_i$ and bulk modulus $\kappa_i = \lambda_i + 2 \mu_i / d$. When an isotropic load $\uu{E} = \uu{I}$ is applied,\footnote{Plane strain and plane isotropic load is assumed for the cylindrical inclusion in dimension $d=2$.} the composite inclusion is neutral provided the matrix bulk modulus $\kappa_3$ is set to the equivalent bulk modulus of the inclusion:
\begin{equation}\label{eq:hashin-bulk}
\kappa^{\hom} = \kappa_2 \parens*{1 - \dfrac{\beta \alpha \phi}{1 + \alpha \phi} }
\end{equation}
where
\[
\quad \phi = (R_1 / R_2)^d,
\quad \alpha = \dfrac{d(\kappa_2 - \kappa_1)}{(d-1)2\mu_2 + d\kappa_1},
\quad \text{and} \quad \beta = 1 + (d - 1) \dfrac{2\mu_2}{d \kappa_2}.
\]
In that case, $\uu{\epsilon} = \uu{I}$ in the matrix, is uniform in the core, and given in each phase $i$ by
\begin{equation}\label{eq:hashin-bulk-field}
\bracks[\big]{(a_i + b_i / r^d) \uu{I} - db_i / r^d \uv{e}_r \otimes \uv{e}_r} \cdot \uu{E}
\quad \text{where}
\quad r = \norm{\uv{x}} \ \text{and} \ \uv{e}_r = \uv{x} /r,
\end{equation}
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^{\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 $\uuuu{C}_0 = \sqrt{\uuuu{C}_1:\uuuu{C}_2}$. The matrix bulk modulus $\kappa_3$ is set to $\kappa^{\hom}$~\eqref{eq:hashin-bulk}. Its shear modulus $\mu_3$ is set such that $\uuuu{C}_3$ is a scalar multiple of~$\uuuu{C}_0$. Figures~\ref{fig:hashin-field-elasticity-3} and~\ref{fig:hashin-convergence-elasticity-3} present results for 3-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~\cite[Appendix~A]{Gehrig2025b}, 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~\cite{Finel2025}. 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.

\begin{figure}[!h]
\begin{subfigure}{.48\textwidth}
\includegraphics[width=\textwidth]{field/hashin_elasticity_3_epsilonxx_soft_core_contrast_1000_hires.pdf}
\caption{$\epsilon_{xx}$}\label{fig:hashin-field-elasticity-3_a}
\end{subfigure}
\hfill
\begin{subfigure}{.48\textwidth}
\includegraphics[width=\textwidth]{field/hashin_elasticity_3_error_epsilon_soft_core_contrast_1000_hires.pdf}
\caption{$\norm{{\uu{\epsilon}}_{\num} - {\uu{\epsilon}}_{\rmref}}$}\label{fig:hashin-field-elasticity-3_b}
\end{subfigure}
\caption{Hashin's composite inclusion in 3D linear elasticity with $\ov{\uu{\epsilon}} = \uu{I}$. \eqref{fig:hashin-field-elasticity-3_a}~$\epsilon_{xx}$. \eqref{fig:hashin-field-elasticity-3_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.}\label{fig:hashin-field-elasticity-3}
\end{figure}

Figure~\ref{fig:hashin-field-elasticity-3} presents strain fields $\epsilon_{xx}$ on the lowest resolution grid ($20^3$). Without grid adaption, large errors appear (locally exceeding 100\% of $\norm{\uu{E}}$), in particular for elements at the boundary which is poorly approximated by cubes. Again, the \emph{boundary node projection transform} (see Section~\ref{sec:eigenvalues-spectrum}) is applied. Errors on the strain field are drastically reduced after transformation. Using eight quadrature points, the error drops below $\norm{\delta \uu{\epsilon}} < 10\% \norm{\uu{E}}$ in all elements, while it can exceed 20\% in elements at the highly contrasted core-shell boundary while using only one 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 shape of elements. Indeed, hexahedral elements with three faces and seven nodes belonging to the boundary are transformed into tetrahedrons (see e.g.~\cite[Figures~11 and~13]{Zecevic2026} and Figure~\ref{fig:hashin-field-elasticity-3}), up to the curvature of the boundary. As a result, out of the two tetrahedrons defined in the stencils by~\cite{Finel2025}, 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~\eqref{eq:transformed-stiffness} is undefined.

Results of Figure~\ref{fig:hashin-field-elasticity-3} are confirmed upon grid refinement in Figure~\ref{fig:hashin-convergence-elasticity-3}\eqref{fig:hashin-convergence-elasticity-3_b}. The RMSE on the strain field, computed as $\sqrt{N_p^{-1}\sum_{\tuple{p} \in \mathcal{I}(\tuple{N})} J_\tuple{p} \delta \uu{\epsilon}_\tuple{p} : \delta \uu{\epsilon}_\tuple{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 (Figure~\ref{fig:hashin-convergence-elasticity-3}\eqref{fig:hashin-convergence-elasticity-3_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~\cite{Babuska1970} 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 2-dimensional conductivity, see Figures~\ref{fig:hashin-field-conduction}--\ref{fig:hashin-convergence-elasticity-2} in Appendix~\ref{sec:hashin-2d}.

\begin{figure}[!h]
\begin{subfigure}{.48\textwidth}
\includegraphics[width=\textwidth]{convergence/error_hom_hashin_Elasticity_3_rev.pdf}
\caption{$\abs[\big]{\kappa^\hom_{\num}-\kappa^\hom_{\rmref}}/\kappa^\hom_{\rmref}$}\label{fig:hashin-convergence-elasticity-3_a}
\end{subfigure}
\hfill
\begin{subfigure}{.48\textwidth}
\includegraphics[width=\textwidth]{convergence/error_field_hashin_Elasticity_3_rev.pdf}
\caption{RMSE on strain field $\uu{\epsilon}$}\label{fig:hashin-convergence-elasticity-3_b}
\end{subfigure}
\caption{Hashin's composite inclusion in 3D linear elasticity. \eqref{fig:hashin-convergence-elasticity-3_a}~Relative error on homogenized bulk modulus. \eqref{fig:hashin-convergence-elasticity-3_b}~Absolute error on strain field $\uu{\epsilon}$ for $\ov{\uv{\epsilon}} = \uu{I}$.}\label{fig:hashin-convergence-elasticity-3}
\end{figure}

The number of iterations for PCG convergence stabilizes under grid refinement for uniform grids, with or without composite voxels, (Figure~\ref{fig:PCG-convergence}\eqref{fig:PCG-convergence_b}), as expected from the properties of the spectrum of the Green preconditioned operator, see Section~\ref{sec:eigenvalues-spectrum}. The \emph{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.


\subsection{Boolean model of spheres}\label{sec:boolean-model}

The last investigated microstructure is a Boolean model of spheres, see Figure~\ref{fig:test-microstructures}. 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 $\mathopen{]}2R, 2R+\delta\mathclose{[}$ where $\delta = \sqrt3 L / 50$. This implies that spheres either overlap, either are separated by at least one voxel at low resolution ($\geq 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 $\uuuu{C}_1 = 0.01 \uuuu{C}_2$, and the reference one $\uuuu{C}_0 = \sqrt{\uuuu{C}_1:\uuuu{C}_2}$. The applied macroscopic strain is $\uu{E} = \uv{e}_x \otimes \uv{e}_x$.

\begin{figure}[b]
\includegraphics[width=0.3\textwidth]{field/boolean_spheres_epsilon_colormap.pdf}

\vskip .5\baselineskip

\begin{subfigure}{0.24\textwidth}
\includegraphics[width=\textwidth]{field/boolean_spheres_50_uni_Q1.png}

\vskip .21\baselineskip

\includegraphics[width=\textwidth]{field/boolean_spheres_250_uni_Q1.png}

\vskip .22\baselineskip

\includegraphics[width=\textwidth]{field/boolean_spheres_250_uni_Q1_zoom.png}
\caption{Uniform grid $N_q=1$}
\end{subfigure}
\hfill
\begin{subfigure}{0.24\textwidth}
\includegraphics[width=\textwidth]{field/boolean_spheres_50_uni_Q8.png}

\vskip .21\baselineskip

\includegraphics[width=\textwidth]{field/boolean_spheres_250_uni_Q8.png}

\vskip .22\baselineskip

\includegraphics[width=\textwidth]{field/boolean_spheres_250_uni_Q8_zoom.png}
\caption{Uniform grid $N_q=8$}
\end{subfigure}
\hfill
\begin{subfigure}{0.24\textwidth}
\includegraphics[width=\textwidth]{field/boolean_spheres_50_trans_Q1.png}

\vskip .21\baselineskip

\includegraphics[width=\textwidth]{field/boolean_spheres_250_trans_Q1.png}

\vskip .22\baselineskip

\includegraphics[width=\textwidth]{field/boolean_spheres_250_trans_Q1_zoom.png}
\caption{Transformed grid $N_q=1$}
\end{subfigure}
\hfill
\begin{subfigure}{0.24\textwidth}
\includegraphics[width=\textwidth]{field/boolean_spheres_50_trans_Q8.png}

\vskip .21\baselineskip

\includegraphics[width=\textwidth]{field/boolean_spheres_250_trans_Q8.png}

\vskip .22\baselineskip

\includegraphics[width=\textwidth]{field/boolean_spheres_250_trans_Q8_zoom.png}
\caption{Transformed grid $N_q=8$}
\end{subfigure}
\caption{Strain field $\epsilon_{xx}$ for $\uu{E} = \uv{e}_x \otimes \uv{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).}\label{fig:boolean-fields-elasticity-3}
\end{figure}

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 Figure~\ref{fig:test-microstructures}. After that procedure, few elements at complex sphere intersection have extremely low 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 Figure~\ref{fig:boolean-fields-elasticity-3} 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.

\begin{figure}[!h]
\begin{subfigure}{0.48\textwidth}
\includegraphics[width=\textwidth]{convergence/convergence_PCG_booleancircles_Elasticity_3.pdf}
\caption{PCG iteration count}\label{fig:boolean-convergence-elasticity-3_a}
\end{subfigure}
\hfill
\begin{subfigure}{0.48\textwidth}
\includegraphics[width=\textwidth]{convergence/hom_booleancircles_Elasticity_3.pdf}
\caption{Homogenized stiffness $C^\hom_{1111}$}\label{fig:boolean-convergence-elasticity-3_b}
\end{subfigure}
\caption{Boolean model of spheres in 3D linear elasticity. \eqref{fig:boolean-convergence-elasticity-3_a}~PCG iteration count. \eqref{fig:boolean-convergence-elasticity-3_b}~Homogenized stiffness.}\label{fig:boolean-convergence-elasticity-3}
\end{figure}

Figure~\ref{fig:boolean-convergence-elasticity-3}\eqref{fig:boolean-convergence-elasticity-3_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 eight or one quadrature points respectively. In turn, Figure~\ref{fig:boolean-convergence-elasticity-3}\eqref{fig:boolean-convergence-elasticity-3_b} shows that the homogenized stiffness $\uuuu{C}^\hom$ has a qualitatively better $h$-convergence on transformed grids (near quadratic) than on uniform ones (linear). Using one 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$ ($250^3$ nodes).

\subsection{Discussion}

Numerical experiments in this section focus on a simple boundary-conforming grid adaptation strategy, the \emph{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.~\cite{Brisard2012,Schneider2015,Kabel2015,Vondrejc2016}. 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~\cite{Zecevic2026} 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 \emph{boundary node projection transform}.

In turn, the grid adaptation strategy based on optimal transport theory introduced in~\cite{Bellis2024} allows to control node density in regions where the solution exhibits high strain gradients. Numerical experiments in~\cite{Bellis2024} indicate that the spatial convergence can be significantly improved with this non boundary-conforming $r$-adaptation strategy. Yet, no $h$-convergence estimate of the homogenized properties is currently available. Indeed, the transformed grid is not boundary-conforming. Thus, a two-step grid adaptation may prove more beneficial to reach both $h^2$ convergence and additional reduction of discretization errors by $r$-adaptivity: 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, to benefit from FEM spatial convergence estimates.

On the down side, local node densification allowed by transforms considered in~\cite{Zecevic2026,Bellis2024} 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.

\section{Conclusion}

\paragraph{Main findings}
\begin{enumerate}
\item ``FFT-based methods'' with FEM discretization on structured grid can be efficiently formulated with isoparametric grid transformation.
\item The Green preconditioning is applied in the reference, or pulled-back, grid.
\item 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 2 to 8.
\item Boundary-conforming (or nearly) grids drastically improve the $h$-convergence of the homogenized property $\uuuu{C}^\hom$ w.r.t.\ uniform grids with or without composite voxels. For problems with sufficiently smooth boundary, a quadratic $h$-convergence of $\uuuu{C}^\hom$ is achieved for bi- or tri-linear $Q_1$ elements.
\item Boundary-conforming transforms significantly reduce errors on local strain fields. Hourglassing of bi- or tri-linear elements with reduced integration can be diminished.
\end{enumerate}

\paragraph{Perspectives}
\begin{enumerate}
\item 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~\cite{Antolin2022}.
\item 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~\cite{Thompson1998,Zecevic2026}, or optimal transport theory~\cite{Bellis2024} with subsequent projection of nodes on boundaries, remain to be investigated. The ideal transform should seek optimization of the accuracy while maintaining a favorable eigenvalue spectrum of the Green preconditioned operator.
\end{enumerate}

\section*{Acknowledgments}

The manuscript was written through contributions of all authors.

\section*{Declaration of interests}

The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.

\section*{Supplementary material}

Supporting information for this article, including example scripts, is available in~\cite{Bignonnet_2026_18458599} and on the journal's website under \printDOI\ or from the authors.

\CDRsupplementaryTwotypes{supplementary-material}{\cdrattach{linear-conductivity-2D.ipynb}}
\CDRsupplementaryTwotypes{supplementary-material}{\cdrattach{linear-elasticity-2D-hashin.ipynb}}

\appendix

\section{Conductivity}\label{sec:conductivity}

For thermal conductivity, the homogenization problem~\eqref{eq:homogenization-problem} is modified as follows. The scalar $u$ denotes the periodic fluctuation temperature and the vector $\uv{\sigma}$ the heat flux. The temperature gradient $\uv{\epsilon}$ is split in $\uv{\epsilon} = \uv{E} + \grad u$ where $\uv{E}$ is the macroscopic temperature gradient. The local conductivity is the second order tensor $\uu{C}$, and $\uv{\sigma} = - \uu{C} \cdot \uv{\epsilon}$. For grid transformation, the pulled-back conductivity is $\tilde{C}_{ij} = J F^{-1}_{ik} C_{k\ell} F^{-T}_{\ell j}$ in place of~\eqref{eq:transformed-stiffness}.

For Hashin's composite inclusion problem in Section~\ref{sec:hashin} and Appendix~\ref{sec:hashin-2d}, the three phases have isotropic conductivity $\uu{C}_i = C_i \uu{I}$. The equivalent conductivity of the inclusion is given by~\eqref{eq:hashin-bulk} up to the replacement of $\kappa$ by $C$ and modification $\alpha = (C_2 - C_1) \big/ \parens[\big]{(d-1)C_2 + C_1}$, $\beta=d$. The temperature gradient is still obtained by~\eqref{eq:hashin-bulk-field}, where $\uu{E}$ is now a vector.

\newpage

\section{Hashin's composite inclusion in dimension~2}\label{sec:hashin-2d}

\begin{figure}[!h]
\begin{subfigure}{0.4\textwidth}
\centering
\includegraphics[width=.84\textwidth]{field/hashin_conduction_2_gradux_rev.pdf}
\caption{Temperature gradient ${\epsilon_x}_{\num}$}
\end{subfigure}
\hspace{.05\textwidth}
\begin{subfigure}{0.4\textwidth}
\centering
\includegraphics[width=.84\textwidth]{field/hashin_conduction_2_error_gradux_rev.pdf}
\caption{Error on temperature gradient ${\epsilon_x}_{\num} - {\epsilon_x}_{\rmref}$}
\end{subfigure}
\caption{Hashin's composite inclusion in 2D conductivity with $\uv{E} = \uv{e}_x$, $C_1 = 0.01$ and $C_2=1$. Only one quarter or one half of fields is shown by symmetry.}\label{fig:hashin-field-conduction}
\end{figure}

\begin{figure}[!h]
\begin{subfigure}{0.48\textwidth}
\includegraphics[width=\textwidth]{convergence/error_hom_hashin_Conduction_2_rev.pdf}
\caption{$\abs{C^\hom_{\num}-C^\hom_{\rmref}}/C^\hom_{\rmref}$}
\end{subfigure}
\hfill
\begin{subfigure}{0.48\textwidth}
\includegraphics[width=\textwidth]{convergence/error_field_hashin_Conduction_2_rev.pdf}
\caption{RMSE on $\uv{\epsilon}$}
\end{subfigure}
\caption{Hashin's composite inclusion in 2D conductivity for $\uv{E} = \uv{e}_x$, $C_1 = 0.01$ and $C_2=1$.}\label{fig:hashin-convergence-conduction}
\end{figure}

\begin{figure}[!h]
\begin{subfigure}{0.48\textwidth}
\includegraphics[width=\textwidth]{field/hashin_elasticity_2_epsilonxx_rev.pdf}
\caption{Strain $\epsilon_{xx}$}
\end{subfigure}
\hfill
\begin{subfigure}{0.48\textwidth}
\includegraphics[width=\textwidth]{field/hashin_elasticity_2_error_epsilonxx_rev.pdf}
\caption{Error on strain ${\epsilon_{xx}}_{\num} - {\epsilon_{xx}}_{\rmref}$}
\end{subfigure}
\caption{Hashin's composite inclusion in 2D elasticity with $\uu{E} =\uu{I}_{2d}$. Only one quarter or one half of fields is shown by symmetry.}\label{fig:hashin-field-elasticity-2}
\end{figure}

\vskip \baselineskip

\begin{figure}[!h]
\begin{subfigure}{0.48\textwidth}
\includegraphics[width=\textwidth]{convergence/error_hom_hashin_Elasticity_2_rev.pdf}
\caption{$\abs[\big]{\kappa^\hom_{\num}-\kappa^\hom_{\rmref}}/\kappa^\hom_{\rmref}$}
\end{subfigure}
\hfill
\begin{subfigure}{0.48\textwidth}
\includegraphics[width=\textwidth]{convergence/error_field_hashin_Elasticity_2_rev.pdf}
\caption{RMSE on $\uu{\epsilon}$}
\end{subfigure}
\caption{Hashin's composite inclusion in 2D linear elasticity for $\uv{E} = \uu{I}_{2d}$.}\label{fig:hashin-convergence-elasticity-2}
\end{figure}

\begin{figure}[t]
\begin{subfigure}{0.48\textwidth}
\includegraphics[width=\textwidth]{convergence/convergence_PCG_hashin_Conduction_2_rev.pdf}
\caption{Conductivity}\label{fig:hashin-elasticity-iterations_a}
\end{subfigure}
\hfill
\begin{subfigure}{0.48\textwidth}
\includegraphics[width=\textwidth]{convergence/convergence_PCG_hashin_Elasticity_2_rev.pdf}
\caption{Elasticity}\label{fig:hashin-elasticity-iterations_b}
\end{subfigure}
\caption{PCG iteration count for Hashin's composite inclusion in 2D for \eqref{fig:hashin-elasticity-iterations_a}~conduction and \eqref{fig:hashin-elasticity-iterations_b}~linear elasticity problems.}\label{fig:hashin-elasticity-iterations}
\end{figure}

\printbibliography

\end{document}