%~Mouliné par MaN_auto v.0.40.4 (a17591c0) 2026-06-10 18:04:50
\documentclass[CRMECA,Unicode,XML,screen,biblatex]{cedram}
\usepackage{subfig}
\usepackage{bm}
\usepackage{upgreek}

\DeclareMathOperator{\Div}{div}
\DeclareMathOperator{\DFT}{DFT}
\DeclareMathOperator{\DFTt}{DFT3}
\DeclareMathOperator{\IDFT}{IDFT}
\DeclareMathOperator{\err}{err}



%~\COI{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.}
%~\COI{Les auteurs ne travaillent pas, ne conseillent pas, ne possèdent pas de parts, ne reçoivent pas de fonds d'une organisation qui pourrait tirer profit de cet article, et n'ont déclaré aucune autre affiliation que leurs organismes de recherche}

\addbibresource{CRMECA_Amadou-Sanoko_20260132.bib}

\newcommand{\fine}{\textrm{fine}}
\newcommand{\rnu}{\textrm{nu}}
\newcommand{\rnum}{\textrm{num}}

%Better widebar

\makeatletter
\let\save@mathaccent\mathaccent
\newcommand*\if@single[3]{%
  \setbox0\hbox{${\mathaccent"0362{#1}}^H$}%
  \setbox2\hbox{${\mathaccent"0362{\kern0pt#1}}^H$}%
  \ifdim\ht0=\ht2 #3\else #2\fi
  }
%The bar will be moved to the right by a half of \macc@kerna, which is computed by amsmath:
\newcommand*\rel@kern[1]{\kern#1\dimexpr\macc@kerna}
%If there's a superscript following the bar, then no negative kern may follow the bar;
%an additional {} makes sure that the superscript is high enough in this case:
\newcommand*\widebar[1]{\@ifnextchar^{{\wide@bar{#1}{0}}}{\wide@bar{#1}{1}}}
%Use a separate algorithm for single symbols:
\newcommand*\wide@bar[2]{\if@single{#1}{\wide@bar@{#1}{#2}{1}}{\wide@bar@{#1}{#2}{2}}}
\newcommand*\wide@bar@[3]{%
  \begingroup
  \def\mathaccent##1##2{%
%Enable nesting of accents:
    \let\mathaccent\save@mathaccent
%If there's more than a single symbol, use the first character instead (see below):
    \if#32 \let\macc@nucleus\first@char \fi
%Determine the italic correction:
    \setbox\z@\hbox{$\macc@style{\macc@nucleus}_{}$}%
    \setbox\tw@\hbox{$\macc@style{\macc@nucleus}{}_{}$}%
    \dimen@\wd\tw@
    \advance\dimen@-\wd\z@
%Now \dimen@ is the italic correction of the symbol.
    \divide\dimen@ 3
    \@tempdima\wd\tw@
    \advance\@tempdima-\scriptspace
%Now \@tempdima is the width of the symbol.
    \divide\@tempdima 10
    \advance\dimen@-\@tempdima
%Now \dimen@ = (italic correction / 3) - (Breite / 10)
    \ifdim\dimen@>\z@ \dimen@0pt\fi
%The bar will be shortened in the case \dimen@<0 !
    \rel@kern{0.6}\kern-\dimen@
    \if#31
      \overline{\rel@kern{-0.6}\kern\dimen@\macc@nucleus\rel@kern{0.4}\kern\dimen@}%
      \advance\dimen@0.4\dimexpr\macc@kerna
%Place the combined final kern (-\dimen@) if it is >0 or if a superscript follows:
      \let\final@kern#2%
      \ifdim\dimen@<\z@ \let\final@kern1\fi
      \if\final@kern1 \kern-\dimen@\fi
    \else
      \overline{\rel@kern{-0.6}\kern\dimen@#1}%
    \fi
  }%
  \macc@depth\@ne
  \let\math@bgroup\@empty \let\math@egroup\macc@set@skewchar
  \mathsurround\z@ \frozen@everymath{\mathgroup\macc@group\relax}%
  \macc@set@skewchar\relax
  \let\mathaccentV\macc@nested@a
%The following initialises \macc@kerna and calls \mathaccent:
  \if#31
    \macc@nested@a\relax111{#1}%
  \else
%If the argument consists of more than one symbol, and if the first token is
%a letter, use that letter for the computations:
    \def\gobble@till@marker##1\endmarker{}%
    \futurelet\first@char\gobble@till@marker#1\endmarker
    \ifcat\noexpand\first@char A\else
      \def\first@char{}%
    \fi
    \macc@nested@a\relax111{\first@char}%
  \fi
  \endgroup
}
\makeatother

\let\oldbar\bar
\renewcommand*{\bar}[1]{{\mathchoice{\widebar{#1}}{\widebar{#1}}{\widebar{#1}}{\oldbar{#1}}}}

\makeatletter
\def\@setafterauthor{%
  \vglue3mm%
%  \hspace*{0pt}%
\begingroup\hsize=12cm\advance\hsize\abstractmarginL\raggedright
\noindent
%\hspace*{\abstractmarginL}\begin{minipage}[t]{10cm}
   \leftskip\abstractmarginL
  \normalfont\Small
  \@afterauthor\par
\endgroup
\vskip2pt plus 3pt minus 1pt
}
\makeatother

\newcommand{\moy}[2][]{{\left\langle{#2}\right\rangle}_{#1}}

\newcommand\VbE{{\bigl\lVert\overline{\mathbf{E}}\bigr\rVert}}
\newcommand{\bs}{\,\mk\overline{\!\mk s\mK}}

\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}

\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{Reformulation and influence of discretization in Fourier-series-based FFT methods for heterogeneous materials}
\alttitle{Reformulation et influence de la discrétisation dans les méthodes FFT basées sur les séries de Fourier pour les matériaux hétérogènes}

\author{\firstname{Abdoul Magid} \lastname{Amadou Sanoko}}
\address{Univ. Bordeaux, CNRS, Bordeaux INP, I2M, UMR 5295, F-33400 Talence, France}
\address{Arts et Metiers Institute of Technology, CNRS, Bordeaux INP, I2M, UMR 5295, F-33400 Talence, France}
\email[A.M. Amadou Sanoko] {abdoul-magid.amadou-sanoko@u-bordeaux.fr}

\author{\firstname{Simon} \lastname{Essongue}}
\address[1]{Univ. Bordeaux, CNRS, Bordeaux INP, I2M, UMR 5295, F-33400 Talence, France}
\address[2]{Arts et Metiers Institute of Technology, CNRS, Bordeaux INP, I2M, UMR 5295, F-33400 Talence, France}
\email[S. Essongue] {simon.essongue-boussougou@u-bordeaux.fr}

\author{\firstname{Joseph} \lastname{Paux}}
\address{Université Sorbonne Paris Nord, Laboratoire des Sciences des Procédés et des Matériaux, LSPM, CNRS, UPR 3407, F-93430, Villetaneuse, France}
\email[J. Paux] {joseph.paux@univ-paris13.fr}

\author{\firstname{Léo} \lastname{Morin}\IsCorresp}
\address[1]{Univ. Bordeaux, CNRS, Bordeaux INP, I2M, UMR 5295, F-33400 Talence, France}
\address[2]{Arts et Metiers Institute of Technology, CNRS, Bordeaux INP, I2M, UMR 5295, F-33400 Talence, France}
\email[L. Morin] {leo.morin@u-bordeaux.fr}

%~\ESM{Supplementary material for this article is supplied as a separate archive available from the journal's website under article's URL or from the author.}

\keywords{\kwd{FFT-based solvers}
\kwd{Strong and weak formulations}
\kwd{Fourier series}
\kwd{Green operators}}

\altkeywords{\kwd{Solveurs basés sur la FFT}
\kwd{formulations forte et faible}
\kwd{séries de Fourier}
\kwd{opérateurs de Green}}

\COI{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.}

\begin{abstract}
The aim of this work is to revisit the formulation of FFT-based methods for heterogeneous materials in the periodic setting. These numerical methods are based on the iterative resolution of an auxiliary problem involving a reference homogeneous material and a polarization tensor. Assuming a description of the fields by Fourier series, we show the equivalence between three discretization approaches based either on the strong and the weak formulation of the problem.~A special emphasis is put on the representation of the local fields described by Fourier series including the possibility of using non-uniform grid. Numerical experiments are performed on a model problem of conductivity with a checkerboard microstructure for which an analytical solution allows to assess the effect of Fourier modes together with the grid discretization. The occurrence of oscillations is finally addressed by studying (generalized) discrete Green operators (still in the context of Fourier series) and interface spreading approaches based on smoothing techniques.
\end{abstract}

\begin{altabstract}
L’objectif de ce travail est de revisiter la formulation des méthodes basées sur la FFT pour les matériaux hétérogènes dans un cadre périodique. Ces méthodes numériques reposent sur la résolution itérative d’un problème auxiliaire impliquant un matériau homogène de référence et un tenseur de polarisation. En supposant une description des champs par des séries de Fourier, nous montrons l’équivalence entre trois approches de discrétisation fondées soit sur la formulation forte, soit sur la formulation faible du problème. Une attention particulière est portée à la représentation des champs locaux décrits par des séries de Fourier, incluant la possibilité d’utiliser un maillage non uniforme. Des expériences numériques sont réalisées sur un problème modèle de conductivité avec une microstructure en damier, pour lequel une solution analytique permet d’évaluer l’effet des modes de Fourier conjointement à la discrétisation du maillage. L’apparition d’oscillations est enfin étudiée en s’appuyant sur des opérateurs de Green discrets (généralisés) toujours dans le cadre des séries de Fourier, ainsi que sur des approches d’étalement d’interfaces basées sur des techniques de lissage.
\end{altabstract}


\begin{document}
\maketitle

\section{Introduction} FFT-based methods for heterogeneous materials, introduced in the seminal article of Moulinec and Suquet (1998) \cite{moulinec_numerical_1998}, are a class of efficient methods that have been developed in the context of periodic homogenization of mechanical properties in heterogeneous microstructures, mainly for \emph{micromechanics} (see the review papers of~\cite{schneider_review_2021} and~\cite{lucarini_fft_2021}). This method is based on an iterative solution of the so-called Lippmann--Schwinger equation using fast Fourier transforms (FFT). The main {characteristics of these methods} are that (i) calculations are based on the computational complexity of fast Fourier transforms in $O(N \log N)$, (ii) it is image-based (obtained from {imaging} techniques such as scanning electron microscopy (SEM) or tomography) so it does not require meshing operations and (iii) it can be suited for parallel implementation (see e.g.~\cite{amitex_2022}). It is generally used to calculate the local and overall fields, and has been applied to various problems in heterogeneous microstructures including notably~\cite{moulinec_numerical_1998}, J2-plasticity~\cite{moulinec_numerical_1998}, crystal (elasto)viscoplasticity~\cite{lebensohn_n-site_2001}, dislocation-mediated plasticity~\cite{brenner_numerical_2014,bertin_fft-based_2015}, conductivity~\cite{eyre_fast_1999}, piezoelectricity~\cite{brenner_numerical_2009}, porous ductile solids~\cite{bilger_effect_2005,paux_plastic_2018}.

The original method developed by Moulinec and Suquet (1998) \cite{moulinec_numerical_1998} is implicitly based on the resolution of the strong formulation (of the equations of elastostatics) using infinite Fourier series. It was then shown by~\cite{vondrejc_fft-based_2014} the formulation of~\cite{moulinec_numerical_1998} is equivalent to a Galerkin method using trigonometric polynomials as interpolation basis functions. The use of a Galerkin method based on a weak formulation instead of a strong formulation is particularly important in extensions of FFT-based methods to non-periodic boundary conditions using sine and cosine series
%for which the auxiliary problem can only be solved using a Galerkin method
\cite{risthaus2024imposing,morin_fast_2024-1,risthaus2024fft,paux_discrete_2025-1,paux_discrete_2025}.

The scheme of~\cite{moulinec_numerical_1998} relies on the introduction of a Green operator which can be expressed easily in Fourier space using Fourier transforms and standard spectral derivative rule, defining the so-called \emph{continuous} Green operator; the terminology ``continuous'' has been notably introduced by~\cite{willot_fourier-based_2014}, to highlight its differences with other Green operators constructed using different derivative rules and/or different formulations. In the work of~\cite{willot_fourier-based_2014}, a finite difference scheme is used to approximate the partial derivatives defining the local elliptic problem (see also~\cite{willot2015fourier}) and leads to the so-called \emph{discrete} Green operator that was shown to improve both the convergence of the iterative procedure and the local fields. It was also shown in several works that the use of uniform finite elements on periodic Cartesian grids leads to similar discrete Green operators depending on the type of elements considered~\cite{zeman2017finite,schneider_fft-based_2017}: the FFT solver can then be seen as a preconditioner of the linear system arising in the FEM discretization.

The aim of this work is to revisit several aspects of the method of~\cite{moulinec_numerical_1998}. First, we compare several techniques of the numerical resolution of the auxiliary problem, either based on strong or weak formulations, by describing explicitly the unknown fields using (infinite or partial) Fourier series following initial ideas of~\cite{vondrejc_fft-based_2014}; this allows to clearly show the contribution of the number of modes, aliasing and underintegration in the procedure of resolution {which are classical topics in the spectral methods literature~\cite{canuto_spectral_2006, trefethen2000spectral}}. As the fields are described by Fourier series, a particular emphasis is put on the representation of solution fields between grid nodes with the possibility of using non-uniform sampling: the occurrence of spurious oscillations is carefully assessed especially between grid nodes. To address the possible remedies for oscillations, we derive, still in the context of Fourier series, generalized Green operators that are based on modified derivative laws which allows to retrieve the so-called continuous and discrete operators~\cite{willot_fourier-based_2014}. The importance of ``spreading'' (i.e.\/ smoothing material properties) \cite{morin_periodic_2021} is also assessed in order to remove these spurious oscillations.

The paper is organized as follows. In Section~\ref{sec:prelim}, the framework of Fourier series and transforms is recalled as well as the physical problem of conductivity considered throughout the article. { The numerical discretization of the auxiliary problem is then done in Section~\ref{sec:formulations}, based either on strong and weak formulations together with numerical approximations.} The representation of solution fields on an arbitrary grid is then addressed in Section~\ref{sec:representation}. Numerical experiments are performed in Section~\ref{sec:num} on a problem of conductivity with a checkerboard microstructure. Finally, remedies for oscillations, i.e.\/ discrete Green operators and smoothing techniques, are discussed in Section~\ref{sec:disc}.


\section{Preliminaries}\label{sec:prelim}
\subsection{Equations of conductivity} We consider a problem of conductivity (as a model problem for FFT-based solvers) which consists in the computation of the electric field $\mathbf{E}(\varphi(\mathbf{x}))$, current $\mathbf{J}(\mathbf{x})$ and electric potential $\varphi(\mathbf{x})$, at each point $\mathbf{x}$ in $\Omega$, for a given conductivity field $ \mathbf{c}(\mathbf{x})$. The unit-cell is a square or cubic domain $\Omega = [-L/2,~L/2]^d$ ($d=2$ or $3$) with $L$ the period of the microstructure. Tensorial components refer to a system of Cartesian coordinates $(\mathbf{e}_1;\mathbf{e}_2)$ in 2D and $(\mathbf{e}_1;\mathbf{e}_2;\mathbf{e}_3)$ in 3D. Periodic boundary conditions are considered and the local electric field $\mathbf{E}$ is split into its average $\bar{\mathbf{E}}=\moy{\mathbf{E}}$ (where $\moy{\cdot}$ denotes the spatial average over $\Omega$) and a fluctuation term $\mathbf{E}(\varphi^*(\mathbf{x}))$:
\begin{equation}
\mathbf{E}(\varphi(\mathbf{x})) = \mathbf{E}(\varphi^*(\mathbf{x})) + \bar{\mathbf{E}} \quad \quad \text{or equivalently}\quad \quad \varphi(\mathbf{x})=\varphi^*(\mathbf{x})-\bar{\mathbf{E}}\cdot\mathbf{x}.
\end{equation}
Periodic boundary conditions imply that the fluctuating part $\varphi^*$ is periodic (notation: $\varphi^*\#$) and that the term $\mathbf{J}\cdot\mathbf{n}$ is anti-periodic on the boundary between two neighboring cells with $\mathbf{n}$ the outer normal along the boundary $\partial \Omega$ of $\Omega$ (notation: $\mathbf{J}\cdot\mathbf{n} - \#$).

The equations of conductivity are
\begin{align}\label{eq:conduct_hetero}
\left\{
\begin{array}{lll}
\Div \mathbf{J}(\mathbf{x}) & = & {0} \\
\mathbf{J}(\mathbf{x}) & = & \mathbf{c}(\mathbf{x}) : \left(\mathbf{E}(\varphi^*(\mathbf{x})) + \bar{\mathbf{E}} \right) \\
\mathbf{E}\left(\varphi^{*}(\mathbf{x})\right) & = & - \mathbf{\nabla}\varphi^{*}(\mathbf{x}),
\end{array}\right.
\end{align}
where $\mathbf{c}(\mathbf{x})$ is the local second-order conductivity tensor. In the particular case of isotropic conductivity, it reads
\begin{equation}
\mathbf{c}(\mathbf{x}) = c(\mathbf{x})\mathbf{I}_d
\end{equation}
with $c(\mathbf{x})$ the local (scalar) conductivity field and $\mathbf{I}_d$ the second-order identity tensor. Hereafter we will consider isotropic constituents for simplicity (although the method is designed for arbitrary anisotropic behaviors).

\subsection{Fourier series and the discrete Fourier transform} We recall important results related to Fourier series and the associated discrete Fourier transforms.

\paragraph*{Fourier series: the 1-d case.} We consider a periodic function $f$ with period $L$. Such function can be described by the Fourier series
\begin{equation}
f(x) = \sum_{i=-\infty}^{+\infty} F_i \exp\left(\imath \xi_i x\right),\quad \quad \quad \forall x\in[0,L],
\end{equation}
where $\imath$ is the imaginary unit and $\xi_i$ is some ``frequency'' parameter that is given by
\begin{equation}\label{eq:xi}
\xi_i = 2\pi \frac{i}{L},\quad \quad \quad \quad \quad i=-\infty,\dots,+\infty,
\end{equation}
and $F_i$ are the Fourier coefficients given by
\begin{equation}\label{eq:coeff_fourier_series}
F_i = \frac{1}{L}\int_0^L f(x) \exp\left(-\imath \xi_i x\right) \mathrm{d}x.
\end{equation}
It is interesting to note that the (first) derivative of $f$ is the Fourier series
\begin{equation}
f'(x) = \sum_{i=-\infty}^{+\infty} \imath \xi_i F_i \exp\left(\imath \xi_i x\right),\quad \quad \quad \forall x\in[0,L].
\end{equation}

In addition, the associated partial (i.e.\/ truncated) series of order $M+1$ (for every positive even integer $M$) is defined by
\begin{equation}
f(x) = \sum_{i=-M/2}^{M/2} F_i \exp\left(\imath \xi_i x\right),\quad \quad \quad \forall x\in[0,L].
\end{equation}

\paragraph*{Discrete Fourier transform: the 1d-case.}

The computation of the Fourier coefficients $F_i$ can be done approximately using the discrete Fourier transform (DFT). An efficient calculation of the DFT is classically performed using fast Fourier transforms (FFT) \cite{frigo_fftw_1998}.

The 1-d domain of size $[0,~L]$ is discretized uniformly with $N+1$ segments so the spatial scale associated with the uniform grid is $\Delta x = L/(N+1)$ and we denote by $x_a=a \Delta x$ (for $a=0,\dots,N$) the grid points. The grid points values of function $f$ are denoted by $f_a=f(x_a)$ (for $a=0,\dots,N$) and may thus be written as a vector $\bm{f}$ of size $N+1$ as
\begin{equation}
\bm{f} =
\begin{pmatrix} f_0 & f_1 & \dots & f_{N-1} & f_N
\end{pmatrix}
\end{equation}

The discrete Fourier transform of $\bm{f}$, denoted by $\DFT(\bm{f})$, is the vector $\widehat{\bm{F}}$
\begin{equation}\label{eq:DFT}
\widehat{\bm{F}}= \DFT(\bm{f}) =
\begin{pmatrix} \widehat{F}_{-\frac{N}{2}} & \widehat{F}_{-\frac{N}{2}+1} & \dots \widehat{F}_{\frac{N}{2}-1} & \widehat{F}_{\frac{N}{2}}
\end{pmatrix},
\end{equation}
where the $\widehat{F}_i$ are defined by
\begin{equation}
\widehat{F}_i = \sum_{j=0}^{N} f_j \exp\left(-\frac{2\pi \imath }{N}ij \right) = \sum_{j=0}^{N} f_j \exp\left(-\imath \xi_i x_j \right),\quad \quad \quad i=-N/2,\dots,N/2.
\end{equation}
The inverse discrete Fourier transform, denoted by $\bm{f}=\IDFT(\widehat{\bm{F}})$, is given by
\begin{equation}
f_i = \frac{1}{N} \sum_{j=-N/2}^{N/2} \widehat{F}_j \exp\left(\frac{2\pi \imath }{N}ij \right).
\end{equation}

The coefficients $\widehat{F}_i$ provide a numerical approximation of the coefficients $F_i$ of the Fourier series given by equation~\eqref{eq:coeff_fourier_series}, using a standard trapezoidal rule~\cite{trefethen2014exponentially}
\begin{equation}\label{eq:dft_fourier_series}
F_i \simeq \frac{1}{N} \widehat{F}_i,\quad \quad \quad i=-N/2,\dots,N/2.
\end{equation}


\paragraph*{Extension to the 3-d case.} The extension of Fourier series and transforms to the 3-d case is straightforward. We consider a periodic function $f$ in the prismatic domain $\Omega=[0,L_x]\times[0,L_y]\times[0,L_z]$; its Fourier series simply reads
\begin{equation}
f(x,y,z)=\sum_{i=-\infty}^{+\infty}\sum_{j=-\infty}^{+\infty}\sum_{k=-\infty}^{+\infty} F_{ijk}\exp\left(\imath \xi_i^x x\right)\exp\left(\imath \xi_j^y y\right)\exp\left(\imath \xi_k^z z\right),\quad\forall(x,y,z)\in\Omega,
\end{equation}
where the frequencies $\xi_i^x$, $\xi_j^y$ and $\xi_k^z$ are given by
\begin{equation}
\xi_i^x = 2\pi \frac{i}{L_x},\quad\quad\quad \xi_j^y = 2\pi \frac{j}{L_y}, \quad\quad\quad\xi_k^z = 2\pi \frac{k}{L_z},
\end{equation}
with $i=-\infty,\dots,+\infty$, $j=-\infty,\dots,+\infty$ and $k=-\infty,\dots,+\infty$. The Fourier coefficients $F_{ijk}$ are given~by
\begin{equation}\label{eq:coeff_fourier_series_3d}
{F_{ijk}} = \frac{1}{L_x L_y L_z}\int_\Omega f(x,y,z) \exp\left(-\imath \xi_i^x x\right)\exp\left(-\imath \xi_j^y y\right)\exp\left(-\imath \xi_k^z z\right)\mathrm{d}x\mathrm{d}y\mathrm{d}z.
\end{equation}
The coefficients $F_{ijk}$ can be approximated using the 3-d discrete Fourier transform. The 3-d domain $\Omega$ is discretized uniformly with $(N_x+1)\times(N_y+1)\times(N_z+1)$ voxels so that the spatial scales are $\Delta x = L_x/(N_x+1)$, $\Delta y = L_y/(N_y+1)$ and $\Delta z = L_z/(N_z+1)$. The coordinates of the grid points are denoted by $x_a=a \Delta x$ (for $a=0,\dots,N_x$), $y_b=b \Delta y$ (for $b=0,\dots,N_y$) and $z_c=c \Delta z$ (for $c=0,\dots,N_z$). The values of the function $f$ at the grid points are denoted by $f_{abc}=f(x_a,y_b,z_c)$ (for $a=0,\dots,N_x$, $b=0,\dots,N_y$ and $c=0,\dots,N_z$) and we denote by $\bm{f}$ the array of size $(N_x+1)\times (N_y+1)\times (N_z+1)$ containing the values of function $f$ at the grid points. In 3-d, the discrete Fourier transform of $\bm{f}$, denoted by $\DFTt(\bm{f})$, is the $(N_x+1)\times (N_y+1)\times (N_z+1)$ array $\widehat{\bm{F}}$
\begin{equation}
\widehat{\bm{F}}= \DFTt(\bm{f}),
\end{equation}
whose components are defined by
\begin{equation}\label{eq:dft_3d}
\widehat{F}_{ijk} = \sum_{l=0}^{N_x}\sum_{m=0}^{N_y}\sum_{n=0}^{N_z} f_{lmn} \exp\left(-\frac{2\pi \imath }{N_x}il \right) \exp\left(-\frac{2\pi \imath }{N_y}jm \right)\exp\left(-\frac{2\pi \imath }{N_z}kn \right).
\end{equation}
As in the 1-d case, the coefficients $\widehat{F}_{ijk}$ provide a numerical approximation of the coefficients ${F}_{ijk}$ using the trapezoidal rule:
\begin{equation}\label{eq:series_transform_equivalence}
{F}_{ijk} \simeq \frac{1}{N_x N_y N_z} \widehat{F}_{ijk},\quad \;\; i=-N_x/2,\dots,N_x/2, \quad j=-N_y/2,\dots,N_y/2, \quad k=-N_z/2,\dots,N_z/2.
\end{equation}


\section{{ Discretization approaches}}\label{sec:formulations}
In the seminal paper of~\cite{moulinec_numerical_1998}, the type of formulation (strong or weak) and the associated numerical discretization leading to the numerical algorithm was not explicitly detailed, although it is generally considered that~\cite{moulinec_numerical_1998}'s scheme was obtained using a strong formulation (see e.g.~\cite{lucarini2019dbfft}). { In this section, we discuss the numerical resolution of the auxiliary problem from different approaches (strong and weak formulations) and discretization methods (truncation and DFT approximation), following ideas of~\cite{vondrejc_fft-based_2014}.}

\subsection{The auxiliary problem} First, the problem of heterogeneous conductivity~\eqref{eq:conduct_hetero} is written as the auxiliary problem,
\begin{align}\label{eq:aux}
\left\{
\begin{array}{lll}
\Div \mathbf{J}(\mathbf{x}) & = & {0} \\
\mathbf{J}(\mathbf{x}) & = & c_0 \mathbf{E}(\varphi^*(\mathbf{x})) + \bm{\uptau}(\mathbf{x}) \\
\mathbf{E}\left(\varphi^{*}(\mathbf{x})\right) & = & - \mathbf{\nabla}\varphi^{*}(\mathbf{x}),
\end{array}\right.
\end{align}
where $c_0$ is the conductivity of some homogeneous reference medium and $\bm{\uptau}$ is the polarization tensor defined by,
\begin{equation}\label{eq:polarization}
\bm{\uptau}(\bm{x})= \left(\mathbf{c}(\mathbf{x}) - c_0 \mathbf{I}_d\right) : \left(\mathbf{E}(\varphi^*(\mathbf{x}))+\bar{\mathbf{E}} \right) + c_0\bar{\mathbf{E}}.
\end{equation}

The principle of resolution is to determine the fluctuation $\varphi^{*}$ solution of the problem~\eqref{eq:aux}, assuming that $\bm{\uptau}$ is known, and then iterate to find $\bm{\uptau}$ corresponding to the solution of the problem~\eqref{eq:conduct_hetero}. From the knowledge of $\bm{\uptau}$, the numerical determination of $\varphi^{*}$ solution of the problem~\eqref{eq:aux} can then be done using several discretization strategies.

\subsection{{Discretization of the strong formulation by DFT}} The strong formulation was implicitly used in the seminal article of~\cite{moulinec_numerical_1998} by assuming that all fields are described by \emph{(full) Fourier series}; the conductivity field and the fluctuation field then read
\begin{align}\label{eq:infinite_strong}
\varphi^{*}(x,y,z)&=\sum_{i=-\infty}^{+\infty}\sum_{j=-\infty}^{+\infty}\sum_{k=-\infty}^{+\infty} \Phi^*_{ijk}\exp\left(\imath \xi_i^x x\right)\exp\left(\imath \xi_j^y y\right)\exp\left(\imath \xi_k^z z\right), \\
c(x,y,z)&=\sum_{i=-\infty}^{+\infty}\sum_{j=-\infty}^{+\infty}\sum_{k=-\infty}^{+\infty} C_{ijk}\exp\left(\imath \xi_i^x x\right)\exp\left(\imath \xi_j^y y\right)\exp\left(\imath \xi_k^z z\right).
\end{align}
As a consequence, the polarization field is also described by (full) Fourier series and can be written as
\begin{align}
\tau_x(x,y,z)&=\sum_{i=-\infty}^{+\infty}\sum_{j=-\infty}^{+\infty}\sum_{k=-\infty}^{+\infty} T^x_{ijk}\exp\left(\imath \xi_i^x x\right)\exp\left(\imath \xi_j^y y\right)\exp\left(\imath \xi_k^z z\right), \\
\tau_y(x,y,z)&=\sum_{i=-\infty}^{+\infty}\sum_{j=-\infty}^{+\infty}\sum_{k=-\infty}^{+\infty} T^y_{ijk}\exp\left(\imath \xi_i^x x\right)\exp\left(\imath \xi_j^y y\right)\exp\left(\imath \xi_k^z z\right), \\
\tau_z(x,y,z)&=\sum_{i=-\infty}^{+\infty}\sum_{j=-\infty}^{+\infty}\sum_{k=-\infty}^{+\infty} T^z_{ijk}\exp\left(\imath \xi_i^x x\right)\exp\left(\imath \xi_j^y y\right)\exp\left(\imath \xi_k^z z\right).
\end{align}
With this description of functions, the auxiliary problem reduces to the equation
\begin{multline}
c_0 \sum_{i=-\infty}^{+\infty}\sum_{j=-\infty}^{+\infty}\sum_{k=-\infty}^{+\infty} \left((\xi_i^x)^2 + (\xi_j^y)^2 + (\xi_k^z)^2 \right) \Phi^*_{ijk}\exp\left(\imath \xi_i^x x\right)\exp\left(\imath \xi_j^y y\right)\exp\left(\imath \xi_k^z z\right) \\
= - \sum_{i=-\infty}^{+\infty}\sum_{j=-\infty}^{+\infty}\sum_{k=-\infty}^{+\infty} \imath \left(\xi_i^x T^x_{ijk} + \xi_j^y T^y_{ijk} + \xi_k^z T^z_{ijk} \right)\exp\left(\imath \xi_i^x x\right)\exp\left(\imath \xi_j^y y\right)\exp\left(\imath \xi_k^z z\right),
\end{multline}
which leads to the relation between the Fourier series coefficients of the fluctuation and the components of the polarization
\begin{equation}\label{eq:strong_fourier_series_coeff}
\Phi^*_{ijk} = - \imath \frac{\xi_i^x T^x_{ijk} + \xi_j^y T^y_{ijk} + \xi_k^z T^z_{ijk} }{c_0\left((\xi_i^x)^2 + (\xi_j^y)^2 + (\xi_k^z)^2\right)}.
\end{equation}
It is important to note that this relation holds for all frequencies, i.e.\/ for $i=-\infty,\dots,+\infty$, $j=-\infty,\dots,+\infty$ and $k=-\infty,\dots,+\infty$ (except in the particular case $i=j=k=0$ where one has ${\Phi^*}_{000}=0$), {and is expressed in terms of the \emph{Fourier series coefficients of $\bm{\uptau}$}. \eqref{eq:strong_fourier_series_coeff} provides the analytical expression of the Fourier coefficients of $\varphi^\ast$, \emph{which need to be evaluated numerically}, leading to a first numerical approximation. Moreover, the numerical computation of the Fourier coefficients $\Phi^*_{ijk}$ from~\eqref{eq:strong_fourier_series_coeff} is only performed for a finite number of indexes $ijk$, leading to a second numerical approximation in the estimation of $\varphi^\ast$ due to the \emph{truncation} of its Fourier series into a partial one.}

In practical applications, the Fourier series coefficients are approximated by the coefficients of the discrete Fourier transform as shown by equation~\eqref{eq:series_transform_equivalence}, which need to be calculated using some discretization of the domain. The prismatic $\Omega$ is therefore discretized uniformly with $(N_x+1)\times(N_y+1)\times(N_z+1)$ voxels so that the spatial scales are $\Delta x = L_x/(N_x+1)$, $\Delta y = L_y/(N_y+1)$ and $\Delta z = L_z/(N_z+1)$. Following equation~\eqref{eq:series_transform_equivalence}, the discrete Fourier transforms of the fields $\varphi^{*}$, $\tau_x$, $\tau_y$ and $\tau_z$ provide a numerical approximation of the Fourier series coefficients:
\begin{equation}\label{eq:phi_T_DFT}
\!\Phi^*_{ijk} \simeq \frac{1}{N_x N_y N_z} \widehat{\Phi^*}_{ijk}, \;\;\: T^x_{ijk} \simeq \frac{1}{N_x N_y N_z} \widehat{T}^x_{ijk}, \;\;\: T^y_{ijk} \simeq \frac{1}{N_x N_y N_z} \widehat{T}^y_{ijk}, \;\;\: T^z_{ijk} \simeq \frac{1}{N_x N_y N_z} \widehat{T}^z_{ijk},\!
\end{equation}
for $i=-N_x/2,\dots,N_x/2$, $j=-N_y/2,\dots,N_y/2$ and $k=-N_z/2,\dots,N_z/2$. This finally leads to
\begin{equation}\label{eq:strong_dft}
\widehat{\Phi^*}_{ijk} = - \imath \frac{\xi_i^x \widehat{T}^x_{ijk} + \xi_j^y \widehat{T}^y_{ijk} + \xi_k^z \widehat{T}^z_{ijk} }{c_0\left((\xi_i^x)^2 + (\xi_j^y)^2 + (\xi_k^z)^2\right)},
\end{equation}
again for $i=-N_x/2,\dots,N_x/2$, $j=-N_y/2,\dots,N_y/2$ and $k=-N_z/2,\dots,N_z/2$, except in the particular case $i=j=k=0$ where one has $\widehat{\Phi^*}_{000}=0$.

It is worth noting that, if the FFT algorithm suggests the computations of the Fourier coefficients for these indexes, it is possible to take fewer Fourier modes, i.e.\/ compute~\eqref{eq:strong_dft} for $i=-M_x/2,\dots,M_x/2$, $j=-M_y/2,\dots,M_y/2$ and $k=-M_z/2,\dots,M_z/2$ with $M_x\leq N_x$, $M_y\leq N_y$ and $M_z\leq N_z$. While this choice seems to discard easy-to-compute coefficients and to reduce the order of the Fourier truncation of $\varphi^\ast$, therefore downgrades its numerical estimation, one must keep in mind that the Fourier coefficients of $\bm{\uptau}$ in~\eqref{eq:phi_T_DFT} are numerically estimated by DFT, therefore they have a numerical error which increases with $i$, $j$ and $k$ (see, for instance, \cite{tasche2001worst}). Considering $M_x\leq N_x$, $M_y\leq N_y$ and $M_z\leq N_z$, therefore suppressing \emph{high frequency} coefficients, discards the coefficients with highest numerical error coming from the DFT of $\bm{\uptau}$. The influence of this approximation will be evaluated in the next section.

Equation~\eqref{eq:strong_dft} provides the coefficients of the discrete Fourier transform of $\varphi^{*}$ as a function of the coefficients of the discrete Fourier transform of $\tau_x$, $\tau_y$ and $\tau_z$ and thus it implicitly defines the so-called Green operator on the electric potential field. Owing to the relation between the electric potential field and the electric field, this leads to the relations
\begin{align}\label{eq:strong_E}
\left\{
\begin{array}{lll}
\widehat{E}^x_{ijk} & = & \displaystyle - \xi_i^x \frac{\xi_i^x \widehat{T}^x_{ijk} + \xi_j^y \widehat{T}^y_{ijk} + \xi_k^z \widehat{T}^z_{ijk} }{c_0\left((\xi_i^x)^2 + (\xi_j^y)^2 + (\xi_k^z)^2\right)} \\[0.5cm]
\widehat{E}^y_{ijk} & = & \displaystyle - \xi_j^y \frac{\xi_i^x \widehat{T}^x_{ijk} + \xi_j^y \widehat{T}^y_{ijk} + \xi_k^z \widehat{T}^z_{ijk} }{c_0\left((\xi_i^x)^2 + (\xi_j^y)^2 + (\xi_k^z)^2\right)} \\[0.5cm]
\widehat{E}^z_{ijk} & = & \displaystyle - \xi_k^z \frac{\xi_i^x \widehat{T}^x_{ijk} + \xi_j^y \widehat{T}^y_{ijk} + \xi_k^z \widehat{T}^z_{ijk} }{c_0\left((\xi_i^x)^2 + (\xi_j^y)^2 + (\xi_k^z)^2\right)},
\end{array}\right.
\end{align}
for $i=-N_x/2,\dots,N_x/2$, $j=-N_y/2,\dots,N_y/2$ and $k=-N_z/2,\dots,N_z/2$ (except in the case $i=j=k=0$). In the case $i=j=k=0$, one has
\begin{equation}\label{eq:0-freq}
\widehat{E}^x_{000} = \bar{E}_x,\quad \widehat{E}^y_{000} = \bar{E}_y, \widehat{E}^z_{000} = \bar{E}_z,
\end{equation}
where $\bar{E}_x$, $\bar{E}_y$ and $\bar{E}_z$ are the components of the average of the electric field $\bar{\mathbf{E}}$. In equation~\eqref{eq:strong_E}, $\widehat{E}^x_{ijk}$, $\widehat{E}^x_{ijk}$ and $\widehat{E}^x_{ijk}$ are the coefficients of the discrete Fourier transforms of the components of the local electric field $\mathbf{E}$:
\begin{equation}\label{eq:DFT_E}
\widehat{\bm{E}}^x= \DFTt(\bm{E}_x),\quad \widehat{\bm{E}}^y= \DFTt(\bm{E}_y),\quad \widehat{\bm{E}}^z= \DFTt(\bm{E}_z).
\end{equation}
Equation~\eqref{eq:strong_E} can be alternatively written
\begin{equation}
\widehat{\bm{E}}(\bm{\xi}) = - \widehat{\bm{\Gamma}}^0(\bm{\xi}):\widehat{\bm{\tau}}(\bm{\xi})\quad \forall \bm{\xi}\neq\bm{0},\quad \widehat{\bm{E}}(\bm{0})=\bar{\mathbf{E}},
\end{equation}
where $\bm{\xi}$ denotes the frequency in Fourier space. This leads to the definition of the so-called Green operator $\widehat{\bm{\Gamma}}^0$ which has an explicit form in Fourier space:
\begin{equation}\label{eq:green_operator}
\widehat{\bm{\Gamma}}^0(\bm{\xi}) = \frac{\bm{\xi}\otimes\bm{\xi}}{c^0 |\bm{\xi}|^2}.
\end{equation}


\subsection{Discretization of the strong formulation by truncation}

In the previous section, the fields have been described initially using Fourier series, i.e.\/ with an infinite number of Fourier coefficients. In practice, this description is not compatible with a standard discretization of functions which requires a finite number of terms. The problem can be alternatively written by considering partial (truncated) Fourier series (of order $(M_x+1)\times\linebreak (M_y+1) \times (M_z+1)$) for the conductivity field and the fluctuation field:
\begin{align}
\varphi^{*}(x,y,z)&=\sum_{i=-M_x/2}^{M_x/2}\sum_{j=-M_y/2}^{M_y/2}\sum_{k=-M_z/2}^{M_z/2} \Phi^*_{ijk}\exp\left(\imath \xi_i^x x\right)\exp\left(\imath \xi_j^y y\right)\exp\left(\imath \xi_k^z z\right) \\
c(x,y,z)&=\sum_{i=-M_x/2}^{M_x/2}\sum_{j=-M_y/2}^{M_y/2}\sum_{k=-M_z/2}^{M_z/2} C_{ijk}\exp\left(\imath \xi_i^x x\right)\exp\left(\imath \xi_j^y y\right)\exp\left(\imath \xi_k^z z\right).
\end{align}
With this description, the local electric field $\mathbf{E}$ is also described by a partial (truncated) Fourier series (of order $(M_x+1)\times (M_y+1) \times (M_z+1)$).

Since the polarization tensor is given by equation~\eqref{eq:polarization}, as the product of the conductivity field and the local electric field, it is necessarily described by partial (truncated) Fourier series of order $(2M_x+1)\times (2M_y+1) \times (2M_z+1)$, i.e.,
\begin{align}
\tau_x(x,y,z)&=\sum_{i=-M_x}^{M_x}\sum_{j=-M_y}^{M_y}\sum_{k=-M_z}^{M_z} T^x_{ijk}\exp\left(\imath \xi_i^x x\right)\exp\left(\imath \xi_j^y y\right)\exp\left(\imath \xi_k^z z\right), \\
\tau_y(x,y,z)&=\sum_{i=-M_x}^{M_x}\sum_{j=-M_y}^{M_y}\sum_{k=-M_z}^{M_z} T^y_{ijk}\exp\left(\imath \xi_i^x x\right)\exp\left(\imath \xi_j^y y\right)\exp\left(\imath \xi_k^z z\right), \\
\tau_z(x,y,z)&=\sum_{i=-M_x}^{M_x}\sum_{j=-M_y}^{M_y}\sum_{k=-M_z}^{M_z} T^z_{ijk}\exp\left(\imath \xi_i^x x\right)\exp\left(\imath \xi_j^y y\right)\exp\left(\imath \xi_k^z z\right).
\end{align}

Consequently, the auxiliary problem reduces to the equation
\begin{multline}\label{eq:alternative_strong}
c_0 \sum_{i=-M_x/2}^{M_x/2}\sum_{j=-M_y/2}^{M_y/2}\sum_{k=-M_z/2}^{M_z/2} \left((\xi_i^x)^2 + (\xi_j^y)^2 + (\xi_k^z)^2 \right) \Phi^*_{ijk}\exp\left(\imath \xi_i^x x\right)\exp\left(\imath \xi_j^y y\right)\exp\left(\imath \xi_k^z z\right) \\
= - \sum_{i=-M_x}^{M_x}\sum_{j=-M_y}^{M_y}\sum_{k=-M_z}^{M_z} \imath \left(\xi_i^x T^x_{ijk} + \xi_j^y T^y_{ijk} + \xi_k^z T^z_{ijk} \right)\exp\left(\imath \xi_i^x x\right)\exp\left(\imath \xi_j^y y\right)\exp\left(\imath \xi_k^z z\right).
\end{multline}
In contrast with the previous case where functions where described by infinite Fourier series, equation~\eqref{eq:alternative_strong} does not possess a direct solution as it involves Fourier series of orders $(M_x+1)\times\linebreak (M_y+1) \times (M_z+1)$ and $(2M_x+1)\times (2M_y+1) \times (2M_z+1)$; the polarization tensor contains contributions at frequencies exceeding the maximum frequencies of the unknown field $\varphi^{*}$, which can be considered as \emph{aliasing}.

To overcome this issue, we will simply remove the high-frequencies of the polarization tensor, i.e.\/ we will assume that it is expressed as
\begin{equation}
\begin{aligned}\label{eq:pola_fourier_partial}
\tau_x(x,y,z)&=\sum_{i=-M_x/2}^{M_x/2}\sum_{j=-M_y/2}^{M_y/2}\sum_{k=-M_z/2}^{M_z/2} T^x_{ijk}\exp\left(\imath \xi_i^x x\right)\exp\left(\imath \xi_j^y y\right)\exp\left(\imath \xi_k^z z\right) \\
\tau_y(x,y,z)&=\sum_{i=-M_x/2}^{M_x/2}\sum_{j=-M_y/2}^{M_y/2}\sum_{k=-M_z/2}^{M_z/2} T^y_{ijk}\exp\left(\imath \xi_i^x x\right)\exp\left(\imath \xi_j^y y\right)\exp\left(\imath \xi_k^z z\right) \\
\tau_z(x,y,z)&=\sum_{i=-M_x/2}^{M_x/2}\sum_{j=-M_y/2}^{M_y/2}\sum_{k=-M_z/2}^{M_z/2} T^z_{ijk}\exp\left(\imath \xi_i^x x\right)\exp\left(\imath \xi_j^y y\right)\exp\left(\imath \xi_k^z z\right).
\end{aligned}
\end{equation}
With this assumption, equation~\eqref{eq:alternative_strong} now reads
\begin{multline}\label{eq:alternative_strong2}
c_0 \sum_{i=-M_x/2}^{M_x/2}\sum_{j=-M_y/2}^{M_y/2}\sum_{k=-M_z/2}^{M_z/2} \left((\xi_i^x)^2 + (\xi_j^y)^2 + (\xi_k^z)^2 \right) \Phi^*_{ijk}\exp\left(\imath \xi_i^x x\right)\exp\left(\imath \xi_j^y y\right)\exp\left(\imath \xi_k^z z\right) \\
= - \sum_{i=-M_x/2}^{M_x/2}\sum_{j=-M_y/2}^{M_y/2}\sum_{k=-M_z/2}^{M_z/2} \imath \left(\xi_i^x T^x_{ijk} + \xi_j^y T^y_{ijk} + \xi_k^z T^z_{ijk} \right)\exp\left(\imath \xi_i^x x\right)\exp\left(\imath \xi_j^y y\right)\exp\left(\imath \xi_k^z z\right),
\end{multline}
which leads to the relation between the Fourier series coefficients of the fluctuation and the components of the polarization
\begin{equation}\label{eq:alt_strong_fourier_series_coeff}
\Phi^*_{ijk} = - \imath \frac{\xi_i^x T^x_{ijk} + \xi_j^y T^y_{ijk} + \xi_k^z T^z_{ijk} }{c_0\left((\xi_i^x)^2 + (\xi_j^y)^2 + (\xi_k^z)^2\right)},
\end{equation}
for $i=-M_x/2,\dots,M_x/2$, $j=-M_y/2,\dots,M_y/2$ and $k=-M_z/2,\dots,M_z/2$ (except in the particular case $i=j=k=0$ where one has ${\Phi^*}_{000}=0$).

Then as previously, these coefficients are approximated by the coefficients of the discrete Fourier transform which need to be calculated using a discretization of the domain. The prismatic $\Omega$ is discretized uniformly with $(N_x+1)\times(N_y+1)\times(N_z+1)$ voxels so that the spatial scales are $\Delta x = L_x/(N_x+1)$, $\Delta y = L_y/(N_y+1)$ and $\Delta z = L_z/(N_z+1)$. The discrete Fourier transforms of the fields $\varphi^{*}$, $\tau_x$, $\tau_y$ and $\tau_z$ provide a numerical approximation of the Fourier series coefficients as shown previously. The only precaution that is required is to have a \emph{minimum number of ``integration'' points for the trapezoidal rule} (i.e.\/ $N_x$, $N_y$ and $N_z$) with respect to the \emph{number of modes} chosen in the description of the partial Fourier series (i.e.\/ $M_x$, $M_y$ and $M_z$), that is, again,
\begin{equation}
N_x \geq M_x, \quad \quad N_y \geq M_y, \quad \quad N_z \geq M_z.
\end{equation}
This finally leads to
\begin{equation}\label{eq:alt_strong_dft}
\widehat{\Phi^*}_{ijk} = - \imath \frac{\xi_i^x \widehat{T}^x_{ijk} + \xi_j^y \widehat{T}^y_{ijk} + \xi_k^z \widehat{T}^z_{ijk} }{c_0\left((\xi_i^x)^2 + (\xi_j^y)^2 + (\xi_k^z)^2\right)},
\end{equation}
for $i=-M_x/2,\dots,M_x/2$, $j=-M_y/2,\dots,M_y/2$ and $k=-M_z/2,\dots,M_z/2$, except in the particular case $i=j=k=0$ where one has $\widehat{\Phi^*}_{000}=0$.



It is worth noting that equation~\eqref{eq:alt_strong_dft} is equivalent to equation~\eqref{eq:strong_dft}, so the formulation in terms of the Green operator~\eqref{eq:green_operator} will naturally follow exactly the same way. The standard strong formulation and the alternative strong formulation with truncation will therefore lead to the same exact algorithm. The two approaches simply differ from their derivation: {
\begin{itemize}
\item In the classical approach, one starts with full Fourier series, the truncation of the Fourier series is a necessary consequence of the numerical estimation of the Fourier coefficients of $\bm{\uptau}$.
\item In the alternative approach, the truncation of the Fourier series is a starting discretization approximation, leading to a second necessary truncation of higher-order modes associated with higher-order frequencies (related to aliasing effects) that need to be disregarded to solve the auxiliary problem;
\end{itemize}
In both cases, the evaluation of the Fourier coefficients is done using discrete Fourier transforms with at least the same number of integration points than the number of modes considered in the description of the fields. It should be noted that one can consider more integration points than the number of modes (i.e.\/ in the case when $N_x > M_x$, $N_y > M_y$ and $N_z > M_z$); this is actually equivalent to disregard Fourier coefficients in the case of high frequencies (i.e.\/ for $|i|> M_x/2$, $|j|> M_y/2$ and $|k|> M_z/2$).} {It is important to note that the natural choice, i.e., $N_x = M_x$, $N_y = M_y$, $N_z = M_z$, is \emph{not necessarily} the optimal choice, in terms of quality of the solution. The influence of the ratio $M/N$ will be investigated in Section~\ref{sec:num} on a (specific) checkerboard microstructure.}

\subsection{Weak formulation with Galerkin approximation} We now consider the weak formulation of the auxiliary problem together with a Galerkin method. First, the unknown fluctuation field is described by a (partial) Fourier series,
\begin{equation}\label{eq:interp_func}
\varphi^{*}(x,y,z)=\sum_{i=-M_x/2}^{M_x/2}\sum_{j=-M_y/2}^{M_y/2}\sum_{k=-M_z/2}^{M_z/2} \Phi^*_{ijk}\exp\left(\imath \xi_i^x x\right)\exp\left(\imath \xi_j^y y\right)\exp\left(\imath \xi_k^z z\right).
\end{equation}
We are thus looking for the Fourier modes $\Phi^*_{ijk}$ for $i=-M_x/2,\dots,M_x/2$, $j=-M_y/2,\dots,M_y/2$ and $k=-M_z/2,\dots,M_z/2$, except in the particular case $i=j=k=0$ where one has ${\Phi^*}_{000}=0$ by construction.

The weak formulation of the local equilibrium~\eqref{eq:aux} reads
\begin{equation}\label{eq:weak_formulation}
\displaystyle \int_{\Omega} \Delta \widetilde{\varphi}^*(\mathbf{x}) v (\mathbf{x}) \mathrm{d}\Omega = \frac{1}{c_0}\int_{\Omega} \Div \bm{\tau}(\mathbf{x})v(\mathbf{x}) \mathrm{d}\Omega,
\end{equation}
where $\Delta$ is the Laplacian operator and $v(\mathbf{x})$ is a testing differentiable function. Following equation~\eqref{eq:interp_func}, these trial functions are the interpolation basis of $\varphi^{*}$, i.e.,
\begin{equation}
v_{lmn} (\mathbf{x}) = \exp\left(\imath \xi_l^x x\right)\exp\left(\imath \xi_m^y y\right)\exp\left(\imath \xi_n^z z\right)
\end{equation}
for $l=-M_x/2,\dots,M_x/2$, $m=-M_y/2,\dots,M_y/2$ and $n=-M_z/2,\dots,M_z/2$. The integrals involved in the weak formulation~\eqref{eq:weak_formulation}, that need to be calculated, are denoted by
\begin{equation}
\displaystyle I_{lmn}=\int_{\Omega} \Delta \widetilde{\varphi}^*(\mathbf{x}) v_{lmn} (\mathbf{x}) \mathrm{d}\Omega, \quad J_{lmn} = \int_{\Omega} \left(\Div \bm{\tau}(\mathbf{x}) \right) v_{lmn} (\mathbf{x}) \mathrm{d}\Omega.
\end{equation}
\paragraph*{Calculation of integral $I_{lmn}$.} The calculation of $I_{lmn}$ simply relies on the property of orthogonality of complex exponentials and is calculated analytically; one has
\begin{equation}
I_{lmn} = - L_x L_y L_z \left((\xi_{-l}^x)^2 + (\xi_{-m}^y)^2 + (\xi_{-n}^z)^2 \right) \Phi^*_{-l-m-n}
\end{equation}
since
\begin{equation}
\int_0^{L_x} \exp\left(\imath \xi_i^x x\right) \exp\left(\imath \xi_l^x x\right) \mathrm{d}x =
\left\{
\begin{array}{lll}
0 & \text{if} & i \neq l \\
L_x & \text{if} & i=-l.
\end{array}
\right.
\end{equation}

The calculation of $J_{lmn}$ can then be done using two different ways.
\paragraph*{Calculation of integral $J_{lmn}$.} Assuming that the polarization field is a periodic function described by Fourier series as in equation~\eqref{eq:pola_fourier_partial}, and using again the property of orthogonality of complex exponentials, one has
\begin{equation}
J_{lmn} = \imath L_x L_y L_z \left(\xi_{-l}^x T^x_{-l-m-n} + \xi_{-m}^y T^y_{-l-m-n} + \xi_{-n}^z T^z_{-l-m-n} \right).
\end{equation}

This leads to
\begin{equation}
\Phi^*_{-l-m-n} = \imath \frac{\xi_{-l}^x T^x_{-l-m-n} + \xi_{-m}^y T^y_{-l-m-n} + \xi_{-n}^z T^z_{-l-m-n}}{c_0 \left((\xi_{-l}^x)^2 + (\xi_{-m}^y)^2 + (\xi_{-n}^z)^2 \right)}
\end{equation}
or equivalently, using the relations $\xi_{-l}^x=-\xi_{l}^x$, $\xi_{-m}^y=-\xi_{m}^y$ and $\xi_{-n}^z=-\xi_{n}^z$,
\begin{equation}
\Phi^*_{lmn} = -\imath \frac{\xi_{l}^x T^x_{lmn} + \xi_{m}^y T^y_{lmn} + \xi_{n}^z T^z_{lmn}}{c_0 \left((\xi_{l}^x)^2 + (\xi_{m}^y)^2 + (\xi_{n}^z)^2 \right)}.
\end{equation}
This leads to the exact same formulation than the strong forms derived previously (equations~\eqref{eq:strong_fourier_series_coeff} and~\eqref{eq:alt_strong_fourier_series_coeff}). Therefore, as previously, one needs to evaluate these coefficients using discrete Fourier transforms (see equation~\eqref{eq:phi_T_DFT}).

\paragraph*{Alternative calculation of integral $J_{lmn}$.} Integral $J_{lmn}$ can be alternatively calculated without assuming that the polarization field is described by Fourier series. In that case, the integral is directly approximated using a trapezoidal rule and will thus lead to the same discrete Fourier transform and thus the exact same formulation.

\subsection{Iterative scheme} The basic scheme of~\cite{moulinec_numerical_1998} then consists in a fixed point in order to determine $\bm{\uptau}(\bm{x})$ solution of the problem. Irrespective of the method considered (strong form, alternative strong form or weak form), the resolution of the auxiliary problem leads to the same relation between the Fourier coefficients of the unknown fluctuation and the Fourier coefficients of the polarization tensor. The following fixed-point algorithm then follows:
\begin{align}
\left\{
\begin{array}{rl}
{\rm Initialization} & \quad \mathbf{E}^0(\mathbf{x}) = \bar{\mathbf{E}} \\
& \quad \mathbf{J}^0(\mathbf{x}) = c(\mathbf{x})\mathbf{E}^0(\mathbf{x}) \\[0.25cm] {\rm Iterate~i+1} & \quad \mathbf{E}^i~{\rm and}~\mathbf{J}^i~{\rm being~known} \\
{\rm (a)} & \quad \bm{\uptau}^i (\mathbf{x}) = \mathbf{J}^i(\mathbf{x})-c^0\mathbf{E}^i(\mathbf{x}) \\
{\rm (b)} & \quad \widehat{\bm{\tau}}^i = \text{DFT3}\left(\bm{\tau}^i \right) \\
{\rm (c)} & \quad {\rm Convergence~test} \\
{\rm (d)} & \quad \widehat{\bm{E}}^{i+1}(\bm{\xi}) = - \widehat{\mathbf{\Gamma}}^0(\bm{\xi}):\widehat{\bm{\tau}}^i(\bm{\xi})~\forall \bm{\xi} \neq \mathbf{0}\quad {\rm and}~\widehat{\bm{E}}^{i+1}(\mathbf{0}) = \bar{\mathbf{E}} \\
{\rm (e)} & \quad{\bm{E}}^{i+1} = \text{IDFT3} \left(\widehat{\bm{E}}^{i+1} \right) \\
& \quad {\mathbf{E}}^{i+1}(\mathbf{x}) \quad \text{expressed using Fourier series} \\
{\rm (f)} & \quad \mathbf{J}^{i+1}(\mathbf{x}) = c(\mathbf{x}){\mathbf{E}}^{i+1}(\mathbf{x})
\end{array}\right.
\end{align}
Convergence is reached when $\mathbf{J}^{i+1}$ is in equilibrium, i.e.\/ its divergence is null. In practice, the iterative procedure is stopped when the error serving to check convergence~\cite{moulinec_numerical_1998},
\begin{equation}\label{eq:conv}
e^i = \frac{ \left(\moy{ \left\Vert \Div(\mathbf{J}^i) \right\Vert^2 } \right)^{1/2} }{ \left\Vert \moy{\mathbf{J}^i} \right\Vert},
\end{equation}
is smaller than a prescribed (small) value.


\section{Representation of solution fields and non-uniform discretization}\label{sec:representation}

\subsection{Representation of solution fields}

The FFT-based method is based on the approximate calculation of Fourier series coefficients (equation~\eqref{eq:coeff_fourier_series_3d}) using the discrete Fourier transform on a \emph{regular} grid. The role of the (uniform) spatial discretization of the domain is therefore to calculate numerically the Fourier coefficients used in the description of the fields (equation~\eqref{eq:dft_3d}), based on a trapezoidal (and thus uniform) quadrature rule (with $N+1$ equi-spaced points) \cite{trefethen2014exponentially}, i.e.,
\begin{equation}\label{eq:dft_unif}
\widehat{\bm{F}}= \DFT(\bm{f}),
\end{equation}
where $\bm{f}$ contains the nodal values of function $f$ on a uniform grid. By taking advantage of the relation~\eqref{eq:dft_fourier_series} between the discrete Fourier coefficients (i.e.\/ $\widehat{\bm{F}}$) and the Fourier series coefficients (i.e.\/ ${\bm{F}}$), the reconstruction of function $f$ then follows \emph{in the whole domain}
\begin{equation}\label{eq:series_nu}
f(x) = \sum_{i=-M/2}^{M/2} F_i \exp\left(\imath \xi_i x\right),\quad \quad \quad \forall x\in[0,L].
\end{equation}

Usually, the number of Fourier modes ($M+1$) is taken equal to the number of spatial nodes ($N+1$), i.e.\/ $N=M$, and then the inverse Fourier transform can then be used to evaluate, \emph{on the same spatial uniform grid used to approximate the integrals}, function $f$ as,
\begin{equation}
\bm{f}= \IDFT(\widehat{\bm{F}}).
\end{equation}
In the literature, the solution fields are then generally represented from these nodal values, either using uniform values in a voxel or a linear interpolation between nodal values (see e.g.~\cite{moulinec_numerical_1998}). However, as the solution fields are described, by construction, by Fourier series, they are defined in the whole domain as shown by equation~\eqref{eq:series_nu}. This implies that representing the field only at the nodes, or uniformly in the voxels, cannot capture variations between spatial nodes, and may therefore hide some spurious numerical effects such as Gibbs effects. {This can be overcome by representing the solution fields on a finer grid, using directly equation~\eqref{eq:series_nu}. In the context of homogenization, such a refined representation does not modify macroscopic quantities, which depend only on spatial averages. Representing the solution field on a finer grid is therefore not motivated by the computation of effective properties, but rather by situations where an accurate evaluation of the solution fields is required (e.g., evaluation of local criteria, stress-assisted diffusion, or phase-field models).}

Let us denote by $x^{\fine}_i$ the grid points of a uniform finer grid (with $N^{\fine}+1$ points) and by $f^{\fine}_i=f(x^{\fine}_i)$ the value of function $f$ on this refined grid. We denote by $\bm{f}^{\fine}$ the array of size $N^{\fine}+1$ containing the values $f^{\fine}_i$. From equation~\eqref{eq:series_nu} one simply has
\begin{equation}
f^{\fine}_i = \sum_{j=-M/2}^{M/2} F_j \exp\left(\imath \xi_j x^{\fine}_i \right)
\end{equation}
which can be written as
\begin{equation}
\bm{f}^{\fine} = \bm{A}\cdot \bm{F},
\end{equation}
where $\bm{A}$ is the $(N^{\fine}+1)\times (M+1)$ array defined by
\begin{equation}
A_{ij} = \exp\left(\imath \xi_j x^{\fine}_i \right).
\end{equation}
The determination of $\bm{f}^{\fine}$ then requires a matrix-vector operation involving matrix $\bm{A}$ for which the numerical complexity scales in $\mathcal{O}(N^{\fine}\times M)$. This evaluation can be considerably reduced by using the inverse discrete Fourier transform with a modified vector for the Fourier coefficients, denoted by $\widehat{\bm{F}}^{\fine}$, taken null for the extra high frequency modes:
\begin{equation}
\widehat{F}^{\fine}_i =
\begin{cases}
\widehat{F}_i & \text{if } i \leq M/2 \\
0 & \text{if } i> M/2.
\end{cases}
\end{equation}
In that case the numerical complexity would scale in $\mathcal{O}(N^{\fine}\log (N^{\fine}))$.

{It is important to note that representing the solution on a finer grid than the resolution grid is \emph{not} a reinterpolation of the solution; it is simply the exact evaluation of the same spectral approximation at new points. No additional interpolation procedure is performed; the underlying trigonometric polynomial remains unchanged.}



\subsection{Non-uniform discretization} {The representation of the local fields on a finer grid entails a non-negligible computational cost, scaling as $\mathcal{O}(N^{\fine}\log N^{\fine})$. In situations where only localized regions of the RVE are of interest, the computational cost can be reduced by resorting to non-uniform FFT evaluations on \emph{selected points}.}

{Indeed, as} function $f$ is described by the Fourier series~\eqref{eq:series_nu}, it can be calculated in any grid point, even in an \emph{non-uniform grid}. Let us denote by $x^{\rnu}_i$ the grid points of a non-uniform grid (with $N^{\rnu}+1$ points $x^{\rnu}_i \in [0,L[$) and by $f^{\rnu}_i=f(x^{\rnu}_i)$ the value of function $f$ on this non-uniform grid. We denote by $\bm{f}^{\rnu}$ the array of size $N^{\rnu}+1$ containing the values $f^{\rnu}_i$. From equation~\eqref{eq:series_nu} one simply has
\begin{equation}
f^{\rnu}_i = \sum_{j=-M/2}^{M/2} F_j \exp\left(\imath \xi_j x^{\rnu}_i \right)
\end{equation}
which can be written as
\begin{equation}\label{eq:A_nufft}
\bm{f}^{\rnu} = \bm{A}\cdot \bm{F},
\end{equation}
where $\bm{A}$ is the $(N^{\rnu}+1)\times (M+1)$ array defined by
\begin{equation}
A_{ij} = \exp\left(\imath \xi_j x^{\rnu}_i \right).
\end{equation}
As previously the numerical complexity of this operation is in $\mathcal{O}(N^{\rnu}\times M)$. However, in the case where $N^{\rnu}=M$, but still with a non-uniform grid, the evaluation of function $f$ can be done with the same numerical complexity of fast Fourier transforms at nonequispaced points, which is reffered to as NFFT, for nonequispaced fast Fourier transform, or nonuniform fast Fourier transform (NUFFT) \cite{duijndam1999nonuniform,greengard2004accelerating}.

Knowing the Fourier series coefficients of function $f$, i.e.\/ $\bm{F}$ (or equivalently the discrete Fourier coefficients $\widehat{\bm{F}}$), the NUFFT defined by equation~\eqref{eq:A_nufft} permits to evaluate \emph{fastly} function $f$ at any points of an arbitrary grid. {It thus can be used simply as a post-treatment to represent the solution field with a fine resolution in areas where strong gradients are expected instead of using a very fine mesh everywhere in the domain (which would entail a non-negligible computational cost).}

In addition, as equation~\eqref{eq:A_nufft} provides the matrix expression of the NFFT, it could be inverted so as to provide the value of the Fourier coefficients as a function of the grid point values of function $f$ and can be calculated efficiently as standard inverse FFT~\cite{kircheis2019direct}. It thus defines an alternative way of calculating the Fourier coefficients on a non-uniform grid which can be interpreted as a modified quadrature rule. However it is expected not to improve the results obtained using the trapezoidal rule which is exponentially convergent for a periodic function~\cite{trefethen2014exponentially}.


\section{Numerical experiments}\label{sec:num}
{The aim of this section is to study, for a specific problem admitting an analytical solution, the influence of the number of Fourier modes and grid discretization on the quality of the solution fields. We investigate (i) the influence of the grid discretization ($N$) for a given number of modes ($M$) and (ii) the influence of the number of modes ($M$) for a given grid discretization ($N$). Indeed, as explained in Section~\ref{sec:formulations}, the number of grid points $N$ can be chosen independently of the number of mode $M$ with the sole condition $N \geq M$; the relevance of the standard choice $M=N$ will be notably studied.} A particular attention will be put on the representation of the fields as explained in Section~\ref{sec:representation}.

\subsection{Description of the problem} We consider the classical problem of conductivity for checkerboard-type microstructures which is of interest for the method since this problem has an analytical exact solution~\cite{dykhne1971conductivity,obnosov_periodic_1999,craster_four-phase_2001} (see also~\cite{bellis_eigendecomposition-based_2020}). This problem was notably investigated in previous studies since it is highly challenging for the accurate description, without spurious oscillations of the local fields~\cite{willot_fourier-based_2014,dorn2019lippmann,morin_periodic_2021}. {This is an extreme case due to the presence of sharp interfaces and corners; it is thus a good benchmark to assess the influence of the ratio $M/N$.}

We consider a 2D periodic square domain $[0,1]\times[0,1]$ discretized with $(N+1) \times (N+1)$ grid points, using the notations of Section~\ref{sec:formulations}. We consider that the isotropic conductivity field $c(x_1,x_2)$ is heterogeneous and forms a checkerboard-type microstructure with values $c_1=1$ and $c_2=100$ as shown in Figure~\ref{fig:th}a and defined as
\begin{equation}
c(x_1,x_2)=
\begin{cases}
c_2 & \text{if }x_1 < 0.5~\text{and}~x_2 < 0.5 \\
c_1 & \text{elsewhere}.
\end{cases}
\end{equation}
{Thus, the grid point at $(0,0)$ is at the left bottom corner of the included square. One might remark that other choice of alignment between the grid and the microstructure might lead to slightly different results. This influence, sometimes referred as the ``egg-box'' effect (see, for instance, \cite{briggs1996real}), is not investigated here.}

The values $\bar{E}_1=1$ and $\bar{E}_2=0$ are considered for the prescribed macroscopic electric field $\bar{E}$. Convergence of the iterative procedure is reached when the error defined by equation~\eqref{eq:conv} is smaller than the prescribed value $10^{-8}$. The normalized electric field $E_1/\VbE$ obtained from the theoretical solution of~\cite{craster_four-phase_2001} is represented in Figure~\ref{fig:th}b.

\begin{figure}[!ht]%fig1
\centering
\subfloat[]{\includegraphics[width=6.5cm]{figures/microstructure.png}}\hspace{0.5cm}
\subfloat[]{\includegraphics[width=6.5cm]{figures/Th_E1.png}}
\caption{Checkerboard problem considered. (a) Distribution of the conductivity field $c(x_1,x_2)$ field and (b) Distribution of the theoretical normalized electric field $E_1/\VbE$.}\label{fig:th}
\end{figure}

For the numerical calculations, several values for the number of grid points will be considered, $N=[1024,~512,~256,~128]$, and several values for the number of Fourier modes will also be considered, $M=[1024,~512,~256,~128]$. The field will be represented on both their ``original'' grid (of size $(N+1) \times (N+1)$) as well as on a refined grid as explained in Section~\ref{sec:representation}, of size $(N^{\fine}+1) \times (N^{\fine}+1)$ with $N^{\fine}=8192$.

Qualitative results are first presented (Section~\ref{sec:gridN} and~\ref{sec:modeM}) to illustrate separately the effect of mesh refinement and Fourier modes refinement on the local fields. Quantitative comparisons are then provided in Section~\ref{sec:quantitative}.


\subsection{Influence of the grid discretization (\texorpdfstring{$N$}{N}) for a given number of modes (\texorpdfstring{$M$}{M})}\label{sec:gridN}

We consider the case $M=128$ Fourier modes and we vary the number of grid nodes ($N=[1024,$ 512, 256, $128]$), which therefore consists in modifying the number of integration points for calculating numerically the Fourier coefficients. First, the distribution of the component $E_1$ of the electric field is represented in Figure~\ref{fig:conductivity_E1_MvsN} for several values of the grid nodes. In each case, the distribution is represented on both the calculation grid (i.e.\/ the coarse grid of size $(N+1)\times (N+1)$) and the fine grid (of size $(N^{\fine}+1)\times (N^{\fine}+1)$).

Overall, the results in all cases are close to the analytical solution, but spurious oscillations are observed. Interestingly, these oscillations are not visible in the case $(M=128,~N=128)$ when the solution is represented in the coarse grid, i.e.\/ on the grid used for the calculation of the FFT (Figure~\ref{fig:conductivity_E1_MvsN}a), but they are visible when the solution is represented in the fine grid (Figure~\ref{fig:conductivity_E1_MvsN}b); this means that, when the solution is only plotted on the nodes of calculation, oscillations between these nodes may not be seen, but between grid points these oscillations do exist. When increasing the number of integration points, i.e.\/ $N$, the solution seems to converge, and the solutions represented on the discretization grid and the fine grid coincide.

\begin{figure}[!ht]%fig2
\centering
\subfloat[]{\includegraphics[width=4.25cm]{figures/Coarse_M_128_N_128.png}}\hspace{0.35cm}
\subfloat[]{\includegraphics[width=4.25cm]{figures/Coarse_M_128_N_256.png}} \hspace{0.35cm}
\subfloat[]{\includegraphics[width=4.25cm]{figures/Coarse_M_128_N_1024.png}} \\
\subfloat[]{\includegraphics[width=4.25cm]{figures/Fine_M_128_N_128.png}} \hspace{0.35cm}
\subfloat[]{\includegraphics[width=4.25cm]{figures/Fine_M_128_N_256.png}}\hspace{0.35cm}
\subfloat[]{\includegraphics[width=4.25cm]{figures/Fine_M_128_N_1024.png}}
\caption{Influence of the grid discretization ($N$) for a given number of modes ($M$) on the distribution of the normalized electric field $E_1/\VbE$. The solution is represented on both the discretization grid (coarse mesh) and the fine grid (fine mesh).}\label{fig:conductivity_E1_MvsN}
\end{figure}


\begin{figure}[!h]%fig3
\centering
\includegraphics[width=8cm]{figures/line_M128_N_varies}
\caption{Line plot of the normalized electric field $E_1/\VbE$ interface between the two phases ($x_2=0$) for a fixed number of modes ($M=128$). The continuous lines corresponds to the solution on the fine grid and the triangles corresponds to the solution on the grid used for the FFT calculations.}\label{fig:line_M128_N_varies}
\end{figure}


Additionally, an enlarged view, in the vicinity of the interface between the two phases, is plotted in Figure~\ref{fig:line_M128_N_varies}, on the line $x_2=0$. In the case $(M=128,~N=128)$, the oscillations are not visible if one only represents the gridpoint values of the field (represented by blue triangles) but they become visible when the solution is represented in the fine grid. It is interesting to note that the solution field stabilizes when $N$ increases. This is consistent with the role played by $N$, which is to calculate approximately integrals: for a fixed number of modes $M$, the increase of $N$ allows a better approximation of the integrals, but this evaluation is expected to reach a stationary value when increasing the number of integration points.


\subsection{Influence of the number of modes (\texorpdfstring{$M$}{M}) for a given grid discretization (\texorpdfstring{$N$}{N})}\label{sec:modeM}

We now consider the case $N=1024$ and we vary the number of Fourier modes ($M=[1024,~256,~128]$), which therefore consists in removing some high frequencies in the solution. First, the distribution of the component $E_1$ of the electric field is represented in Figure~\ref{fig:conductivity_E1_NvsM} for several values of the grid nodes. In each case, the distribution is represented on both the calculation grid (i.e.\/ the coarse grid of size $(N+1)\times (N+1)$) and the fine grid (of size $(N^{\fine}+1)\times (N^{\fine}+1)$).

First, the distribution seems quite similar irrespective of the resolution considered for the representation. As previously, spurious oscillations are observed on the solution fields, which tend to decrease when the number of Fourier modes increases.

\begin{figure}[!ht]
\centering
\subfloat[]{\includegraphics[width=4.25cm]{figures/Coarse_M_128_N_1024.png}}\hspace{0.35cm}
\subfloat[]{\includegraphics[width=4.25cm]{figures/Coarse_M_256_N_1024.png}} \hspace{0.35cm}
\subfloat[]{\includegraphics[width=4.25cm]{figures/Coarse_M_1024_N_1024.png}} \\
\subfloat[]{\includegraphics[width=4.25cm]{figures/Fine_M_128_N_1024.png}} \hspace{0.35cm}
\subfloat[]{\includegraphics[width=4.25cm]{figures/Fine_M_256_N_1024.png}}\hspace{0.35cm}
\subfloat[]{\includegraphics[width=4.25cm]{figures/Fine_M_1024_N_1024.png}}
\caption{Influence of number of modes ($M$) for a given grid discretization ($N$) on the distribution of the normalized electric field $E_1/\VbE$. The solution is represented on both the discretization grid (coarse mesh) and the fine grid (fine mesh).}\label{fig:conductivity_E1_NvsM}
\end{figure}

An enlarged view of the solution field, at the vicinity of the interface between the two phases, is plotted in Figure~\ref{fig:line_N1024_M_varies}, on the line $x_2=0$. This representation allows to see the oscillatory nature of all solutions, with a period that depends on the number of Fourier modes ($M$). This is consistent with the interpretation of the Fourier modes which are associated with increasing frequencies. Increasing the number of modes allows to restrict the Gibbs phenomenon close to the interface but it may increase the amplitude of the oscillations locally. Therefore, even though the numerical solution seems to get closer to the analytical one when increasing $M$, oscillations with higher frequencies appear, so one cannot conclude on the improvement of quality of the solution field with the increase of the modes.

\begin{figure}[!h]
\centering
\includegraphics[width=8cm]{figures/line_N1024_M_varies_zoom}
\caption{Line plot of the normalized electric field $E_1/\VbE$ interface between the two phases ($x_2=0$) for a fixed grid discretization ($N=1024$). The continuous lines corresponds to the solution on the fine grid and the triangles corresponds to the solution on the grid used for the FFT calculations.}\label{fig:line_N1024_M_varies}
\end{figure}

\subsection{Quantitative comparison with the analytical solution}\label{sec:quantitative}

A quantitative comparison between numerical solutions and the analytical one is made by studying the relative error in energy norm defined as~\cite{kelly1983posteriori}
\begin{equation}\label{eq:err}
\err = \frac{\displaystyle \int_{\Omega} (\mathbf{E}^{\rm th}-\mathbf{E}^{\rnum}):\mathbf{c}:(\mathbf{E}^{\rm th}-\mathbf{E}^{\rnum}) {\mathrm{d} \Omega} }{\displaystyle \int_{\Omega} \mathbf{E}^{\rm th}:\mathbf{c}:\mathbf{E}^{\rm th} {\mathrm{d}\Omega }},
\end{equation}
where $\mathbf{E}^{\rm th}$ is the theoretical electric field and $\mathbf{E}^{\rnum}$ is the numerical electric field obtained by FFT. Following Galerkin's framework, this error should be minimal in the chosen function subspace, i.e.\/ for the number of modes $M$.

The error is represented in Figure~\ref{fig:err}a as a function of the number of grid nodes $N$ (for a given number of modes $M$). For a given number of modes $M$, the error decreases with the number of nodes (which plays here the role of integration points) and tends to an asymptotic value, in agreement with the results of Figure~\ref{fig:line_M128_N_varies}. This confirms that the use of FFT (approximate calculation of integrals using a trapezoidal rule, which can be seen as a reduced integration) induces an error on the solution fields. The stabilization of the error indicates that the reduced-integration-induced error has been suppressed. Residual error comes from the truncated Fourier series approximation, and decreases with the number of modes $M$. This confirms that the approximate calculation of integrals using a trapezoidal rule, for a given function (imposed by the number of modes), reaches a stationary value when increasing the number of integration points $N$.

Then, the error is represented in Figure~\ref{fig:err}b as a function of the number of modes $M$ (for a given number of nodes $N$). Overall, the error tends to decrease when the number of modes $M$ increases (for a fixed $N$) until a limit value $M_{\rm lim}$ above which the error increases. This value $M_{\rm lim}$ depends on the number of nodes $N$ and its order of magnitude is of about $M_{\rm lim} \sim 3N/4$. This peculiar result indicates that, for this problem, the optimal number of modes $M$ is not necessarily equal to the number of integration points $N$: a truncation of the Fourier series (i.e.\/ $M < N$), which corresponds to disregard Fourier coefficients for high frequencies, can lead to a better solution field. This is equivalent to filtering high frequency components known to degrade solution fields near regions of discontinuity (see e.g.~\cite{djaka2015numerical}).

\begin{figure}[!ht]
\centering
\subfloat[]{\includegraphics[width=6.5cm]{figures/error_M}}\hspace{0.5cm}
\subfloat[]{\includegraphics[width=6.5cm]{figures/error_N}}
\caption{Error analysis on the local fields. (a) Influence of the number of integration points $N$ and (b) Influence of the number of modes $M$. }\label{fig:err}
\end{figure}

In addition, since FFT-based methods are generally used to calculate the overall properties of composites, it is also interesting to study the influence of the discretization on the macroscopic conductivity $\bar{c}$ which arises in the relation
\begin{equation}
\left<\mathbf{J}\right>=\bar{c}\left<\mathbf{E}\right>,
\end{equation}
where $\left<\mathbf{J}\right>$ and $\left<\mathbf{E}\right>$ are the volume averages of the electric current and electric field, respectively. The macroscopic conductivity $\bar{c}$, obtained from the numerical calculation, denoted by $\bar{c}^{\rnum}$, is simply obtained using the Fourier coefficients of the electric current and electric field associated to the first Fourier frequencies, since the first Fourier frequencies are associated to the average of the considered field. For the loading considered, one thus has (using the notation of equation~\eqref{eq:0-freq})
\begin{equation}\label{eq:num_conduct}
\bar{c}^{\rnum} = \frac{\widehat{J}^x_{00}}{\widehat{E}^x_{00}}.
\end{equation}
The theoretical macroscopic conductivity, denoted by $\bar{c}^{\rm th}$, is given by~\cite{craster_four-phase_2001}:
\begin{equation}
\bar{c}^{\rm th} = \sqrt{\frac{c_1+3c_2}{3c_1+c_2}}.
\end{equation}
Then we define the relative error for the macroscopic conductivity as
\begin{equation}
\err_{\rm macro} = \frac{\left| \bar{c}^{\rm th} - \bar{c}^{\rnum} \right|}{\bar{c}^{\rm th}}.
\end{equation}

First, the error for the macroscopic conductivity is represented in Figure~\ref{fig:err_macro}a as a function of the number of grid nodes $N$ (for a fixed number of modes $M$). A stationary value of the error is observed when increasing the number of nodes, which act here as integration points. This behaviour is consistent with the convergence of the numerical integration as the number of grid points increases. For the considered test case, the macroscopic error is minimal when the integration grid matches the spectral truncation, i.e.\/ $(N = M)$. This behaviour is attributed to a partial compensation between truncation and numerical integration errors. It is therefore problem-dependent and not expected to hold in general.

\begin{figure}[!ht]
\centering
\subfloat[]{\includegraphics[width=6.5cm]{figures/error_macro_M}}\hspace{0.5cm}
\subfloat[]{\includegraphics[width=6.5cm]{figures/error_macro_N}}
\caption{Error analysis on the macroscopic conductivity. (a) Influence of the number of integration points $N$ and (b) Influence of the number of modes $M$. }\label{fig:err_macro}
\end{figure}

Then, the error for the macroscopic conductivity is represented in Figure~\ref{fig:err_macro}b as a function of the number of modes $M$ (for a given number of nodes $N$). Contrary to the error related to the local fields, the error for the macroscopic conductivity keeps decreasing while increasing the number of modes. Therefore, the estimation of macroscopic properties is not affected by the high frequencies, in contrast to the local fields, which may be affected by Gibbs phenomenon, while macroscopic properties correspond to the zero Fourier mode of the solution, as shown by equation~\eqref{eq:num_conduct}. Nevertheless, the zero Fourier mode results from the solution of the Lippmann--Schwinger equation, which involves a convolution in Fourier space and therefore depends on the full spectral content of the fields. Thus, this result is a numerical observation, which might be case-dependent.


\subsection{{Influence of the contrast on the optimal ratio \texorpdfstring{$M/N$}{M/N}}}

The above results show that the choice $M=N$ is not necessarily the optimal choice in terms of quality of the local fields, which is a peculiar and unexpected emerging result. For this specific case, i.e., checkerboard microstructure with a contrast $c_2/c_1=100$, the optimal value $M/N$ associated to the relative error in energy norm is almost constant and lies between 0.67 and 0.72. This emerging result is not universal as it is expected to depend on both the microstructure and the contrast. Indeed, in the limit case of a homogeneous material, i.e., $c_2/c_1=1$, no oscillations are expected and thus the improvement of the solution field obtained by taking $M$ lower than $N$ is not expected to hold in that case. To assess quantitatively the effect of the contrast on the optimal value $M/N$, we consider several values of the contrast $c_2/c_1=[2,~10,~100,~1000]$ and only one value of the number of grid points $N=512$ is considered.

The error is represented in Figure~\ref{fig:contrast} as a function of the number of modes $M$ (for a given number of nodes $N=512$). The error decreases when the number of modes $M$ increases and, depending on the value of the contrast, it can reach a limit value $M_{\rm lim}$ above which the error increases. As expected the optimal value of the ratio $M/N$ is close to 1 for low-contrasted materials, and then it decreases when the contrast increases. For this specific problem this is attributed to the presence of oscillations, whose occurrence depends on the contrast.

\begin{figure}[!h]
\centering
\includegraphics[width=8cm]{figures/error_N_contrast}
\caption{{Influence of the contrast upon the error analysis on the local fields.}}\label{fig:contrast}
\end{figure}

The emergence of this optimal value $M/N$ cannot be seen as a universal result and is expected to be problem-dependent: the only conclusion that can be drawn with this analysis is that the choice $M=N$ is not necessarily the optimal choice for local fields.


\section{Discussion on spurious oscillations}\label{sec:disc}

In this section, we discuss the occurrence of spurious oscillations, which are the main source of discrepancies between the numerical and analytical solutions, and present remedies based either on discrete Green operators or on smoothing techniques.

\subsection{On generalized Green operators}\label{sec:green}

The three approaches considered in Section~\ref{sec:formulations} all rely on Fourier series together with standard derivative rules which leads to the so-called ``continuous'' Green operator defined by equation~\eqref{eq:green_operator}. It has been shown in the literature that modified ``discrete'' Green operators, obtained using finite differences or finite elements approximations, generally improve the convergence rate of the iterative procedure as well as the quality of the local fields, by suppressing most of the spurious oscillations~\cite{willot_fourier-based_2014,zeman2017finite,schneider_fft-based_2017}. \emph{It should be noted that the developments of Green operators are thus not compatible as is with the description of fields using Fourier series}; the use of discrete Fourier transforms simply arise in the resolution of the linear system obtained using either finite differences or finite elements (which defines a circulant matrix in the case of periodic boundary conditions), leading to the so-called preconditioned FFT-accelerated finite element solvers.

In this section we revisit the construction of the discrete Green operator by introducing \emph{modified derivative rules} into the framework of Section~\ref{sec:formulations}, still in the framework of Fourier series. Indeed, instead of using the standard spectral derivative rules, we will consider a difference (continuous) operator with a given spacing (not necessarily equal to the grid size), based on difference quotient. In order to retrieve~\cite{willot_fourier-based_2014}'s scheme and generalize it, we will consider a forward-type derivative scheme for the electric potential and a backward-type derivative scheme for the current. Other choice for these derivative rules would lead to another modified Green operators. This enables the use of discrete Green operators within the Fourier series framework (i.e.\/ using spectral interpolation) and allows to determine whether oscillations originate from the derivative rule or the interpolation (spectral) functions.

\paragraph*{Modified derivative rules for the electric potential.} We consider the following \emph{continuous} forward derivative scheme for the electric potential
\begin{equation}
\begin{aligned}\label{eq:forward_scheme}
\delta^{x+}_{h_x}\left[\varphi^{*} \right](x,y,z) &= \frac{\varphi^{*}(x+h_x,y,z)-\varphi^{*}(x,y,z)}{h_x}\\
\delta^{y+}_{h_y}\left[\varphi^{*} \right](x,y,z) &= \frac{\varphi^{*}(x,y+h_y,z)-\varphi^{*}(x,y,z)}{h_y}\\
\delta^{z+}_{h_z}\left[\varphi^{*} \right](x,y,z) &= \frac{\varphi^{*}(x,y,z+h_z)-\varphi^{*}(x,y,z)}{h_z}
\end{aligned}
\end{equation}
where $h_x$, $h_y$ and $h_z$ are the spatial steps considered in each spatial directions.

\paragraph*{Modified derivative rules for the electric current.} We then introduce a \emph{continuous} backward derivative scheme
\begin{equation}
\begin{aligned}\label{eq:backward_scheme}
\delta^{x-}_{h_x}\left[J_x \right](x,y,z) &= \frac{J_x(x,y,z)-J_x(x-h_x,y,z)}{h_x}\\
\delta^{y-}_{h_y}\left[J_y \right](x,y,z) &= \frac{J_y(x,y,z)-J_y(x,y-h_y,z)}{h_y}\\
\delta^{z-}_{h_z}\left[J_z \right](x,y,z) &= \frac{J_z(x,y,z)-J_z(x,y,z-h_z)}{h_z}.
\end{aligned}
\end{equation}

Some comments are in order:
\begin{itemize}
\item In the case $h_x=\Delta x$, $h_y=\Delta y$ and $h_z=\Delta z$, the derivative rules are very similar to the finite difference scheme of~\cite{willot_fourier-based_2014}, but the difference is that \emph{it is defined even between nodal values} in the present formulation.
\item In the case $h_x\to 0$, $h_y\to 0$, $h_z\to 0$ the derivative rules considered correspond to the standard derivative rules considered previously.
\end{itemize}

Assuming that $\varphi^{*}$ is described by a Full Fourier series, the forward derivative scheme~\eqref{eq:forward_scheme} in the $x-$direction reads
\begin{equation}
\begin{aligned}
\delta^{x+}_{h_x}\left[\varphi^{*} \right](x,y,z)& = \sum_{i=-\infty}^{+\infty}\sum_{j=-\infty}^{+\infty}\sum_{k=-\infty}^{+\infty} \Phi^*_{ijk} \frac{\exp\left(\imath \xi_i^x (x+h_x)\right) - \exp\left(\imath \xi_i^x x\right) }{h_x} \exp\left(\imath \xi_j^y y\right)\exp\left(\imath \xi_k^z z\right) \\
& = \sum_{i=-\infty}^{+\infty}\sum_{j=-\infty}^{+\infty}\sum_{k=-\infty}^{+\infty} \Phi^*_{ijk} \frac{\exp\left(\imath \xi_i^x h_x\right) - 1}{h_x} \exp\left(\imath \xi_i^x x\right) \exp\left(\imath \xi_j^y y\right)\exp\left(\imath \xi_k^z z\right).
\end{aligned}
\end{equation}

Assuming that $\tau_x$ is described by a full Fourier series, the backward derivative scheme~\eqref{eq:backward_scheme} in the $x-$direction reads
\begin{equation}
\begin{aligned}
\delta^{x-}_{h_x}\left[\tau_x \right](x,y,z)& = \sum_{i=-\infty}^{+\infty}\sum_{j=-\infty}^{+\infty}\sum_{k=-\infty}^{+\infty} T^x_{ijk} \frac{\exp\left(\imath \xi_i^x x \right) - \exp\left(\imath \xi_i^x (x-h_x)\right) }{h_x} \exp\left(\imath \xi_j^y y\right)\exp\left(\imath \xi_k^z z\right) \\
&= \sum_{i=-\infty}^{+\infty}\sum_{j=-\infty}^{+\infty}\sum_{k=-\infty}^{+\infty} T^x_{ijk} \frac{ 1- \exp\left(-\imath \xi_i^x h_x\right) }{h_x} \exp\left(\imath \xi_i^x x\right) \exp\left(\imath \xi_j^y y\right)\exp\left(\imath \xi_k^z z\right).
\end{aligned}
\end{equation}

Then, considering these derivative rules and the description of the fields using full Fourier series as in Section~\ref{sec:formulations}, the auxiliary problem reduces, after a few straightforward calculations involving trigonometric relations, to
\begin{multline}
c_0 \sum_{i=-\infty}^{+\infty}\sum_{j=-\infty}^{+\infty}\sum_{k=-\infty}^{+\infty} \left((\alpha_i^x)^2 + (\alpha_j^y)^2 + (\alpha_k^z)^2 \right) \Phi^*_{ijk}\exp\left(\imath \xi_i^x x\right)\exp\left(\imath \xi_j^y y\right)\exp\left(\imath \xi_k^z z\right) \\
= - \sum_{i=-\infty}^{+\infty}\sum_{j=-\infty}^{+\infty}\sum_{k=-\infty}^{+\infty} \imath \left(\beta_i^x T^x_{ijk} + \beta_j^y T^y_{ijk} + \beta_k^z T^z_{ijk} \right)\exp\left(\imath \xi_i^x x\right)\exp\left(\imath \xi_j^y y\right)\exp\left(\imath \xi_k^z z\right)
\end{multline}
where $\alpha_i^x$, $\alpha_j^y$, and $\alpha_k^z$ are given by
\begin{equation}
\alpha_i^x = \frac{2\sin\left(\displaystyle \frac{\xi_i^x h_x}{2} \right)}{h_x},\quad \quad \alpha_j^y = \frac{2\sin\left(\displaystyle \frac{\xi_j^y h_y}{2} \right)}{h_y}, \quad \quad \alpha_k^z = \frac{2\sin\left(\displaystyle \frac{\xi_k^z h_z}{2} \right)}{h_z},
\end{equation}
and $\beta_i^x$, $\beta_j^y$, and $\beta_k^z$ are given by
\begin{equation}
\beta_i^x = \frac{ 1- \exp\left(-\imath \xi_i^x h_x\right) }{\imath h_x},\quad \quad \beta_j^y = \frac{ 1- \exp\left(-\imath \xi_j^y h_y\right) }{\imath h_y},\quad \quad \beta_k^z = \frac{ 1- \exp\left(-\imath \xi_k^z h_z\right) }{\imath h_z}.
\end{equation}

The relation between the Fourier series coefficients of the fluctuation and the components of the polarization then follows
\begin{equation}\label{eq:strong_fourier_series_coeff_generalized}
\Phi^*_{ijk} = - \imath \frac{\beta_i^x T^x_{ijk} + \beta_j^y T^y_{ijk} + \beta_k^z T^z_{ijk} }{c_0\left((\alpha_i^x)^2 + (\alpha_j^y)^2 + (\alpha_k^z)^2\right)}.
\end{equation}
As previously this relation holds for all frequencies, i.e.\/ for $i=-\infty,\dots,+\infty$, $j=-\infty,\dots,+\infty$ and $k=-\infty,\dots,+\infty$ (except in the particular case $i=j=k=0$ where one has ${\Phi^*}_{000}=0$), and is expressed in terms of the Fourier series coefficients which need to be evaluated. Using discrete Fourier transforms to evaluate these coefficients with $(N_x+1)\times(N_y+1)\times(N_z+1)$ voxels (so that the spatial scales are $\Delta x = L_x/(N_x+1)$, $\Delta y = L_y/(N_y+1)$ and $\Delta z = L_z/(N_z+1)$), one has
\begin{equation}\label{eq:strong_dft_generalized}
\widehat{\Phi^*}_{ijk} = - \imath \frac{\beta_i^x \widehat{T}^x_{ijk} + \beta_j^y \widehat{T}^y_{ijk} + \beta_k^z \widehat{T}^z_{ijk} }{c_0\left((\alpha_i^x)^2 + (\alpha_j^y)^2 + (\alpha_k^z)^2\right)},
\end{equation}
but now for $i=-N_x/2,\dots,N_x/2$, $j=-N_y/2,\dots,N_y/2$ and $k=-N_z/2,\dots,N_z/2$, except in the particular case $i=j=k=0$ where one has $\widehat{\Phi^*}_{000}=0$. Equation~\eqref{eq:strong_dft_generalized} provides the coefficients of the discrete Fourier transform of $\varphi^{*}$ as a function of the coefficients of the discrete Fourier transform of $\tau_x$, $\tau_y$ and $\tau_z$. It then can be used to define a \emph{generalized Green operator} on the eletric field as
\begin{align}\label{eq:strong_E_gen}
\left\{
\begin{array}{lll}
\widehat{E}^x_{ijk} & = & \displaystyle - \beta_i^{x*}\frac{\beta_i^x \widehat{T}^x_{ijk} + \beta_j^y \widehat{T}^y_{ijk} + \beta_k^z \widehat{T}^z_{ijk} }{c_0\left((\alpha_i^x)^2 + (\alpha_j^y)^2 + (\alpha_k^z)^2\right)} \\[0.5cm]
\widehat{E}^y_{ijk} & = & \displaystyle - \beta_j^{y*} \frac{\beta_i^x \widehat{T}^x_{ijk} + \beta_j^y \widehat{T}^y_{ijk} + \beta_k^z \widehat{T}^z_{ijk} }{c_0\left((\alpha_i^x)^2 + (\alpha_j^y)^2 + (\alpha_k^z)^2\right)} \\[0.5cm]
\widehat{E}^z_{ijk} & = & \displaystyle - \beta_k^{z*} \frac{\beta_i^x \widehat{T}^x_{ijk} + \beta_j^y \widehat{T}^y_{ijk} + \beta_k^z \widehat{T}^z_{ijk} }{c_0\left((\alpha_i^x)^2 + (\alpha_j^y)^2 + (\alpha_k^z)^2\right)},
\end{array}\right.
\end{align}
for $i=-N_x/2,\dots,N_x/2$, $j=-N_y/2,\dots,N_y/2$ and $k=-N_z/2,\dots,N_z/2$ (except in the case $i=j=k=0$). In the case $i=j=k=0$, one has
\begin{equation}
\widehat{E}^x_{000} = \bar{E}_x,\quad \widehat{E}^y_{000} = \bar{E}_y, \widehat{E}^z_{000} = \bar{E}_z,
\end{equation}
where $\bar{E}_x$, $\bar{E}_y$ and $\bar{E}_z$ are the components of the average of the electric field $\bar{\mathbf{E}}$. In equation~\eqref{eq:strong_E_gen}, $\beta_i^{x*}$, $\beta_j^{y*}$ and $\beta_k^{z*}$ are respectively the complex conjugate of $\beta_i^{x}$, $\beta_j^{y}$ and $\beta_k^{z}$ and $\widehat{E}^x_{ijk}$, $\widehat{E}^y_{ijk}$ and $\widehat{E}^z_{ijk}$ are the coefficients of the discrete Fourier transforms of the components of the local electric field $\mathbf{E}$. In addition, it is worth noting that
\begin{equation}
\beta_i^{x*}\beta_i^{x} = (\alpha_i^x)^2,
\end{equation}
so one finally has
\begin{equation}
\widehat{\bm{E}}(\bm{\xi}) = - \widehat{\bm{\Gamma}}^{0}_G(\bm{\xi}):\widehat{\bm{\tau}}(\bm{\xi})\quad \forall \bm{\xi}\neq\bm{0},\quad \widehat{\bm{E}}(\bm{0})=\bar{\mathbf{E}},
\end{equation}
where $\widehat{\bm{\Gamma}}^0_G$ is the generalized Green operator which has an explicit form in Fourier space:
\begin{equation}\label{eq:green_operator_generalized}
\widehat{\bm{\Gamma}}^0_G(\bm{\xi}) = \frac{\bm{\beta}^*(\bm{\xi})\otimes\bm{\beta}(\bm{\xi})}{c^0 |\bm{\beta}^*(\bm{\xi})\cdot\bm{\beta}(\bm{\xi})|^2}.
\end{equation}

Some comments are in order:
\begin{itemize}
\item In the case $h_x\to 0$, $h_y\to 0$, $h_z\to 0$, it is easy to show that $\bm{\beta}(\bm{\xi})\to \bm{\xi}$ and therefore $\widehat{\bm{\Gamma}}^0_G(\bm{\xi})\to \widehat{\bm{\Gamma}}^0(\bm{\xi})$. This result is expected as in this case, the modified deritative rules reduce to standard derivative rules, so the generalized Green operator reduces to the standard ``continuous'' Green operator.
\item In the case $h_x=\Delta x$, $h_y=\Delta y$ and $h_z=\Delta z$, the generalized Green operator reduces to the so-called ``discrete'' Green operator of~\cite{willot_fourier-based_2014}, which is also expected as the derivation rules are the same. The main difference is that~\cite{willot_fourier-based_2014}'s scheme is based on finite differences, and therefore the solution field is only defined as the nodes, while the present scheme provides the same nodal values but the solution fields are described by Fourier series and therefore are defined in the whole domain.
\item In the cases $0<h_x<\Delta x$, $0<h_y<\Delta y$ and $0<h_z<\Delta z$, the associated Green operator generalizes the ``continuous'' and ``discrete'' Green operators.
\end{itemize}

The effect of the Green operator, obtained using modified derivation rule, is now investigated on the checkerboard problem considered in Section~\ref{sec:num}. We consider the case $N=M=1024$ and several values for the spatial steps used for the derivative rules:
\begin{itemize}
\item Case 1: $h_x=h_y \to 0$, which corresponds to the (standard) continuous Green operator.
\item Case 2: $h_x=h_y=\Delta x/2=\Delta y/2$.
\item Case 3: $h_x=h_y=\Delta x=\Delta y$, which corresponds to the discrete Green operator of~\cite{willot_fourier-based_2014};
\end{itemize}
The solution field at the vicinity of the interface between the two phases is plotted in Figure~\ref{fig:Green_comparison}, on the line $x_2=0$.

\begin{figure}[!ht]
\centering
\subfloat[]{\includegraphics[width=6.5cm]{figures/line_N1024_M_varies_Green}} \\
\subfloat[]{\includegraphics[width=6.5cm]{figures/line_N1024_M_varies_Green_zoom_cont}} \hspace{0.5cm}
\subfloat[]{\includegraphics[width=6.5cm]{figures/line_N1024_M_varies_Green_zoom_disc}}
\caption{Effect of the Green operator on the normalized electric field $E_1/\VbE$. (a) Line plot at the interface between the two phases ($x_2=0$), (b) Enlarged view, (c) Representation with a linear interpolation between grid node values.}\label{fig:Green_comparison}
\end{figure}

The very first observation is that the predictions obtained with the three considered Green operators lead to the almost exact same results when the solution is represented on the fine grid (Figure~\ref{fig:Green_comparison}a). Interestingly, an enlarged view (Figure~\ref{fig:Green_comparison}b) shows that the grid nodes values calculated with the continuous and discrete Green operators are not the same, although the continuous functions (projected on the fine grid) almost coincide. It is worth noting that, if the solution is interpolated linearly between grid nodes (Figure~\ref{fig:Green_comparison}c), the oscillatory nature of the solution is very different using the continuous or discrete operator: the solution calculated with the discrete operator does not exhibit oscillations using such representation. It is important to remark that this representation (using linear interpolation) is not consistent with the framework considered in this work but it is consistent with a finite element discretization of the fields.

Therefore, the partial conclusion that can be made is that discrete Green operators do not improve intrinsically the local fields; it is the underlying discretization, based on finite elements and leading to discrete Green operators, which may improve the quality of the local fields.

\subsection{Discussion on smoothing techniques}\label{sec:smooth}

Another solution for removing oscillation is to consider a smoother. Indeed, the presentation of the results between spatial nodes clearly shows the oscillatory nature of functions described by Fourier series in presence of discontinuities, which is related to the Gibbs phenomenon. We explore in this section the effect of material properties smoothing using the periodic smoother proposed by~\cite{morin_periodic_2021}. This smoother is based on a modification of the Fourier coefficients of the function to be smoothed, which is briefly recalled hereafter in 1-d. We consider periodic unidimensional function $f$ with period $L$, defined in the interval $[0,~L] $ discretized with $N+1$ segments so the spatial scale is $\Delta x = L/(N+1)$. The discrete values of function $f$ denoted $f_j= f(x_j)$ may thus be written as an array $\bm{f}$ of size $N$. The smooth estimate is denoted by $\tilde{f}$, and we are thus looking for its discrete values $\tilde{f}_i$ which consists in finding the array $\tilde{\bm{f}}=\left(\tilde{f}_0,\tilde{f}_1,\dots,\tilde{f}_{N}\right)$. The smoother proposed by~\cite{morin_periodic_2021} is based on smoothing splines solution of the minimization of a functional which balances the fidelity to the data, through the residual sum-of-squares (RSS), and the smoothness of the estimate, through some penalty term. In the particular case of the first derivative considered in the penalty term, the smoother reads
\begin{equation}\label{eq:smoother}
\tilde{\bm{f}} = \IDFT \left(\bm{\Gamma}(\bs) \circ \DFT (\bm{f}) \right),
\end{equation}
where $\circ$ denoted Hadamard product (pointwise product) and $\bm{\Gamma}(\bs)$ is an array of size $N$, whose components are given by
\begin{equation}
\Gamma_i = \frac{1}{1+ \bs\,L^{2} \xi_i^{2}}.
\end{equation}
In this equation, $\bs$ is a parameter that controls the smoothing, i.e.\/ the tradeoff between the ``smoothness'' of the solution and the fidelity to the data. A given value of $\bs$ will create some interphase between material discontinuities: the increase of $\bs$ increases the thickness of this interphase.

The objective is now to study the effect of this smoother on the local fields, especially on the occurrence of oscillations. Several spatial discretizations are considered; in order to assess the ability of the smoother to reproduce and get closer to the analytical solution, the thickness of the interphase is chosen so as to decrease with the spatial discretization. We thus take the following form for the smoothing parameter
\begin{equation}
\bs = \frac{a}{N^2}.
\end{equation}
This allows to have a finite interphase for small values of $N$ and the size of this interphase tends to zero when the spatial discretization increases. The microstructure obtained increasing the number of nodes $N$ tends to the initial microstructure (with discontinuities): the local fields are thus expected to get closer to the analytical solution, when $N$ increases, without numerical artifacts. Several values for the parameter $a$ will be consider.

We begin with the case $\bs = {1}/{N^2}$ (i.e.\/ $a=1$). First, the microstructure is represented in Figure~\ref{fig:smooth_comparisons}a for several values of the discretization (and thus the smoothing parameter); the node values are represented by discrete symbols and the associated continuous representation is also shown (based on a Fourier interpolation using the node values to calculate the Fourier coefficients and then the Fourier coefficients to reconstruct the function everywhere). Since the smoothing parameter has been taken to decrease when the number of nodes $N$ increases, the thickness of the interphase decreases when $N$ increases and the smooth microstructure tends to the theoretical one.

The line plot of the solution field $E_1/\VbE$ at the interface between the two phases ($x_2=0$) is then represented in Figure~\ref{fig:smooth_comparisons}b. As expected by the effect of smoothing, (i) the solution is spread due to the presence of an interphase of properties and (ii) spurious oscillations are considerably decreased by comparison with the solution calculated without smoothing (see e.g. Figure~\ref{fig:line_M128_N_varies}). The smoothing parameter has been chosen so as to remove most of the oscillations without introducing too much spreading, for a given value of $N$; an increase of the smoothing parameter would necessarily decrease the magnitude of the oscillations but increase the spreading of the solution: a compromise between oscillations and spreading is inevitable. It is worth noting that the numerical solution seems to tend to the analytical one, without numerical artifacts, with an increase of the spatial resolution as well as the introduction of a sufficient (but small enough) smoothing parameter.


\begin{figure}[!ht]
\centering
\subfloat[]{\includegraphics[width=6.5cm]{figures/line_c_smooth}} \hspace{0.5cm}
\subfloat[]{\includegraphics[width=6.5cm]{figures/line_E1_smooth}}
\caption{Effect of smoothing on the local fields. (a) Line plot of the conductivity field $c$ at the center of the cell ($x_2=0.5$). (b) Line plot of the normalized electric field $E_1/\VbE$ at the interface between the two phases ($x_2=0$). The continuous lines corresponds to the solution on the fine grid and the triangles corresponds to the solution on the grid used for the FFT calculations.}\label{fig:smooth_comparisons}
\end{figure}

A quantitative comparison between the numerical and analytical solutions is then made but studying the error defined by equation~\eqref{eq:err}, for $\bs = {1}/{N^2}$ and $\bs = 1/{(200\times N^2)}$. The error is represented in Figure~\ref{fig:err_smooth}a as a function of the number of modes $M$ (for a given number of nodes $N$) and is compared to that obtained without smooth. In contrast to the original case where the error starts increasing for high frequencies (when $M > M_{\rm lim}$), the error decreases and reaches a stationary value for a smooth microstructure. This confirms that oscillations due to high frequencies do not degrade the solution field because they are greatly decreased in the smooth case. It is worth noting that the error, in this case, is greater than the error of the non-smooth case. This is due to the fact that the smooth microstructure is not strictly equivalent to the non-smooth case because of the thin interface induced by smoothing; therefore the problem solved is slightly different which necessarily induces errors in comparison with the analytical solution. {Moreover, this result was expected as the non-smoothed solution, without integration error ($N\gg M$), is optimal in the sense of the error defined by~\eqref{eq:err}. In a way, the non-smoothed solution optimized this error while disregarding the local error, therefore the spurious oscillation.} This is confirmed by considering a smaller value for the smoothing parameter, $\bs = 1/{(200\times N^2)}$. In this case, the error is represented in Figure~\ref{fig:err_smooth}b, again as a function of the number of modes $M$ (for a given number of nodes $N$). For this smaller value of the smoothing parameter, the error keeps decreasing when increasing the number of modes and is smaller than in the case without smoothing. This indicates that decreasing the interphase improves the local fields because (i) the microstructure is closer to the analytical one and (ii) oscillations are still removed by the smoother. It should be noted that, compared to the truncation of high frequencies considered previously (i.e.\/ for $M \sim M_{\rm lim}$), the error defined by~\eqref{eq:err} is not significantly improved by the smoother; this can be explained by the fact that filtering high frequencies and smoothing material discontinuities play somehow a similar role, as both acts as a regularization of spurious oscillations due to high frequencies in presence of discontinuities.


\begin{figure}[!ht]
\centering
\subfloat[]{\includegraphics[width=6.25cm]{figures/error_N_smooth1}}\hspace{0.75cm}
\subfloat[]{\includegraphics[width=6.25cm]{figures/error_N_smooth2}}
\caption{Influence of material properties smoothing on error analysis {on the local fields}. (a) Case $\bs = {1}/{N^2}$ and (b) Case $\bs = 1/{(200\times N^2)}$.}\label{fig:err_smooth}
\end{figure}

For completeness, the influence of the smoother on the macroscopic conductivity is represented in Figure~\ref{fig:err_smooth_macro}. As for the error on the local fields, the error on the macroscopic conductivity decreases with the increase of $M$ and reaches a stationary value for a smooth microstructure. Overall, independently of the value of the smoothing parameter, the use of the smoother slightly alters the predicted effective properties. This behaviour is attributed to the modification of the microstructure induced by the smoothing operation, meaning that the resulting computation corresponds to a slightly different physical problem.

\begin{figure}[!ht]
\centering
\subfloat[]{\includegraphics[width=6.25cm]{figures/error_macro_N_smooth1}}\hspace{0.75cm}
\subfloat[]{\includegraphics[width=6.25cm]{figures/error_macro_N_smooth2}}
\caption{Influence of material properties smoothing on error analysis on the macroscopic conductivity. (a) Case $\bs = {1}/{N^2}$ and (b) Case $\bs = 1/{(200\times N^2)}$.}\label{fig:err_smooth_macro}
\end{figure}


\section{Conclusion}

The aim of this work was to revisit the formulation of FFT-based methods for heterogeneous materials in the periodic setting, originally proposed by the pioneering work of~\cite{moulinec_numerical_1998}. This class of methods is based on the iterative resolution of an auxiliary problem involving a reference homogeneous material and a polarization tensor. Assuming a description of the fields by (infinite or partial) Fourier series, we show the equivalence between { several discretization approaches including strong and weak formulations and different numerical approximations}, following ideas from~\cite{vondrejc_fft-based_2014}. The formulation proposed in this work notably decouples the number of Fourier modes to the spatial discretization grid which is shown to be used to calculate approximately integrals with the trapezoidal rule. The relation between the discretization grid and the representation of the continuous fields is discussed including the possibility of using non-uniform grids. Several numerical experiments are conducted on the checkerboard problem for which an analytical solution is used; we notably investigate the effect of the number of Fourier modes together with the grid discretization. The importance of representing the fields between grid nodes is emphasized as spurious oscillations can emerged between grid nodes. Increasing the number of grid nodes, for a given number of Fourier modes, has only a moderate effect as it simply improves the calculations of the Fourier coefficients. Conversely, {for this specific checkerboard problem,} increasing the number of Fourier modes, for a given grid discretization, improves the solution up to a limit values for the number of modes, above which the quality of the solution is degraded due to the occurrence of spurious oscillations at high frequencies. In this continuous setting, i.e.\/ by considering Fourier series to describe fields, generalized Green operators are derived using modified derivative rules, which leads to a generalization of both continuous and discrete Green operators relying on a single spatial parameter. It was shown that this does not improve the local fields, as oscillations are mainly due to the description of solution fields using Fourier series. On the other hand, the effect of material discontinuity smoothing is studied by using~\cite{morin_periodic_2021}'s interface smoother: it is shown that a combination of mesh refinement and suitable smoothing improves the numerical solution without numerical artifacts even between grid nodes. {Overall, the standard choice $M=N$ is not necessarily optimal.}

It is important to remark that this work was concerned with a description of fields based on Fourier series, even when (generalized) discrete Green operators were used; this is mainly why spurious oscillations between nodes were observed in that case. Standard approaches making used of discrete Green operators are generally based either on finite differences~\cite{willot_fourier-based_2014,willot2015fourier} or a finite element discretization, generally using linear elements~\cite{zeman2017finite,schneider_fft-based_2017}. Therefore in that case, the obtained fields should be post-treated according to the interpolation basis functions used and may not show any oscillations between the node grid used for evaluating the Fourier transforms as shown using a linear interpolation between grid node values. Oscillations are thus expected to be less important using these approaches because of the representation of the solution field.

This work provides a new perspective on solution fields obtained by FFT-based solvers, especially on the quality of the local fields. In presence of material discontinuities, spurious oscillations inherent of the Fourier framework degrades the solution field. Simple remedies consist in (i) filtering the high frequencies, by simply removing some Fourier coefficients as a post-treatment to reconstruct the solution field or (ii) removing material discontinuities using an interface spreading smoother.



\printCOI

\printbibliography
\end{document}