# A discrete sine-cosine based method for the elasticity of heterogeneous materials with arbitrary boundary conditions

joseph paux, Léo Morin, Lionel Gélébart, Abdoul Magid Amadou Sanoko

► To cite this version:

joseph paux, Léo Morin, Lionel Gélébart, Abdoul Magid Amadou Sanoko. A discrete sine-cosine based method for the elasticity of heterogeneous materials with arbitrary boundary conditions. Computer Methods in Applied Mechanics and Engineering, 2025, 433, pp.117488. 10.1016/j.cma.2024.117488. hal-04946333

HAL Id: hal-04946333

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

Submitted on 13 Feb 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Archive Ouverte - Open Repository

Science Arts & Métiers (SAM)

is an open access repository that collects the work of Arts et Métiers Institute of
Technology researchers and makes it freely available over the web where possible.

This is an author-deposited version published in: https://sam.ensam.eu
Handle ID: http://hdl.handle.net/10985/25775

To cite this version :

Paux JOSEPH, Léo MORIN, Lionel GELEBART, Abdoul Magid AMADOU SANOKO - A discrete sine-cosine based method for the elasticity of heterogeneous materials with arbitrary boundary conditions - 2024
<figure><img src="image_3.png" /></figure>

ELSEVIER

Contents lists available at [ScienceDirect](https://www.sciencedirect.com)

Comput. Methods Appl. Mech. Engrg.

journal homepage: [www.elsevier.com/locate/cma](https://www.elsevier.com/locate/cma)

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

# A discrete sine-cosine based method for the elasticity of heterogeneous materials with arbitrary boundary conditions

Joseph Paux<sup>a,b</sup>, Léo Morin<sup>c,d,*</sup>, Lionel Gélébart<sup>e</sup>, Abdoul Magid Amadou Sanoko<sup>c,d</sup>

## ARTICLE INFO

**Keywords:**

FFT-based solvers

Dirichlet/Neumann boundary conditions

Discrete sine-cosine transforms

Elasticity

## ABSTRACT

The aim of this article is to extend Moulinec and Suquet (1998)'s FFT-based method for heterogeneous elasticity to non-periodic Dirichlet/Neumann boundary conditions. The method is based on a decomposition of the displacement into a known term verifying the boundary conditions and a fluctuation term, with no contribution on the boundary, and described by appropriate sine-cosine series. A modified auxiliary problem involving a polarization tensor is solved within a Galerkin-based method, using an approximation space spanned by sine-cosine series. The elementary integrals emerging from the weak formulation of the equilibrium are approximated by discrete sine-cosine transforms, which makes the method relying on the numerical complexity of Fourier transforms. The method is finally assessed in several problems including kinematic uniform, static uniform and arbitrary Dirichlet/Neumann boundary conditions.

## 1. Introduction

This work is concerned with the development of efficient numerical methods for solving problems of micromechanics in heterogeneous materials. Fast Fourier transforms (FFT) based methods, introduced in the seminal paper of Moulinec and Suquet [1], are a class of solvers that are commonly used in micromechanics to compute the local and overall fields of heterogeneous (composite) materials and constitute an alternative to the finite element method (FEM) in problems governed by elliptic equations. The three main advantages of FFT-based solvers are that (i) the numerical complexity of the method is in $O(N \log N)$ (and is therefore lower than the complexity of matrix-assembled methods), (ii) it does not require meshing operations since it is based on voxelized images (that can be directly used from imagery techniques such as scanning electron microscopy (SEM) or tomography), and (iii) it is adapted to parallel implementation as most FFT packages are designed for distributed-memory parallel machines [2]. The method proposed by Moulinec and Suquet [1] concerns periodic microstructures and is based on an iterative solution to an integral Lippmann-Schwinger type equation, which can be solved using a fixed-point method (as done by Moulinec and Suquet [1] and referred as the so-called basic scheme) or accelerated iterative schemes [3-7], which are more suitable for highly-contrasted materials. In the context of periodic homogenization, the method has been successfully applied to a large class of (non-)linear

<footnote><sup>a</sup> Institut de Thermique, Mécanique, Matériaux (ITheMM, EA 7548), Université de Reims Champagne-Ardenne, Campus Sup Ardenne, 08000 Charleville-Mézières, France</footnote>

<footnote><sup>b</sup> Université Sorbonne Paris Nord, Laboratoire des Sciences des Procédés et des Matériaux, LSPM, CNRS, UPR 3407, F-93430, Villetaneuse, France</footnote>

<footnote><sup>c</sup> Univ. Bordeaux, CNRS, Bordeaux INP, I2M, UMR 5295, F-33400 Talence, France</footnote>

<footnote><sup>d</sup> Arts et Metiers Institute of Technology, CNRS, Bordeaux INP, I2M, UMR 5295, F-33400 Talence, France</footnote>

<footnote><sup>e</sup> Université Paris-Saclay, CEA, SRMA, 91191, Gif/Yvette, France</footnote>

<footnote>* Corresponding author at: Univ. Bordeaux, CNRS, Bordeaux INP, I2M, UMR 5295, F-33400 Talence, France.
E-mail address: leo.morin@u-bordeaux.fr (L. Morin).</footnote>
problems, including, among others, anisotropic elasticity [1], J2-plasticity [1], crystal viscoplasticity [8], dislocation-mediated plasticity [9,10], electrical conduction [3], piezoelectricity [11], porous ductile solids [12,13] (see the papers of Schneider [14], Lucarini et al. [15] for a comprehensive review of FFT-based method applied to homogenization).

Periodic boundary conditions are intrinsic of the method as it is based on a description of the fields through Fourier series. As a consequence, this method is not adapted in several problems of interest in micromechanics such as (i) materials with non-periodic microstructures (such as fibrous networks and particular situations of polycrystalline materials) which requires to periodicize the microstructures obtained from imagery techniques, (ii) higher-order homogenization<fnref n="1" /> with higher-order boundary conditions [18,19], and (iii) transient problems such as elastodynamic with wave propagation [20] for which periodic boundary conditions induce two traveling waves.

The extension of FFT-based solvers to non-periodic boundary conditions has been first proposed in the work of Gélébart [21], devoted to elasticity problems subjected to Dirichlet boundary conditions, through the introduction of a buffer zone surrounding the domain of interest. The complete cell, composed of the domain of interest plus the buffer zone, can be periodic so this technique is based on standard FFT packages with periodic boundary conditions applied on the exterior of the buffer zone. Despite its ability to apply correctly Dirichlet and also Neumann (see [22]) boundary conditions, this method requires (i) a larger unit cell than the microstructure studied because of the (elastic) buffer zone and (ii) additional internal iterations to find the eig displacement field at the boundary of the buffer zone. Mention has to be made to the works of Grimm-Strele and Kabel [23] and Monchiet and Bonnet [24] which are based on mirror unit cells and can be used for *uniform* Dirichlet and/or Neumann boundary conditions. This method, initially proposed by Wiegmann [25], takes advantage of the equivalence between particular symmetries of the unit cell and boundary conditions on a portion of the unit cell, and can be related to discrete sine-cosine transforms.

Recently, an alternative method has been proposed for imposing non-periodic boundary conditions using explicitly discrete sine-cosine transforms for uniform [26] or arbitrary boundary conditions [27–29], for problems of heterogeneous conductivity. The idea of this approach is to split the solution field (e.g. the temperature field in thermal conductivity problems) of the elliptic problem into a known term verifying the Dirichlet/Neumann boundary conditions and an unknown fluctuation term described by an appropriate trigonometric series with no contribution on the boundary. The work by Paux et al. [29] concerns the conductivity of heterogeneous materials subjected to mixed Dirichlet/Neumann boundary conditions<fnref n="2" /> and is based on a Galerkin discretization of the auxiliary problem, using sine-cosine series as trial functions. The work of Gélébart [27], on the other hand, concerns the conductivity of heterogeneous materials subjected to mixed Dirichlet/Neumann (and also periodic) boundary conditions and is based on a finite difference scheme. Both approaches (Galerkin and finite differences based) rely on the use of discrete sine and cosine transforms which can be computed using FFT packages and have thus a computational complexity in $O(N \log N)$. It has been shown that both methods lead to very similar results in terms of local and overall fields but they differ in terms of convergence; the convergence rate of the Galerkin-based schemes scales linearly with the contrast, while the number of iterations to convergence is bounded or slowly increases with the contrast using a discrete scheme [29] (see also [30]). Mention has to be made to the recent work of Risthaus and Schneider [31], who extends the approaches for conductivity [26,28] to elasticity problems with kinematic uniform boundary conditions using sine series (as a particular case of Dirichlet boundary conditions). It must be noted that discrete sine-cosine transforms have also been used in the literature to solve efficiently the Poisson equation (discretized by a finite difference scheme) instead of using matrix-inversion (see e.g. [32–34]).

The aim of this work is to develop a numerical method based on discrete sine-cosine transforms for heterogeneous elasticity problems subjected to general non-periodic mixed Dirichlet/Neumann boundary conditions (including, as subcases, kinematic uniform boundary conditions (KUBC) and static uniform boundary conditions (SUBC)). Thus this work extends the approach of Morin and Paux [28], Gélébart [27], Paux et al. [29], restricted to conductivity problems, to elasticity. The main difference between elasticity and conductivity equations is that elasticity involves higher-order tensors which will induce extra difficulties in the resolution of the auxiliary problem, mainly due to the presence of cross derivatives in the equilibrium equations (see e.g. [25,26]). The present method thus constitutes a direct extension of Moulinec and Suquet [1]'s method to general non-periodic boundary conditions.

The paper is organized as follows. In Section 2, important results related to sine-cosine series and their associated transforms are recalled. The numerical method for elasticity problems subjected to Dirichlet/Neumann boundary conditions is then presented in Section 3. Finally, the numerical scheme is applied in Section 4 to several problems of heterogeneous elasticity including homogenization under kinematic uniform or static uniform boundary conditions.

## 2. Sine-cosine series and associated discrete transforms

The resolution of the heterogeneous elasticity problem with arbitrary Dirichlet/Neumann boundary conditions will be based, as in conductivity problems [27–29], on a decomposition of the solution field in a known (given) field verifying the boundary conditions and a fluctuation field having no contribution on the boundary conditions. The approach followed in this work is based on a Galerkin method, which requires the definition of a continuous trial function for the unknown field, described by half-range or quarter-range sine-cosine series depending on the boundary conditions. Moreover, as in the periodic case [1,35], the elementary

<footnote n="1">¹ Higher-order homogenization can be handled in the periodic setting in several cases, as done by Gélébart [16] by adding strain gradient components for plates of beams or by Tran et al. [17] using the asymptotic expansion method.</footnote>

<footnote n="2">² The work of Paux et al. [29] extends Morin and Paux [28]’s approach, which was restricted to Dirichlet boundary conditions, to mixed Dirichlet/Neumann boundary conditions.</footnote>
<table><caption>**Table 1**
Functions and coefficients defining the four mixed Dirichlet/Neumann boundary conditions on the faces $x = 0$ and $x = L_x$.</caption><thead><tr><th></th><th>Boundary conditions</th><th>$g$ and $\bar{g}$</th><th>$k_i$</th><th>$\alpha_i$</th><th>$\xi_i$</th></tr></thead><tbody><tr><td rowspan="2">DD</td><td>$f(x, y, z) = 0 \quad \forall x = 0$</td><td>$g(x) = \sin(x)$</td><td rowspan="2">$k_i = \frac{i}{L_x} \pi$</td><td rowspan="2">$\alpha_i = 1 \quad \forall i \ge 0$</td><td rowspan="2">$\xi_i = k_i$</td></tr><tr><td>$f(x, y, z) = 0 \quad \forall x = L_x$</td><td>$\bar{g}(x) = \cos(x)$</td></tr><tr><td rowspan="2">NN</td><td>$\frac{\partial f}{\partial x}(x, y, z) = 0 \quad \forall x = 0$</td><td>$g(x) = \cos(x)$</td><td rowspan="2">$k_i = \frac{i}{L_x} \pi$</td><td>$\alpha_0 = 1/2$</td><td rowspan="2">$\xi_i = -k_i$</td></tr><tr><td>$\frac{\partial f}{\partial x}(x, y, z) = 0 \quad \forall x = L_x$</td><td>$\bar{g}(x) = \sin(x)$</td><td>$\alpha_i = 1 \quad \forall i \ge 1$</td></tr><tr><td rowspan="2">DN</td><td>$f(x, y, z) = 0 \quad \forall x = 0$</td><td>$g(x) = \sin(x)$</td><td rowspan="2">$k_i = \frac{2i+1}{2L_x} \pi$</td><td rowspan="2">$\alpha_i = 1 \quad \forall i \ge 0$</td><td rowspan="2">$\xi_i = k_i$</td></tr><tr><td>$\frac{\partial f}{\partial x}(x, y, z) = 0 \quad \forall x = L_x$</td><td>$\bar{g}(x) = \cos(x)$</td></tr><tr><td rowspan="2">ND</td><td>$\frac{\partial f}{\partial x}(x, y, z) = 0 \quad \forall x = 0$</td><td>$g(x) = \cos(x)$</td><td rowspan="2">$k_i = \frac{2i+1}{2L_x} \pi$</td><td rowspan="2">$\alpha_i = 1 \quad \forall i \ge 0$</td><td rowspan="2">$\xi_i = -k_i$</td></tr><tr><td>$f(x, y, z) = 0 \quad \forall x = L_x$</td><td>$\bar{g}(x) = \sin(x)$</td></tr></tbody></table>

integrals arising in the weak formulation will be approximated using discrete sine-cosine transforms, which are intimately related to the sine-cosine series [36,37]. It must be noted that discrete sine transforms (DST) and discrete cosine transforms (DCT) are related to standard discrete Fourier transforms (DFT) through symmetry extensions; a given non-periodic signal can be periodized by symmetry extension (that can be periodic or anti-periodic, and located an ending point or a half-element after) and the DFT of this periodized signal corresponds to a DST or DCT depending on the symmetry extension (see [27]). We thus recall, following the presentation of Paux et al. [29], important results related to sine and cosine series and their associated discrete transforms in a generic framework.

## 2.1. 1-d case

### 2.1.1. Sine-cosine series in 1-d

We consider a function $f$ defined in the interval $[0, L_x]$ and verifying either null Dirichlet or null Neumann boundary conditions. Such function can be described by half-range or quarter-range sine-cosine series and written under the generic form

$$f(x) = \sum_{i=0}^{+\infty} \alpha_i F_i g(k_i x), \quad \forall x \in [0, L_x], \qquad (1)$$

with $k_i$ some “frequency” parameter, $g$ a cosine or a sine function (depending on the type of boundary conditions) and $F_i$, $i \ge 0$, the sine-cosine series coefficients given by

$$F_i = \frac{2}{L_x} \int_{0}^{L_x} f(x)g(k_i x)dx, \quad \forall i \ge 0. \qquad (2)$$

The associated partial series of order $N$ (for every positive integer $N$) is then defined as

$$f_N(x) = \sum_{i=0}^{N} F_i g(k_i x), \quad \forall x \in [0, L_x]. \qquad (3)$$

A summary of all functions and coefficients is given in Table 1 to cover all types of boundary conditions on $x = 0$ and $x = L_x$ (Dirichlet-Dirichlet (DD), Neumann-Neumann (NN), Dirichlet-Neumann (DN) and Neumann-Dirichlet (ND)).

Derivative rules associated with (1) can be written as

$$\frac{d f}{d x}(x) = \sum_{i=0}^{+\infty} \xi_i \alpha_i F_i \bar{g}(k_i x), \quad \forall x \in [0, L_x], \qquad (4)$$

$$\frac{d^2 f}{d x^2}(x) = \sum_{i=0}^{+\infty} -(\xi_i)^2 \alpha_i F_i g(k_i x), \quad \forall x \in [0, L_x]. \qquad (5)$$

In these expressions, $\bar{g}$ is the dual function of $g$ (i.e. a sine function if the initial function is a cosine function and vice versa) and $\xi_i$ is the coefficient arising from the derivation (see Table 1).

## 2.2. Discrete sine-cosine transforms in 1-d

A function that is described by a sine-cosine series given by Eq. (1) relies on the computation of the associated sine-cosine coefficients given in (2). For an arbitrary function $f$, there is no analytical formula for calculating the integral in Eq. (2). The discrete sine transform (DST) and the discrete cosine transform (DCT) provide the numerical integration<fnref n="3" /> of the coefficients defined by Eq. (2), depending on the type of boundary conditions. An efficient computation of the coefficients of the DFT, DST and DCT is classically done using fast Fourier transforms (FFT) [38].

<footnote n="3">³ This numerical integration is similar to the standard Discrete Fourier Transform (DFT) used for the computation of the coefficients of the Fourier series for periodic functions.</footnote>
<table><caption>**Table 2**

Computational details of the Discrete sine–cosine Transforms.</caption><thead><tr><th>BC</th><th>Transform</th><th>Î<sub>i</sub> (i = 0, ..., N + 1)</th><th>f<sub>a</sub> (a = 0, ..., N + 1)</th></tr></thead><tbody><tr><td rowspan="3">DD</td><td>DST-I</td><td>Î<sub>i</sub> = &sum;<sub>a=0</sub><sup>N+1</sup> &beta;<sub>a</sub>f<sub>a</sub> sin(&pi;a<i>i</i> / (N + 1))</td><td>f<sub>a</sub> = 2 / (N + 1) &sum;<sub>i=0</sub><sup>N+1</sup> &gamma;<sub>i</sub>Î<sub>i</sub> sin(&pi;a<i>i</i> / (N + 1))</td></tr><tr><td>Î&#x0302; = D<sub>xI</sub>(f)</td><td>&beta;<sub>0</sub> = &beta;<sub>N+1</sub> = 1/2</td><td>&gamma;<sub>0</sub> = &gamma;<sub>N+1</sub> = 1/2</td></tr><tr><td>f = D<sub>sI</sub><sup>-1</sup>(Î&#x0302;)</td><td>&beta;<sub>a</sub> = 1 for a = 1, ..., N</td><td>&gamma;<sub>i</sub> = 1 for i = 1, ..., N</td></tr><tr><td rowspan="3">NN</td><td>DCT-I</td><td>Î&#x0302; = &sum;<sub>a=0</sub><sup>N+1</sup> &beta;<sub>a</sub>f<sub>a</sub> cos(&pi;a<i>i</i> / (N + 1))</td><td>f<sub>a</sub> = 2 / (N + 1) &sum;<sub>i=0</sub><sup>N+1</sup> &gamma;<sub>i</sub>Î&#x0302; cos(&pi;a<i>i</i> / (N + 1))</td></tr><tr><td>Î&#x0302; = D<sub>cI</sub>(f)</td><td>&beta;<sub>0</sub> = &beta;<sub>N+1</sub> = 1/2</td><td>&gamma;<sub>0</sub> = &gamma;<sub>N+1</sub> = 1/2</td></tr><tr><td>f = D<sub>cI</sub><sup>-1</sup>(Î&#x0302;)</td><td>&beta;<sub>a</sub> = 1 for a = 1, ..., N</td><td>&gamma;<sub>i</sub> = 1 for i = 1, ..., N</td></tr><tr><td rowspan="3">DN</td><td>DST-III</td><td>Î<sub>i</sub> = &sum;<sub>a=0</sub><sup>N+1</sup> &beta;<sub>a</sub>f<sub>a</sub> sin(&pi;a(i + 1/2) / (N + 1))</td><td>f<sub>a</sub> = 2 / (N + 1) &sum;<sub>i=0</sub><sup>N+1</sup> &gamma;<sub>i</sub>Î<sub>i</sub> sin(&pi;a(i + 1/2) / (N + 1))</td></tr><tr><td>Î&#x0302; = D<sub>sIII</sub>(f)</td><td>&beta;<sub>0</sub> = &beta;<sub>N+1</sub> = 1/2</td><td>&gamma;<sub>i</sub> = 1 for i = 0, ..., N</td></tr><tr><td>f = D<sub>sIII</sub><sup>-1</sup>(Î&#x0302;)</td><td>&beta;<sub>a</sub> = 1 for a = 1, ..., N</td><td>&gamma;<sub>N+1</sub> = 0</td></tr><tr><td rowspan="3">ND</td><td>DCT-III</td><td>Î<sub>i</sub> = &sum;<sub>a=0</sub><sup>N+1</sup> &beta;<sub>a</sub>f<sub>a</sub> cos(&pi;a(i + 1/2) / (N + 1))</td><td>f<sub>a</sub> = 2 / (N + 1) &sum;<sub>i=0</sub><sup>N+1</sup> &gamma;<sub>i</sub>Î<sub>i</sub> cos(&pi;a(i + 1/2) / (N + 1))</td></tr><tr><td>Î&#x0302; = D<sub>cIII</sub>(f)</td><td>&beta;<sub>0</sub> = &beta;<sub>N+1</sub> = 1/2</td><td>&gamma;<sub>i</sub> = 1 for i = 0, ..., N</td></tr><tr><td>f = D<sub>cIII</sub><sup>-1</sup>(Î&#x0302;)</td><td>&beta;<sub>a</sub> = 1 for a = 1, ..., N</td><td>&gamma;<sub>N+1</sub> = 0</td></tr></tbody></table>

We consider the 1-d domain of size [0, L<sub>x</sub>], discretized uniformly with N<sub>x</sub>+2 points. The spatial scale associated with the uniform grid is thus Δx = L<sub>x</sub>/(N<sub>x</sub> + 1) we denote by x<sub>a</sub> = aΔx (for a = 0, ..., N + 1) the grid points. The grid points values of function f are denoted by f<sub>a</sub> = f(x<sub>a</sub>) (for a = 0, ..., N + 1). The discrete transform associated to each boundary conditions<sup>4</sup> can be written under the generic form as follows

$$ \hat{F} = D_x(f), \qquad (6) $$

where f is the array of size N + 2 containing the values of function f at the grid points x<sub>a</sub> and F is the array of size N + 2 containing the values of the associated discrete transform coefficients F<sub>i</sub>. The notation D<sub>x</sub>(·) corresponds to the discrete transform operator (with x = sI for the DST-I, x = cI for the DCT-I, x = sIII for the DST-III and x = cIII for the DCT-III, depending on the boundary conditions, see footnote 3). The inverse discrete transform can be formally written as

$$ f = D_x^{-1}(\hat{F}), \qquad (7) $$

where the notation D<sub>x</sub><sup>-1</sup>(·) corresponds to the associated inverse discrete transform operator. In addition, the inverse transforms of type I (sine and cosine) discrete transforms are type I (sine and cosine) discrete transforms multiplied by a factor 2/(N + 1), while the inverse transforms of type III (sine and cosine) discrete transforms are type II (sine and cosine) discrete transforms multiplied by a factor 2/(N + 1).

The expressions of all discrete transforms is given in Table 2. Interestingly, one can notice that the coefficients F<sub>i</sub> provide a numerical approximation of the coefficients F<sub>i</sub> (of the associated cosine or sine series) through

$$ \hat{F}_i \approx \frac{N+1}{2} F_i, \quad \forall i \ge 0. \qquad (8) $$

Finally, it is interesting to note that the four discrete transforms can be put under the following generic form

$$ \hat{F}_i = \sum_{a=0}^{N+1} \beta_a f_a g(k_i x_a), \quad \forall i \ge 0, \qquad (9) $$

where β<sub>a</sub>, k<sub>i</sub>, α<sub>i</sub> and g(k<sub>i</sub>x) depend on the boundary conditions. In addition, the inverse transforms can be written as

$$ f_a = \sum_{i=0}^{N+1} \gamma_i \hat{F}_i g(k_a x_i), \quad \forall a \ge 0. \qquad (10) $$

## 2.3. The 3-d case

### 2.3.1. 3-d sine–cosine series

We consider a function f defined in the prismatic domain Ω = [0, L<sub>x</sub>] × [0, L<sub>y</sub>] × [0, L<sub>z</sub>] and verifying null Dirichlet boundary condition on faces ∂Ω<sub>D</sub>

$$ f(x, y, z) = 0, \quad \forall (x, y, z) \in \partial\Omega_D, \qquad (11) $$

<footnote n="4">4 DST-I encodes Dirichlet–Dirichlet, DST-I encodes Neumann–Neumann, DST-III encodes Dirichlet–Neumann and DCT-III encodes Neumann–Dirichlet.</footnote>
and null Neumann boundary condition on faces $\partial\Omega_N$

$$ \nabla f(x, y, z) \cdot \mathbf{n} = 0, \quad \forall(x, y, z) \in \partial\Omega_N, \qquad (12) $$

where $\mathbf{n}$ is the outward normal to the boundary $\partial\Omega_N$ (with $\partial\Omega_N \cup \partial\Omega_D = \partial\Omega$). Function $f$ can be described by a generalized 3-d sine-cosine series as

$$ f(x, y, z) = \sum_{i=0}^{+\infty} \sum_{j=0}^{+\infty} \sum_{k=0}^{+\infty} \alpha_i^x \alpha_j^y \alpha_k^z F_{ijk} g^x(k_i^x x) g^y(k_j^y y) g^z(k_k^z z), \quad \forall (x, y, z) \in \Omega, \qquad (13) $$

where the functions $g^x$, $g^y$, $g^z$, and the coefficients $\alpha_i^x$, $\alpha_j^y$, $\alpha_k^z$, $k_i^x$, $k_j^y$, $k_k^z$ depend on the type of boundary conditions on each couple of opposite faces<fnref n="5" /> (see Table 1), and $F_{ijk}$ are the sine-cosine series coefficients, given by

$$ F_{ijk} = \frac{8}{L_x L_y L_z} \int_{\Omega} f(x, y, z) g^x(k_i^x x) g^y(k_j^y y) g^z(k_k^z z) \, dx dy dz, \quad \forall i, j, k \ge 0. \qquad (14) $$

The partial derivatives of $f$ read

$$ \frac{\partial f}{\partial x}(x, y, z) = \sum_{i=0}^{+\infty} \sum_{j=0}^{+\infty} \sum_{k=0}^{+\infty} \overline{g^x(k_i^x x) g^y(k_j^y y) g^z(k_k^z z)}, \quad \forall (x, y, z) \in \Omega, \qquad (15) $$

$$ \frac{\partial^2 f}{\partial x^2}(x, y, z) = \sum_{i=0}^{+\infty} \sum_{j=0}^{+\infty} \sum_{k=0}^{+\infty} -\overline{(g^x)^2 \alpha_i^x \alpha_j^y \alpha_k^z F_{ijk} g^x(k_i^x x) g^y(k_j^y y) g^z(k_k^z z)}, \quad \forall (x, y, z) \in \Omega, \qquad (16) $$

where $\overline{g^x}$ is the dual function of $g^x$ (as defined in Table 1).

### 2.3.2. Discrete sine-cosine transforms in 3-d

The 3-d prismatic domain $\Omega$ is discretized uniformly with $(N_x + 2) \times (N_y + 2) \times (N_z + 2)$ points 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 + 1$), $y_b = b\Delta y$ (for $b = 0, \dots, N_y + 1$) and $z_c = c\Delta z$ (for $c = 0, \dots, N_z + 1$). 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 + 1$, $b = 0, \dots, N_y + 1$ and $c = 0, \dots, N_z + 1$).

Discrete sine-cosine transforms in 3-d are simply a composition of 1-d transforms in each spatial directions. The coefficients of the discrete sine-cosine transform ($\hat{F}_{ijk}$) of function $f$ are then given by (following the conventions of Table 2)

$$ \hat{F}_{ijk} = \sum_{a=0}^{N_x+1} \sum_{b=0}^{N_y+1} \sum_{c=0}^{N_z+1} \beta_a^x \beta_b^y \beta_c^z f_{abc} g^x(k_i^x x_a) g^y(k_j^y y_b) g^z(k_k^z z_c), \qquad (17) $$

where the functions $g^x$, $g^y$, $g^z$, and the coefficients $k_i^x$, $k_j^y$, $k_k^z$, $\beta_a^x$, $\beta_b^y$ and $\beta_c^z$ depend on the type of boundary conditions. Let us denote by $\hat{F}$ the array of size $(N_x + 2) \times (N_y + 2) \times (N_z + 2)$ containing the discrete sine-cosine coefficients and by $f$ the array of size $(N_x + 2) \times (N_y + 2) \times (N_z + 2)$ containing the values of function $f$ at the grid points. The discrete transform (17) is then written as

$$ \hat{F} = D_{xyz}(f), \qquad (18) $$

where $x \in \{sI, cI, sIII, cIII\}$, $y \in \{sI, cI, sIII, cIII\}$ and $z \in \{sI, cI, sIII, cIII\}$, depending on the boundary conditions applied on the x-faces, y-faces and z-faces, respectively.

The associated inverse transform reads

$$ f_{abc} = \frac{8}{(N_x + 1)(N_y + 1)(N_z + 1)} \sum_{i=0}^{N_x+1} \sum_{j=0}^{N_y+1} \sum_{k=0}^{N_z+1} \gamma_i^x \gamma_j^y \gamma_k^z \hat{F}_{ijk} g^x(k_i^x x_i) g^y(k_j^y y_j) g^z(k_k^z z_k), \qquad (19) $$

and it can be written as

$$ f = D_{xyz}^{-1}(\hat{F}). \qquad (20) $$

A useful notation is finally introduced to account for dual mixed sine-cosine transforms. Assuming that the transform $D_{xyz}(\cdot)$ corresponds to Eq. (17), then we denote by

$$ \hat{F} = D_{\bar{xyz}}(f), \qquad (21) $$

the transform defined by

$$ \hat{F}_{ijk} = \sum_{a=0}^{N_x+1} \sum_{b=0}^{N_y+1} \sum_{c=0}^{N_z+1} \beta_a^x \beta_b^y \beta_c^z f_{abc} \overline{g^x(k_i^x x_i) g^y(k_j^y y_j) g^z(k_k^z z_k)}. \qquad (22) $$

The notation $D_{\bar{xyz}}$ defines the dual transform of $D_{xyz}$ in the x-direction. For example, if $D_{xyz}$ corresponds to the transform $D_{sIsIsI}$, then $D_{\bar{xyz}}$ corresponds to $D_{cIsIsI}$.

<footnote n="5">5 There is a total of $2^6 = 64$ different sine-cosine series in the 3-d case (corresponding to the 2 types of boundary conditions on the 6 faces).</footnote>
# 3. A fast numerical method for elasticity problems in heterogeneous media subjected to non-periodic boundary conditions

## 3.1. Equations of elasticity

We consider a problem of linear elasticity in a heterogeneous finite medium. The finite cell is a prismatic domain denoted by $\Omega = [0, L_x] \times [0, L_y] \times [0, L_z]$ in 3-d and a rectangle denoted by $\Omega = [0, L_x] \times [0, L_y]$ in 2-d. Tensorial components refer to a system of Cartesian coordinates $(e_x, e_y, e_z)$ in 3-d (and $(e_x, e_y)$ in 2-d).

The problem of elasticity consists in the computation of the elastic strain $\epsilon(\mathbf{u}(\mathbf{x}))$, stress $\sigma(\mathbf{x})$ and displacement field $\mathbf{u}(\mathbf{x})$, at each point $\mathbf{x}$ in $\Omega$, for a given elasticity tensor field $C(\mathbf{x})$. Since the principle of the method is to prescribe the value of the solution field or its normal derivative [29], we first consider mixed Dirichlet/Neumann boundary conditions on the components of the displacement field $\mathbf{u}(\mathbf{x})$, i.e.

$$ 
\left\{ 
\begin{array}{lll}
\forall \mathbf{x} \in \partial \Omega_{D_x}, & u_x(\mathbf{x}) = u_{\partial \Omega_{D_x}}(\mathbf{x}) & \text{(Dirichlet on } u_x) \\
\forall \mathbf{x} \in \partial \Omega_{N_x}, & \nabla u_x(\mathbf{x}) \cdot \mathbf{n} = \mathbf{f}_{\partial \Omega_{N_x}}(\mathbf{x}) & \text{(Neumann on } u_x) \\
\forall \mathbf{x} \in \partial \Omega_{D_y}, & u_y(\mathbf{x}) = u_{\partial \Omega_{D_y}}(\mathbf{x}) & \text{(Dirichlet on } u_y) \\
\forall \mathbf{x} \in \partial \Omega_{N_y}, & \nabla u_y(\mathbf{x}) \cdot \mathbf{n} = \mathbf{f}_{\partial \Omega_{N_y}}(\mathbf{x}) & \text{(Neumann on } u_y) \\
\forall \mathbf{x} \in \partial \Omega_{D_z}, & u_z(\mathbf{x}) = u_{\partial \Omega_{D_z}}(\mathbf{x}) & \text{(Dirichlet on } u_z) \\
\forall \mathbf{x} \in \partial \Omega_{N_z}, & \nabla u_z(\mathbf{x}) \cdot \mathbf{n} = \mathbf{f}_{\partial \Omega_{N_z}}(\mathbf{x}) & \text{(Neumann on } u_z),
\end{array}
\right. 
\qquad (23) $$

where $u_{\partial \Omega_{D_x}}$, $u_{\partial \Omega_{D_y}}$ and $u_{\partial \Omega_{D_z}}$ are the prescribed components of the displacement field (corresponding to the Dirichlet condition), and $\mathbf{f}_{\partial \Omega_{N_x}}$, $\mathbf{f}_{\partial \Omega_{N_y}}$ and $\mathbf{f}_{\partial \Omega_{N_z}}$ are the prescribed values for the normal derivative of components of the displacement field (corresponding to the Neumann condition on the displacement field). The boundaries associated to the field $u_x$ ($\partial \Omega_{D_x}$ and $\partial \Omega_{N_x}$) are necessary the union of entire faces with of course $\partial \Omega_{D_x} \cup \partial \Omega_{N_x} = \partial \Omega$ (the same properties also hold for the boundaries associated to the fields $u_y$ and $u_z$).

It must be noted that standard Neumann boundary conditions in elasticity involve non-normal derivatives of the displacement components, as they are applied through the normal stress vector, i.e.

$$ \sigma(\mathbf{x}) \cdot \mathbf{n} = C : \nabla \mathbf{u}(\mathbf{x}) \cdot \mathbf{n} = \mathbf{F}_{\partial \Omega_N}, \quad \forall \mathbf{x} \in \partial \Omega_N \qquad (24) $$

where $\mathbf{F}_{\partial \Omega_N}$ is the prescribed surface force and we used $C : \epsilon = C : \nabla \mathbf{u}$, due to the minor symmetry of the stiffness tensor $C$. Therefore, they cannot be treated by the present method as is. To overcome this issue, a thin external layer, with appropriate elastic properties that remove non-normal derivative of the displacement field in the expression of the stress-vector, is added to the domain. Expanding (24) on the boundary $x = L_x$ (without loss of generality) leads to

$$ 
\left\{ 
\begin{array}{ll}
C_{11kl}(L_x, y, z) \frac{\partial u_k}{\partial x_l}(L_x, y, z) = F_x(y, z) & \\
C_{12kl}(L_x, y, z) \frac{\partial u_k}{\partial x_l}(L_x, y, z) = F_y(y, z) & \\
C_{13kl}(L_x, y, z) \frac{\partial u_k}{\partial x_l}(L_x, y, z) = F_z(y, z). & 
\end{array}
\right. 
\qquad (25) $$

The non-normal derivatives are then artificially removed by considering a modified elasticity law in the external layer, leading to

$$ 
\left\{ 
\begin{array}{ll}
C_{11k1}(L_x, y, z) \frac{\partial u_k}{\partial x}(L_x, y, z) = F_x(y, z) & \\
C_{12k1}(L_x, y, z) \frac{\partial u_k}{\partial x}(L_x, y, z) = F_y(y, z) & \\
C_{13k1}(L_x, y, z) \frac{\partial u_k}{\partial x}(L_x, y, z) = F_z(y, z). & 
\end{array}
\right. 
\qquad (26) $$

This modification is made by considering a modified elasticity behavior in the extra-layer, denoted by $C^{\text{layer}}$ and defined by the following tensor

$$ C_{ijkl}^{\text{layer}} = \begin{cases} C_{ijkl} & \text{if } i = l \\ 0 & \text{if } i \neq l. \end{cases} \qquad (27) $$

One might remark that, contrary to the classical elasticity tensor, this modified tensor does not verify the minor and major symmetries. Still, it is a positive definite tensor, which ensures that the problem is well-posed under Neumann boundary conditions. For an infinitely thin extra layer, the modified problem is equivalent to the classical elasticity problem under Neumann boundary conditions. In practice, a one voxel layer will be used.

The local problem to be solved is

$$ 
\left\{ 
\begin{array}{ll}
\forall \mathbf{x} \in \Omega, & \operatorname{div} \sigma(\mathbf{x}) = 0 \\
\forall \mathbf{x} \in \Omega, & \sigma(\mathbf{x}) = C(\mathbf{x}) : \epsilon(\mathbf{u}(\mathbf{x})) \\
\forall \mathbf{x} \in \Omega, & \epsilon(\mathbf{u}(\mathbf{x})) = \frac{1}{2}(\nabla \mathbf{u} + \nabla^T \mathbf{u}),
\end{array}
\right. 
\qquad (28) $$

where the boundary conditions considered are given by (23). In Eq. (28), $C(\mathbf{x})$ is the fourth-order stiffness tensor. In the particular case of isotropic elasticity, it reads (in the case of a 3-d medium)

$$ C(\mathbf{x}) = 3\kappa(\mathbf{x})\mathbb{I} + 2\mu(\mathbf{x})\mathbb{K}, \qquad (29) $$
with $κ(x)$ and $μ(x)$ are respectively the local bulk and shear moduli. $J$ and $K$ are linearly independent isotropic tensors defined by

$$J = \frac{1}{3} I_3 \otimes I_3, \quad K = I - J, \quad (30)$$

with $I$ the fourth order identity tensor and $I_3$ is the second-order identity tensor.

The principle of resolution (see [29] for conductivity problems) is to

1. Write the local problem (28) as an auxiliary problem involving a reference homogeneous material and a polarization term (following the idea of Moulinec and Suquet [1] in the periodic setting);

2. Split each component of the displacement field $u$ in two contributions, a known term $u^0$ verifying the boundary conditions (23) (by construction) and an unknown fluctuation term $\tilde{u}$, defined such that it is null on the Dirichlet boundaries and it has normal derivatives null on the Neumann boundaries, to be determined;

3. Write the weak (Galerkin) formulation by (i) describing the fluctuation term by a partial 3-d mixed sine-cosine series (depending on the type of boundary conditions) and (ii) using trial functions based on the same 3-d mixed sine-cosine series;

4. Solve approximately the linear problem arising from the weak formulation (corresponding to the auxiliary problem for a given polarization tensor) using discrete sine and cosine transforms (to calculate, and in some cases approximate, the integrals defined in the weak formulation);

5. Find the polarization term solution of the problem using an iterative scheme (as done by Moulinec and Suquet [1] in the periodic setting with a fixed-point method).

## 3.2. The classical auxiliary problem and the associated difficulties

Following [1], the local elasticity problem (28), subjected to the boundary conditions (23), is written as the following auxiliary problem

$$\left\{
\begin{array}{l}
\forall x \in \Omega, \quad \operatorname{div} \sigma(x) = 0 \\
\forall x \in \Omega, \quad \sigma(x) = C_0 : \varepsilon(\mathbf{u}(x)) + \tau(x) \\
\forall x \in \Omega, \quad \varepsilon(\mathbf{u}(x)) = \frac{1}{2} (\nabla \mathbf{u}(x) + \nabla^T \mathbf{u}(x)),
\end{array}
\right.
\quad (31)$$

where $C^0$ is the stiffness tensor of a homogeneous reference material and

$$\tau(\mathbf{x}) = (C(\mathbf{x}) - C^0) : \varepsilon(\mathbf{u}(\mathbf{x})) \quad (32)$$

is the polarization tensor. At this stage, the homogeneous reference material is classically taken isotropic (with Lamé's coefficients $λ_0$ and $μ_0$) as in most works related to FFT (see e.g. [1]):

$$C^0_{ijkl} = \mu_0 (\delta_{ik}\delta_{jl} + \delta_{il}\delta_{kj}) + \lambda_0 \delta_{ij}\delta_{kl}, \quad \forall i, j, k, l = 1, 2, 3, \quad (33)$$

where $δ$ is the Kronecker symbol. Despite the fact that $C^0$ is isotropic one should remark that the (possible) anisotropy of the phases is still contained in the expression (32) of the polarization tensor.

Assuming (momentarily) that the polarization field $τ$ is known, the auxiliary problem defines a system of three equations on the unknown field $u$:

$$\left\{
\begin{array}{l}
(\lambda_0 + 2\mu_0) \frac{\partial^2 u_x}{\partial z^2} + \mu_0 \left( \frac{\partial^2 u_x}{\partial y^2} + \frac{\partial^2 u_x}{\partial z^2} \right) + (\lambda_0 + \mu_0) \left( \frac{\partial^2 u_y}{\partial x \partial y} + \frac{\partial^2 u_y}{\partial x \partial z} \right) + \frac{\partial \tau_{xx}}{\partial x} + \frac{\partial \tau_{xy}}{\partial y} + \frac{\partial \tau_{xz}}{\partial z} = 0 \\
(\lambda_0 + 2\mu_0) \frac{\partial^2 u_y}{\partial y^2} + \mu_0 \left( \frac{\partial^2 u_x}{\partial x^2} + \frac{\partial^2 u_x}{\partial z^2} \right) + (\lambda_0 + \mu_0) \left( \frac{\partial^2 u_x}{\partial y \partial x} + \frac{\partial^2 u_x}{\partial y \partial z} \right) + \frac{\partial \tau_{xy}}{\partial x} + \frac{\partial \tau_{yy}}{\partial y} + \frac{\partial \tau_{yz}}{\partial z} = 0 \\
(\lambda_0 + 2\mu_0) \frac{\partial^2 u_z}{\partial z^2} + \mu_0 \left( \frac{\partial^2 u_z}{\partial x^2} + \frac{\partial^2 u_z}{\partial y^2} \right) + (\lambda_0 + \mu_0) \left( \frac{\partial^2 u_x}{\partial z \partial x} + \frac{\partial^2 u_y}{\partial z \partial y} \right) + \frac{\partial \tau_{xz}}{\partial x} + \frac{\partial \tau_{yz}}{\partial y} + \frac{\partial \tau_{zz}}{\partial z} = 0.
\end{array}
\right.
\quad (34)$$

The boundary conditions are taken into account by considering the displacement field under the form

$$\mathbf{u}(\mathbf{x}) = \mathbf{u}_0(\mathbf{x}) + \tilde{\mathbf{u}}(\mathbf{x}), \quad \forall \mathbf{x} \in \Omega, \quad (35)$$

with $u_0$ a known field verifying the boundary conditions (see [29] for examples of reconstruction of the function, depending on the type of boundary conditions) and $\tilde{u}$ a fluctuation vector field which does not interfere with the boundary conditions. Each component of $\tilde{u}$ is thus described by an appropriate sine-cosine series as explained in Section 2:

$$\left\{
\begin{array}{l}
\tilde{u}_x(x, y, z) = \sum_{i=0}^{N_x+1} \sum_{j=0}^{N_y+1} \sum_{k=0}^{N_z+1} \alpha_i^{xx} \alpha_j^{xy} \alpha_k^{xz} U_{ijk}^x g^{xx}(k_i^{xx} x) g^{xy}(k_j^{xy} y) g^{xz}(k_k^{xz} z), \quad \forall (x, y, z) \in \Omega, \\
\\
\tilde{u}_y(x, y, z) = \sum_{i=0}^{N_x+1} \sum_{j=0}^{N_y+1} \sum_{k=0}^{N_z+1} \alpha_i^{yx} \alpha_j^{yy} \alpha_k^{yz} U_{ijk}^y g^{yx}(k_i^{yx} x) g^{yy}(k_j^{yy} y) g^{yz}(k_k^{yz} z), \quad \forall (x, y, z) \in \Omega, \\
\\
\tilde{u}_z(x, y, z) = \sum_{i=0}^{N_x+1} \sum_{j=0}^{N_y+1} \sum_{k=0}^{N_z+1} \alpha_i^{zx} \alpha_j^{zy} \alpha_k^{zz} U_{ijk}^z g^{zx}(k_i^{zx} x) g^{zy}(k_j^{zy} y) g^{zz}(k_k^{zz} z), \quad \forall (x, y, z) \in \Omega.
\end{array}
\right.
\quad (36)$$
In the description of the fluctuation (36), we consider a priori different possible boundary conditions on the components $\tilde{u}_x, \tilde{u}_y$ and $\tilde{u}_z$. Therefore, each component is separately described with its own combination of functions, which implies that nine different sine-cosine type functions are required and are denoted by $g^{xx}, g^{xy}, g^{xz}, g^{yx}, g^{yz}, g^{zx}, g^{zy}$ and $g^{zz}$; the first superscript corresponds to the component of $\tilde{u}$ that is modeled, and the second superscript corresponds to the spatial direction described by the function. This convention will also be adopted for the parameters $\alpha, k$ and $\xi$. The discrete sine-cosine transforms, respectively associated with this description of $\tilde{u}_x, \tilde{u}_y$ and $\tilde{u}_z$, are denoted by $D_{xyz}^x, D_{xyz}^y$ and $D_{xyz}^z$.

As explained in [29], the last coefficients $U_{ijk}^x, U_{ijk}^y$ and $U_{ijk}^z$ may be poorly estimated for $i = N_x + 1$ or $j = N_y + 1$ or $k = N_z + 1$ using discrete sine-cosine transforms. Therefore, the “high frequencies” are taken null:

$$U_{N_x+1jk}^x = U_{iN_y+1k}^x = U_{ijN_z+1}^x = U_{N_x+1jk}^y = U_{iN_y+1k}^y = U_{ijN_z+1}^y = U_{N_x+1jk}^z = U_{iN_y+1k}^z = U_{ijN_z+1}^z = 0. \quad (37)$$

The resolution thus consists in finding the “Fourier” modes associated with the three components of the displacement fluctuation, i.e., $U_{ijk}^x, U_{ijk}^y$ and $U_{ijk}^z$ (for $i = 0, ..., N_x, j = 0, ..., N_y$ and $k = 0, ..., N_z$). The principle of resolution is to assume that the polarization field $\tau$ is known and to solve, for this given value of the polarization, the equilibrium Eqs. (34) with $\mathbf{u}$ described by Eqs. (35) and (36). The resolution can then be done using a Galerkin-type method. The 3-d prismatic domain is discretized uniformly with $(N_x+2) \times (N_y+2) \times (N_z+2)$ points and the spatial scales associated with the uniform grid 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, ..., N_x+1$), $y_b = b\Delta y$ (for $b = 0, ..., N_y+1$) and $z_c = c\Delta z$ (for $c = 0, ..., N_z+1$).

The weak formulation of Eq. (31) reads

$$\int_{\Omega} (\mathrm{div} \, \sigma(\mathbf{x})) \cdot \mathbf{v}_{lmn}^i(\mathbf{x}) d\Omega = 0, \quad (38)$$

with $\mathbf{v}_{lmn}^i$, $i = 1, 2, 3$, $l = 0, ..., N_x$, $m = 0, ..., N_y$, $n = 0, ..., N_z$, the trial functions, reading

$$\mathbf{v}_{lmn}^1(\mathbf{x}) = \begin{pmatrix} v_{lmn}^x \\ 0 \\ 0 \end{pmatrix}; \quad \mathbf{v}_{lmn}^2(\mathbf{x}) = \begin{pmatrix} 0 \\ v_{lmn}^y \\ 0 \end{pmatrix}; \quad \mathbf{v}_{lmn}^3(\mathbf{x}) = \begin{pmatrix} 0 \\ 0 \\ v_{lmn}^z \end{pmatrix}. \quad (39)$$

with

$$v_{lmn}^x = g^{xx}(k_l^{xx})g^{xy}(k_m^{xy})g^{xz}(k_n^{xz})z, \quad (40)$$

$$v_{lmn}^y = g^{yx}(k_l^{yx})g^{yy}(k_m^{yy})g^{yz}(k_n^{yz})z, \quad (41)$$

$$v_{lmn}^z = g^{zx}(k_l^{zx})g^{zy}(k_m^{zy})g^{zz}(k_n^{zz})z. \quad (42)$$

This leads to a system of three equations which can be written under the generic form

$$\left\{ \begin{aligned} I_{lmn}^x + J_{lmn}^x &= -K_{lmn}^x - I_{lmn}^{x0} - J_{lmn}^{x0} \\ I_{lmn}^y + J_{lmn}^y &= -K_{lmn}^y - I_{lmn}^{y0} - J_{lmn}^{y0} \\ I_{lmn}^z + J_{lmn}^z &= -K_{lmn}^z - I_{lmn}^{z0} - J_{lmn}^{z0} \end{aligned} \right. \quad (43)$$

with

$$\left\{ \begin{aligned} I_{lmn}^x &= \int_{\Omega} \left( (\lambda_0 + 2\mu_0) \frac{\partial^2 \tilde{u}_x}{\partial x^2}(\mathbf{x}) + \mu_0 \left( \frac{\partial^2 \tilde{u}_x}{\partial y^2}(\mathbf{x}) + \frac{\partial^2 \tilde{u}_x}{\partial z^2}(\mathbf{x}) \right) \right) v_{lmn}^x(\mathbf{x}) \, d\Omega, \\ I_{lmn}^{x0} &= \int_{\Omega} \left( (\lambda_0 + 2\mu_0) \frac{\partial^2 u_x^0}{\partial x^2}(\mathbf{x}) + \mu_0 \left( \frac{\partial^2 u_x^0}{\partial y^2}(\mathbf{x}) + \frac{\partial^2 u_x^0}{\partial z^2}(\mathbf{x}) \right) \right) v_{lmn}^x(\mathbf{x}) \, d\Omega, \\ J_{lmn}^x &= (\lambda_0 + \mu_0) \int_{\Omega} \left( \frac{\partial^2 \tilde{u}_y}{\partial x \partial y}(\mathbf{x}) + \frac{\partial^2 \tilde{u}_z}{\partial x \partial z}(\mathbf{x}) \right) v_{lmn}^x(\mathbf{x}) \, d\Omega, \\ J_{lmn}^{x0} &= (\lambda_0 + \mu_0) \int_{\Omega} \left( \frac{\partial^2 u_y^0}{\partial x \partial y}(\mathbf{x}) + \frac{\partial^2 u_z^0}{\partial x \partial z}(\mathbf{x}) \right) v_{lmn}^x(\mathbf{x}) \, d\Omega, \\ K_{lmn}^x &= \int_{\Omega} \left( \frac{\partial \tau_{xx}}{\partial x}(\mathbf{x}) + \frac{\partial \tau_{xy}}{\partial y}(\mathbf{x}) + \frac{\partial \tau_{zx}}{\partial z}(\mathbf{x}) \right) v_{lmn}^x(\mathbf{x}) \, d\Omega. \end{aligned} \right. \quad (44)$$

Similar expressions for $I_{lmn}^y, J_{lmn}^y, K_{lmn}^y, I_{lmn}^{y0}, J_{lmn}^{y0}, I_{lmn}^z, J_{lmn}^z, K_{lmn}^z, I_{lmn}^{z0}$ and $J_{lmn}^{z0}$ can be written using Eq. (34). One can notice that only $I_{lmn}^x, J_{lmn}^x, I_{lmn}^y, J_{lmn}^y, I_{lmn}^z, J_{lmn}^z$ are unknown terms. The calculation of these integrals will be done analytically or numerically, depending on the integral:

* The integrals of type $I_{lmn}^x$ and $J_{lmn}^x$ are calculated analytically as they are based on the calculation of elementary integrals involving different type of sine and cosine functions (related to the orthogonality of the modes) which are provided in Appendix A.

* The integrals of type $I_{lmn}^{x0}, J_{lmn}^{x0}$ and $K_{lmn}^x$ involve functions that are defined (a priori) numerically ($\tau$ and $\mathbf{u}_0$); therefore these integrals will be calculated approximately using discrete sine-cosine transforms (defined in Section 2); this is similar to the standard numerical integration of the Moulinec and Suquet [1] scheme in the periodic setting using discrete Fourier transforms (see [35]).
Hereafter, we discuss the integrals $J_{lmn}^x$ and $J_{lmn}^y$, as this will be very similar for the other quantities. The detailed calculation of $J_{lmn}^x$ is given in Appendix B and leads to

$$I_{lmn}^x = -\frac{L_x L_y L_z}{8} \left( (\lambda_0 + 2\mu_0)(\xi_l^{xy})^2 + \mu_0 \left( (\xi_m^{xy})^2 + (\xi_n^{xz})^2 \right) \right) U_{lmn}^x \quad (45)$$

The computation of $J_{lmn}^x$ is more troublesome as it implies the functions $\tilde{u}_y, \tilde{u}_z$ and $v_{lmn}^x$, which, depending on the boundary conditions, can be of different types of sine–cosine series; this leads to two cases:

(1) If $\partial^2 u_y / \partial x \partial y$, $\partial^2 u_z / \partial x \partial z$ and $v_{lmn}^x$ are of the same type of sine–cosine series, $J_{lmn}^x$ is given by

$$J_{lmn}^x = \frac{L_x L_y L_z}{8} (\lambda_0 + \mu_0) (\xi_l^{yx} \xi_m^{yy} U_{lmn}^y + \xi_l^{zx} \xi_n^{zz} U_{lmn}^z), \quad (46)$$

and similar relations are obtained for $J_{lmn}^y$ and $J_{lmn}^z$. Thus, the three equations defining (43) for a triplet $(l, m, n)$ only involve the three coefficients $U_{lmn}^x, U_{lmn}^y$ and $U_{lmn}^z$, leading to an easy-to-solve system of three equations with three unknowns, which has to be solved for every values of $l, m$ and $n$. One might remark that these situations correspond to the periodicity compatible mixed uniform boundary conditions (PMUBC) treated in [23] using mirror unit cells.

(2) Otherwise, the computation of $J_{lmn}^x$ relies on “cross-mode” integrals (see Appendix A), implying that $J_{lmn}^x$ will depend on terms $U_{opq}^x, U_{opq}^y$ and $U_{opq}^z$ with $(o, p, q) \neq (l, m, n)$. As a consequence, the weak formulation (38) will lead to a hard-to-solve (non diagonal, non diagonal by blocks, non sparse) linear system of equations, which will require high memory and computational time to solve. Thus, in the general case, the resolution of the standard auxiliary problem (31) leads to a computational cost higher than finite element method (and the main advantage of FFT-based methods is lost).

The reason why this auxiliary problem is not suitable is the presence of mixed derivatives in Eq. (34) which leads to a non-diagonal linear system. A modified auxiliary problem is therefore required to keep the computational complexity of FFT.

### 3.3. The modified auxiliary problem

To resolve the issue explained in Section 3.2, the mixed derivatives involved in Eq. (34) are “transferred” in the polarization tensor. To do so, the elasticity problem is first expressed as a function of the gradient of the displacement field (unsymmetrized). Eq. (28) is rewritten as

$$\left\{ \begin{array}{ll} \mathrm{div} \, \sigma(\mathbf{x}) = 0 & \forall \mathbf{x} \in \Omega \\ \sigma(\mathbf{x}) = C(\mathbf{x}) : \nabla \mathbf{u} & \forall \mathbf{x} \in \Omega, \end{array} \right. \qquad (47)$$

since $C : \varepsilon = C : \nabla \mathbf{u}$ due to the minor symmetry of the stiffness tensor $C$. Then, we introduce a modified polarization tensor as

$$\mathbf{\tau}(\mathbf{x}) = (C(\mathbf{x}) - \mathbb{B}^0) : \nabla \mathbf{u}(\mathbf{x}), \quad \forall \mathbf{x} \in \Omega, \qquad (48)$$

where $\mathbb{B}^0$ is a fourth-order tensor reading

$$B_{ijkl}^{0} = \mu_{0}\delta_{ik}\delta_{jl} + (\lambda_{0} + \mu_{0})\delta_{ik}\delta_{jl}\delta_{kl}, \quad \forall i, j, k, l = 1, 2, 3, \quad (\text{no summation on the indices}), \qquad (49)$$

which has been obtained by putting to zero the components related to cross-derivatives of $C^0$ (given by Eq. (33)). One might remark that, contrary to $C^0$, $\mathbb{B}^0$ does not verify the minor and major symmetries; it is thus only a mathematical tool serving the numerical resolution of the auxiliary problem.

The modified auxiliary problem thus reads

$$\left\{ \begin{array}{ll} \mathrm{div} \, \sigma(\mathbf{x}) = \mathbf{0} & \forall \mathbf{x} \in \Omega \\ \sigma(\mathbf{x}) = \mathbb{B}^0 : \nabla \mathbf{u}(\mathbf{x}) + \mathbf{\tau}(\mathbf{x}) & \forall \mathbf{x} \in \Omega. \end{array} \right. \qquad (50)$$

One must remark that $\mathbb{B}^0$ is a positive definite tensor, since

$$\nabla \mathbf{u} : \mathbb{B}^0 : \nabla \mathbf{u} = (\lambda_0 + 2\mu_0) \sum_{i=1}^{3} \left( \frac{\partial u_i}{\partial x_i} \right)^2 + \mu_0 \sum_{i \neq j} \left( \frac{\partial u_i}{\partial x_j} \right)^2, \qquad (51)$$

which ensures the existence and the uniqueness of a solution for the modified auxiliary problem. As previously, assuming (momentarily) that the polarization field $\mathbf{\tau}$ is known, the modified auxiliary problem defines a system of three equations on the unknown field $\mathbf{u}$, which is freed from mixed derivatives:

$$\left\{ \begin{array}{l} (\lambda_0 + 2\mu_0) \frac{\partial^2 u_x}{\partial x^2} + \mu_0 \left( \frac{\partial^2 u_x}{\partial y^2} + \frac{\partial^2 u_x}{\partial z^2} \right) + \frac{\partial \tau_{xx}}{\partial x} + \frac{\partial \tau_{xy}}{\partial y} + \frac{\partial \tau_{xz}}{\partial z} = 0 \\ (\lambda_0 + 2\mu_0) \frac{\partial^2 u_y}{\partial y^2} + \mu_0 \left( \frac{\partial^2 u_y}{\partial x^2} + \frac{\partial^2 u_y}{\partial z^2} \right) + \frac{\partial \tau_{xy}}{\partial x} + \frac{\partial \tau_{yy}}{\partial y} + \frac{\partial \tau_{yz}}{\partial z} = 0 \\ (\lambda_0 + 2\mu_0) \frac{\partial^2 u_z}{\partial z^2} + \mu_0 \left( \frac{\partial^2 u_z}{\partial x^2} + \frac{\partial^2 u_z}{\partial y^2} \right) + \frac{\partial \tau_{xz}}{\partial x} + \frac{\partial \tau_{yz}}{\partial y} + \frac{\partial \tau_{zz}}{\partial z} = 0. \end{array} \right. \qquad (52)$$

The weak formulation associated with the system of partial differential equations reduces to

$$\left\{ \begin{array}{l} I_{lmn}^x = -K_{lmn}^x - I_{lmn}^{x0} \\ I_{lmn}^y = -K_{lmn}^y - I_{lmn}^{y0} \\ I_{lmn}^z = -K_{lmn}^z - I_{lmn}^{z0} \end{array} \right. \qquad (53)$$
where $I_{lmn}^x$, $K_{lmn}^x$ and $I_{lmn}^{x0}$ are still given by Eq. (44). Considering Eq. (45) (and similar relations for $I_{lmn}^y$ and $I_{lmn}^z$), the resolution of the weak form of the auxiliary problem read

$$ egin{cases} U_{lmn}^x = \frac{8}{L_x L_y L_z} \left( (\lambda_0 + 2\mu_0) (\xi_l^{xy})^2 + \mu_0 ((\xi_m^{xy})^2 + (\xi_n^{xz})^2) \right)^{-1} (K_{lmn}^x + I_{lmn}^{x0}) \\ U_{lmn}^y = \frac{8}{L_x L_y L_z} \left( (\lambda_0 + 2\mu_0) (\xi_m^{yy})^2 + \mu_0 ((\xi_l^{yx})^2 + (\xi_n^{yz})^2) \right)^{-1} (K_{lmn}^y + I_{lmn}^{y0}) \\ U_{lmn}^z = \frac{8}{L_x L_y L_z} \left( (\lambda_0 + 2\mu_0) (\xi_n^{zz})^2 + \mu_0 ((\xi_l^{zx})^2 + (\xi_m^{zy})^2) \right)^{-1} (K_{lmn}^z + I_{lmn}^{z0}) \end{cases}, \quad \forall l, m, n. \quad (54) $$

The detailed calculations of $K_{lmn}^x$ and $I_{lmn}^{x0}$ are given in Appendix B and rely on the use of discrete sine-cosine transforms (similar relations can be obtained for $K_{lmn}^y$, $I_{lmn}^{y0}$, $K_{lmn}^z$ and $I_{lmn}^{z0}$).

Let us denote by $U^x$, $U^y$ and $U^z$ the arrays of size $(N_x + 2) \times (N_y + 2) \times (N_z + 2)$ containing the values of the sine-cosine series coefficients $U_{lmn}^x$, $U_{lmn}^y$ and $U_{lmn}^z$ respectively. Eq. (54) can be written alternatively

$$ egin{cases} U^x = A^x \odot \left[ D_{xyz}^x (\Delta \mathbf{u}_x^0) + \xi^{xx} \odot D_{\bar{xyz}}^x (\tilde{\tau}_{xx}) + \xi^{xy} \odot D_{x\bar{yz}}^x (\tilde{\tau}_{xy}) + \xi^{xz} \odot D_{xy\bar{z}}^x (\tilde{\tau}_{xz}) + S^x \right] \\ U^y = A^y \odot \left[ D_{xyz}^y (\Delta \mathbf{u}_y^0) + \xi^{yx} \odot D_{\bar{xyz}}^y (\tilde{\tau}_{yx}) + \xi^{yy} \odot D_{x\bar{yz}}^y (\tilde{\tau}_{yy}) + \xi^{yz} \odot D_{xyz}^y (\tilde{\tau}_{yz}) + S^y \right] \\ U^z = A^z \odot \left[ D_{xyz}^z (\Delta \mathbf{u}_z^0) + \xi^{zx} \odot D_{\bar{xyz}}^z (\tilde{\tau}_{zx}) + \xi^{zy} \odot D_{x\bar{yz}}^z (\tilde{\tau}_{yy}) + \xi^{zz} \odot D_{xyz}^z (\tilde{\tau}_{zz}) + S^z \right] \end{cases} \quad (55) $$

where $A^x$, $A^y$ and $A^z$ are the $(N_x + 2) \times (N_y + 2) \times (N_z + 2)$ arrays given by

$$ egin{cases} (A)_{ijk}^x = \frac{1}{(\lambda_0 + 2\mu_0) (\xi_l^{xx})^2 + \mu_0 ((\xi_j^{xy})^2 + (\xi_k^{xz})^2)} \times \frac{8}{(N_x + 1)(N_y + 1)(N_z + 1)} \\ (A)_{ijk}^y = \frac{1}{(\lambda_0 + 2\mu_0) (\xi_j^{yy})^2 + \mu_0 ((\xi_l^{yx})^2 + (\xi_k^{yz})^2)} \times \frac{8}{(N_x + 1)(N_y + 1)(N_z + 1)} \\ (A)_{ijk}^z = \frac{1}{(\lambda_0 + 2\mu_0) (\xi_k^{zz})^2 + \mu_0 ((\xi_l^{zx})^2 + (\xi_j^{zy})^2)} \times \frac{8}{(N_x + 1)(N_y + 1)(N_z + 1)}, \end{cases} \quad (56) $$

and $\xi^{xx}$, $\xi^{xy}$, $\xi^{xz}$, $\xi^{yx}$, $\xi^{yy}$, $\xi^{yz}$, $\xi^{zx}$, $\xi^{zy}$, $\xi^{zz}$ are the $(N_x + 2) \times (N_y + 2) \times (N_z + 2)$ arrays given by

$$ egin{cases} (\xi^{xx})_{ijk} = \xi_l^{xx}, \quad (\xi^{xy})_{ijk} = \xi_j^{xy}, \quad (\xi^{xz})_{ijk} = \xi_k^{xz}, \\ (\xi^{yx})_{ijk} = \xi_l^{yx}, \quad (\xi^{yy})_{ijk} = \xi_j^{yy}, \quad (\xi^{yz})_{ijk} = \xi_k^{yz}, \\ (\xi^{zx})_{ijk} = \xi_l^{zx}, \quad (\xi^{zy})_{ijk} = \xi_j^{zy}, \quad (\xi^{zz})_{ijk} = \xi_k^{zz}. \end{cases} \quad (57) $$

The $(N_x + 2) \times (N_y + 2) \times (N_z + 2)$ arrays $S^x$, $S^y$ and $S^z$ contain the values $S_{lmn}^x$, $S_{lmn}^y$ and $S_{lmn}^z$ (see Appendix B) and $\odot$ denotes Hadamard product (pointwise product) given by

$$ (A \odot B)_{ijk} = (A)_{ijk}(B)_{ijk} \quad (\text{no summation on the indices}). \quad (58) $$

It should be noted that Eq. (55) allows the calculation of $U_{lmn}^x$, $U_{lmn}^y$ and $U_{lmn}^z$ for all frequencies except for $l = N_x + 1$, $m = N_y + 1$ and $n = N_z + 1$ where Eq. (37) holds. Once the arrays $U^x$, $U^y$ and $U^z$ are known, each component of the displacement fluctuation $\tilde{u}$ is then known using Eq. (36).

### 3.4. Iterative scheme

As in the periodic case [1], the determination of the polarization tensor $\tau$ solution of the heterogeneous problem is based on an iterative scheme. We consider a fixed-point iterative scheme based on the modified auxiliary problem. One might remark that, due
to the modification of the auxiliary problem, the convergence of the iterative scheme is not guaranteed a priori and must be proven
by a theoretical proof, and the convergence rate should be assessed by numerical tests. The iterative scheme reads:

$$
\left\{
\begin{array}{ll}
\text{Initialization} & \tilde{\mathbf{u}}^0(\mathbf{x}) = \mathbf{0} \\
& \sigma^0(\mathbf{x}) = \mathbb{C}(\mathbf{x}) : \nabla \mathbf{u}_0(\mathbf{x}) \\
\text{Iterate } n+1 & \tilde{\mathbf{u}}^n \text{ and } \sigma^n \text{ being known} \\
\text{(a)} & (\nabla \mathbf{u}_0(\mathbf{x}) + \nabla \tilde{\mathbf{u}}^n(\mathbf{x}))^n = \sigma^n(\mathbf{x}) - \mathbb{B}^0 : (\nabla \mathbf{u}_0(\mathbf{x}) + \nabla \tilde{\mathbf{u}}^n(\mathbf{x})) \\
\text{(b)} & (U^x)^{n+1} = \mathbf{A}^x \odot \left[ D_{xyz}^x (\Delta \mathbf{u}_x^0) + \xi^{xx} \odot D_{\bar{xy}z}^x (\tilde{\mathbf{r}}_{xx}^n) + \xi^{xy} \odot D_{\bar{xy}z}^x (\tilde{\mathbf{r}}_{xy}^n) + \xi^{xz} \odot D_{xyz}^x (\tilde{\mathbf{r}}_{xz}^n) + \mathbf{S}^x \right] \\
& (U^y)^{n+1} = \mathbf{A}^y \odot \left[ D_{xyz}^y (\Delta \mathbf{u}_y^0) + \xi^{yx} \odot D_{\bar{xy}z}^y (\tilde{\mathbf{r}}_{yx}^n) + \xi^{yy} \odot D_{\bar{xy}z}^y (\tilde{\mathbf{r}}_{yy}^n) + \xi^{yz} \odot D_{xyz}^y (\tilde{\mathbf{r}}_{yz}^n) + \mathbf{S}^y \right] \\
& (U^z)^{n+1} = \mathbf{A}^z \odot \left[ D_{xyz}^z (\Delta \mathbf{u}_z^0) + \xi^{zx} \odot D_{\bar{xy}z}^z (\tilde{\mathbf{r}}_{zx}^n) + \xi^{zy} \odot D_{\bar{xy}z}^z (\tilde{\mathbf{r}}_{zy}^n) + \xi^{zz} \odot D_{xyz}^z (\tilde{\mathbf{r}}_{zz}^n) + \mathbf{S}^z \right] \\
\text{(c)} & \tilde{\mathbf{u}}_x^{n+1}(\mathbf{x}, \mathbf{y}, \mathbf{z}) = \sum_{l=0}^{N_x+1} \sum_{j=0}^{N_y+1} \sum_{k=0}^{N_z+1} \alpha_l^{xx} \alpha_j^{xy} \alpha_k^{xz} (U_{ijk}^x)^{n+1} g^{xx}(k_l^{xx} \mathbf{x}) g^{xy}(k_j^{xy} \mathbf{y}) g^{xz}(k_k^{xz} \mathbf{z}) \\
& \tilde{\mathbf{u}}_y^{n+1}(\mathbf{x}, \mathbf{y}, \mathbf{z}) = \sum_{l=0}^{N_x+1} \sum_{j=0}^{N_y+1} \sum_{k=0}^{N_z+1} \alpha_l^{yx} \alpha_j^{yy} \alpha_k^{yz} (U_{ijk}^y)^{n+1} g^{yx}(k_l^{yx} \mathbf{x}) g^{yy}(k_j^{yy} \mathbf{y}) g^{yz}(k_k^{yz} \mathbf{z}) \\
& \tilde{\mathbf{u}}_z^{n+1}(\mathbf{x}, \mathbf{y}, \mathbf{z}) = \sum_{l=0}^{N_x+1} \sum_{j=0}^{N_y+1} \sum_{k=0}^{N_z+1} \alpha_l^{zx} \alpha_j^{zy} \alpha_k^{zz} (U_{ijk}^z)^{n+1} g^{zx}(k_l^{zx} \mathbf{x}) g^{zy}(k_j^{zy} \mathbf{y}) g^{zz}(k_k^{zz} \mathbf{z}) \\
\text{(d)} & \sigma^{n+1}(\mathbf{x}) = \mathbb{C}(\mathbf{x}) : (\nabla \mathbf{u}_0(\mathbf{x}) + \nabla \tilde{\mathbf{u}}^{n+1}(\mathbf{x})) \\
\text{(e)} & \text{Convergence test.}
\end{array}
\right.
\quad (59)
$$

The convergence test consists in verifying the local equilibrium. Then, following Parseval's identity, the $L^2$ norm of $\operatorname{div} \sigma(\mathbf{x})$
reads

$$
\| \operatorname{div} \sigma \|_{L^2}^2 = \frac{1}{L_x L_y L_z} \int_{\Omega} (\operatorname{div} \sigma(\mathbf{x}))^2 d\Omega = \sum_{l=0}^{N_y} \sum_{j=0}^{N_y} \sum_{k=0}^{N_z} \left( (L_{ijk}^x)^2 + (L_{ijk}^y)^2 + (L_{ijk}^z)^2 \right) \quad (60)
$$

where $L_{ijk}^x$, $L_{ijk}^y$ and $L_{ijk}^z$ are given by

$$
\left\{
\begin{aligned}
L_{ijk}^x &= \frac{8}{L_x L_y L_z} \int_{\Omega} \left( \frac{\partial \sigma_{xx}}{\partial x} + \frac{\partial \sigma_{xy}}{\partial y} + \frac{\partial \sigma_{xz}}{\partial z} \right) v_{ijk}^x(\mathbf{x}) d\Omega \\
L_{ijk}^y &= \frac{8}{L_x L_y L_z} \int_{\Omega} \left( \frac{\partial \sigma_{xy}}{\partial x} + \frac{\partial \sigma_{yy}}{\partial y} + \frac{\partial \sigma_{yz}}{\partial z} \right) v_{ijk}^y(\mathbf{x}) d\Omega \\
L_{ijk}^z &= \frac{8}{L_x L_y L_z} \int_{\Omega} \left( \frac{\partial \sigma_{xz}}{\partial x} + \frac{\partial \sigma_{yz}}{\partial y} + \frac{\partial \sigma_{zz}}{\partial z} \right) v_{ijk}^z(\mathbf{x}) d\Omega.
\end{aligned}
\right.
\qquad (61)
$$

The integrals (61) are computed as described in Appendix B for div **τ**, and the convergence is considered achieved when [27]

$$
\frac{\Delta x \sqrt{\sum_{l=0}^{N_x} \sum_{j=0}^{N_y} \sum_{k=0}^{N_z} \left( (L_{ijk}^x)^2 + (L_{ijk}^y)^2 + (L_{ijk}^z)^2 \right)}}{\sqrt{\sum_{l=0}^{N_x} \sum_{j=0}^{N_y} \sum_{k=0}^{N_z} (\sigma_{xx,ijl}^2 + \sigma_{yy,ijl}^2 + \sigma_{zz,ijl}^2 + 2\sigma_{xy,ijl}^2 + 2\sigma_{xz,ijl}^2 + 2\sigma_{yz,ijl}^2)}} \le \text{tol}, \quad (62)
$$

where tol is the tolerance taken at 10⁻⁸ in the following.

3.5. *Proof of convergence associated with the modified auxiliary problem*

As the proposed modification of the polarization tensor may alter the convergence of the iterative scheme, a proof of convergence
is required to (i) ensure the convergence of the method and (ii) provide guidelines for the determination of an appropriate
comparison material. Following the convergence proof for the classical FFT-based method proposed by Michel et al. [39], we look
into the eigenvalues of the operator leading from an iteration to another. In the periodic framework proposed by Moulinec and
Suquet [1], which operates on the strain field, it is referred as the Green operator and noted I⁰. It is slightly different in this work
as we operate on the displacement field and with different boundary conditions.

For given boundary conditions (either Dirichlet and/or Neumann), let us consider the linear operator Λ⁰ giving the solution **u**
of the auxiliary problem (50) under null boundary conditions and for a given polarization tensor **τ**. Then,

$$
\mathbf{u} = \mathbf{\Lambda}^0(\mathbf{\tau}) \tag{63}
$$
is the unique solution of (50) under null boundary conditions and therefore verifies

$$ \mathrm{div}(\mathbb{B}^0 : \nabla \mathbf{u}(\mathbf{x})) = -\mathrm{div}(\tau(\mathbf{x})), \quad \forall \mathbf{x}. \qquad (64) $$

For further use, we remark that, for a polarization field given by

$$ \tau(\mathbf{x}) = -\mathbb{B}^0 : \nabla \mathbf{u}(\mathbf{x}), \quad \forall \mathbf{x}. \qquad (65) $$

with $\mathbf{u}$ a given displacement field with null boundary conditions, $\mathbf{u}$ is the solution of (64), leading to the identity

$$ \mathbf{u} = \Lambda^0(-\mathbb{B}^0 : \nabla \mathbf{u}). \qquad (66) $$

Considering the algorithm presented in Section 3.4, the fluctuation field at iteration $n+1$ is solution of

$$ \mathrm{div}(\mathbb{B}^0 : \nabla \tilde{\mathbf{u}}^{n+1}(\mathbf{x})) = -\mathrm{div}(\tau''(\mathbf{x}) + \mathbb{B}^0 : \nabla \mathbf{u}^0) \quad \forall \mathbf{x}, \qquad (67) $$

under null boundary conditions. Then, the iterative scheme can be expressed as

$$ \tilde{\mathbf{u}}^{n+1} = \Lambda^0(\tau^n + \mathbb{B}^0 : \nabla \mathbf{u}^0) = \Lambda^0((\mathbb{C} - \mathbb{B}^0) : \nabla(\mathbf{u}^0 + \tilde{\mathbf{u}}^n) + \mathbb{B}^0 : \nabla \mathbf{u}^0). \qquad (68) $$

Let us split the fluctuation field $\tilde{\mathbf{u}}^n$ at iteration $n$ into the fluctuation of the solution field of the actual problem, denoted as $\tilde{\mathbf{u}}_{\mathrm{sol}}$, and the error field at iteration $n$, denoted as $\tilde{\mathbf{u}}_{\mathrm{err}}$. It reads

$$ \tilde{\mathbf{u}}^n = \tilde{\mathbf{u}}_{\mathrm{sol}} + \tilde{\mathbf{u}}_{\mathrm{err}}^n. \qquad (69) $$

Following (68), the solution of the actual problem must verify

$$ \tilde{\mathbf{u}}_{\mathrm{sol}} = \Lambda^0 ((\mathbb{C} - \mathbb{B}^0) : \nabla(\mathbf{u}^0 + \tilde{\mathbf{u}}_{\mathrm{sol}}) + \mathbb{B}^0 : \nabla \mathbf{u}^0). \qquad (70) $$

Considering (68), (69) and (70), the evolution of the error field with the iterations can be expressed as

$$ \tilde{\mathbf{u}}_{\mathrm{err}}^{n+1} = \Lambda^0 ((\mathbb{C} - \mathbb{B}^0) : \nabla \tilde{\mathbf{u}}_{\mathrm{err}}^n). \qquad (71) $$

Thus, one can ensure the convergence of the iterative scheme by proving that the operator

$$ \tilde{\mathbf{u}}_{\mathrm{err}} \rightarrow \Lambda^0 ((\mathbb{C} - \mathbb{B}^0) : \nabla \tilde{\mathbf{u}}_{\mathrm{err}}), \qquad (72) $$

is a contracting operator, i.e. that its eigenvalues lie between minus one and one. Let us consider an eigenvector $\tilde{\mathbf{u}}^a$ of the operator associated with an eigenvalue $\alpha$, i.e. such that

$$ \alpha \tilde{\mathbf{u}}^a = \Lambda^0 ((\mathbb{C} - \mathbb{B}^0) : \nabla \tilde{\mathbf{u}}^a). \qquad (73) $$

Following the reasoning presented in [39], and taking advantage of the identity (66) to obtain $\Lambda^0((1 - \alpha)\mathbb{B}^0 : \nabla \tilde{\mathbf{u}}^a) = -(1 - \alpha)\tilde{\mathbf{u}}^a$, one can rewrite (73) into

$$ \tilde{\mathbf{u}}^a = \Lambda^0 (((\mathbb{C} - (1 - \alpha)\mathbb{B}^0) - \mathbb{B}^0) : \nabla \tilde{\mathbf{u}}^a), \qquad (74) $$

which means that $\tilde{\mathbf{u}}^a$ is solution of

$$ \begin{cases} \mathrm{div} \, \sigma(\mathbf{x}) = 0 & \forall \mathbf{x} \in \Omega \\ \sigma(\mathbf{x}) = (\mathbb{C}(\mathbf{x}) - (1 - \alpha)\mathbb{B}^0) : \nabla \mathbf{u}(\mathbf{x}) & \forall \mathbf{x} \in \Omega, \end{cases} \qquad (75) $$

under null boundary conditions. If $\mathbb{C}(\mathbf{x}) - (1 - \alpha)\mathbb{B}^0$ is positive definite (resp. negative definite) $\forall \mathbf{x} \in \Omega$, Eq. (75) has a unique solution, which is $\tilde{\mathbf{u}}^a = 0$. Thus, for every eigenvalues $\alpha$, there are $\mathbf{x}$ such that $\mathbb{C}(\mathbf{x}) - (1 - \alpha)\mathbb{B}^0$ has negative eigenvalues and $\mathbf{x}$ such that $\mathbb{C}(\mathbf{x}) - (1 - \alpha)\mathbb{B}^0$ has positive eigenvalues. To study the eigenvalues of $\mathbb{C}(\mathbf{x}) - (1 - \alpha)\mathbb{B}^0$, we consider, as done by Michel et al. [39], the particular case of $\mathbb{C}$ isotropic (i.e. that $C_{ijkl}(\mathbf{x}) = \lambda(\mathbf{x})\delta_{ij}\delta_{kl} + \mu(\mathbf{x})(\delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk})$). Eq. (75)$_2$ is then written under a matrix form, reading

$$ \left( \begin{array}{c} \sigma_{xx} \\ \sigma_{yy} \\ \sigma_{zz} \\ \sigma_{xy} \\ \sigma_{yx} \\ \sigma_{zx} \\ \sigma_{zx} \\ \sigma_{yz} \\ \sigma_{zy} \end{array} \right) = \left( \begin{array}{ccccccccc} \Delta_{\mu\lambda} & \lambda & \lambda & 0 & 0 & 0 & 0 & 0 & 0 \\ \lambda & \Delta_{\mu\lambda} & \lambda & 0 & 0 & 0 & 0 & 0 & 0 \\ \lambda & \lambda & \Delta_{\mu\lambda} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \Delta_{\mu} & \mu & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \mu & \Delta_{\mu} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & \Delta_{\mu} & \mu & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & \Delta_{\mu} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & \Delta_{\mu} & \mu \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & \mu & \Delta_{\mu} \end{array} \right) \left( \begin{array}{c} \partial u_x / \partial x \\ \partial u_y / \partial y \\ \partial u_z / \partial z \\ \partial u_y / \partial x \\ \partial u_x / \partial y \\ \partial u_z / \partial x \\ \partial u_x / \partial z \\ \partial u_z / \partial y \\ \partial u_x / \partial z \end{array} \right) \qquad (76) $$

with

$$ \Delta_{\mu\lambda} = 2\mu(\mathbf{x}) + \lambda(\mathbf{x}) - (1 - \alpha)(2\mu^0 + \lambda^0) ; \quad \Delta_{\mu} = \mu(\mathbf{x}) - (1 - \alpha)\mu^0. \qquad (77) $$
The eigenvalues of the matrix in Eq. (76) read

$$ 
\left\{ 
\begin{aligned}
e_1 &= 2\mu(\mathbf{x}) - (1-\alpha)(2\mu_0 + \lambda_0) \\
e_2 &= 2\mu(\mathbf{x}) + 3\lambda(\mathbf{x}) - (1-\alpha)(2\mu_0 + \lambda_0) \\
e_3 &= 2\mu(\mathbf{x}) - (1-\alpha)\mu_0 \\
e_4 &= -(1-\alpha)\mu_0.
\end{aligned}
\right. 
\qquad (78) 
$$

First, it is easy to remark that $\alpha \ge 1$ leads to $e_i \ge 0$ for $i = 1, 2, 3, 4$, which is excluded. In the case $\alpha < 1$, one has $e_4 < 0$. Therefore, the existence of an eigenvalue $\alpha$ requires $e_1 \ge 0$ or $e_2 \ge 0$ or $e_3 \ge 0$, leading to:

$$ \alpha \ge 1 - \frac{2\mu(\mathbf{x})}{2\mu^0 + \lambda^0} \quad \text{or} \quad \alpha \ge 1 - \frac{2\mu(\mathbf{x}) + 3\lambda(\mathbf{x})}{2\mu^0 + \lambda^0} \quad \text{or} \quad \alpha \ge 1 - \frac{2\mu(\mathbf{x})}{\mu^0}. \qquad (79) $$

Then, the following conditions,

$$ 2\mu^0 + \lambda^0 > \mu(\mathbf{x}) \quad \text{and} \quad 2\mu^0 + \lambda^0 > \mu(\mathbf{x}) + \frac{3}{2}\lambda(\mathbf{x}) \quad \text{and} \quad \mu^0 > \mu(\mathbf{x}), \qquad (80) $$

enforce $\alpha > -1$, and therefore ensure the convergence of the iterative scheme.

Following [39], one might try to estimate the rate of convergence thanks to the spectral radius of the operator, i.e. the minimum and maximum values of the eigenvalues $\alpha$. An eigenvalue $\alpha$ can be expressed with respect to $e_i$, $i = 1, 2, 3, 4$ using (78). It reads

$$ \alpha = 1 + \frac{e_1(\mathbf{x}) - 2\mu(\mathbf{x})}{2\mu_0 + \lambda_0} = 1 + \frac{e_2(\mathbf{x}) - (2\mu(\mathbf{x}) + 3\lambda(\mathbf{x}))}{2\mu_0 + \lambda_0} = 1 + \frac{e_3(\mathbf{x}) - 2\mu(\mathbf{x})}{\mu_0} = 1 + \frac{e_4(\mathbf{x})}{\mu_0}, \quad \forall \mathbf{x}. \qquad (81) $$

Let us denote as $\alpha_i(\mathbf{x})$ the following expressions, obtained by putting $e_i(\mathbf{x})$ to zero in (81),

$$ 
\left\{ 
\begin{aligned}
\alpha_1(\mathbf{x}) &= 1 - \frac{2\mu(\mathbf{x})}{2\mu_0 + \lambda_0}, \\
\alpha_2(\mathbf{x}) &= 1 - \frac{2\mu(\mathbf{x}) + 3\lambda(\mathbf{x})}{2\mu_0 + \lambda_0}, \\
\alpha_3(\mathbf{x}) &= 1 - \frac{2\mu(\mathbf{x})}{\mu_0}, \\
\alpha_4(\mathbf{x}) &= 1.
\end{aligned}
\right. 
\qquad (82) 
$$

Considering that, for any eigenvalue $\alpha$, there are $i$ and $x$ such that $e_i(\mathbf{x}) \ge 0$ and $i$ and $x$ such that $e_i(\mathbf{x}) \le 0$, equation (81) leads to the following bounds

$$ \min_{i,x} \alpha_i(\mathbf{x}) \le \alpha \le \max_{i,x} \alpha_i(\mathbf{x}) = 1. \qquad (83) $$

Thus, this approach only proves that the spectral radius is inferior to one, which ensures convergence of the iterative scheme (at least on the finite dimension spaces considered in the practical applications), but does not give any further information about the rate of convergence. Considering the classical scheme and the proof of convergence proposed by Michel et al. [39], in which the spectral radius is bounded by values related to the comparison material and strictly inferior to one, it is interesting to note that the modification of the auxiliary problem leads to different eigenvalues of the local operator (matrix (76)), and particularly to the $e_4$ eigenvalue, whose eigenvector is related to the antisymmetric part of $\nabla u$, not considered in the classical framework. This particularity prevents from a better bound of the spectral radius, and therefore, contrary to the periodic case considered in [39], it does not provide an estimate of the rate of convergence. Further numerical assessments will give a better appreciation of the actual rate of convergence of the proposed method.

# 4. Applications

## 4.1. Description of the simulations

The discrete sine-cosine based method is now applied to the study of the local and overall response of a composite material made of two isotropic phases,<fnref n="6" /> a matrix of elasticity properties ($\kappa_1, \mu_1$) and an inclusion of elasticity properties ($\kappa_2, \mu_2$). In the applications considered, we fix Poisson's ratios in each phase to $\nu_1 = \nu_2 = 0.25$, so the contrast between the phases reduce to a single parameter and is defined by $\kappa_2/\kappa_1 = \mu_2/\mu_1 = \lambda_2/\lambda_1 = E_2/E_1$. For the application, we consider a 2-d plane strain case and the 2-d domain $[0, 1] \times [0, 1]$ is discretized with 512 x 512 grid points, comprising a centered square inclusion of side 0.5. Following the convergence conditions given by (80), we consider the values $\lambda_0 = 1.5 \times \max(\lambda_1, \lambda_2)$ and $\mu_0 = 1.5 \times \max(\mu_1, \mu_2)$ for the homogeneous reference material used to define the auxiliary problem.

Several types of boundary conditions will be considered:

* Kinematic Uniform Boundary Conditions (KUBC) (corresponding to a particular Dirichlet case) associated with the volume average (macroscopic) of the strain field $\langle \epsilon \rangle$, that is

$$ \mathbf{u}(\mathbf{x}) = \langle \epsilon \rangle \cdot \mathbf{x}, \quad \forall \mathbf{x} \in \partial\Omega. \qquad (84) $$

<footnote n="6">⁶ We consider, for illustrative purposes, an isotropic behavior, but the method is naturally written in the anisotropic setting as shown by Eq. (32).</footnote>
<figure><img src="image_5.png" /><figcaption>Fig. 1. Distribution of the components of the displacement field in the KUBC case. (a) Component $u_1$ and (b) Component $u_2$.</figcaption></figure>

• Static Uniform Boundary Conditions (SUBC) (corresponding to a particular Neumann case) associated with the volume average of the stress field $\langle \sigma \rangle$, that is

$$ \sigma(\mathbf{x}) \cdot \mathbf{n}(\mathbf{x}) = \langle \sigma \rangle \cdot \mathbf{n}(\mathbf{x}), \quad \forall \mathbf{x} \in \partial \Omega, \qquad (85) $$

where $\mathbf{n}(\mathbf{x})$ is the outward normal.

• Arbitrary Dirichlet–Neumann Boundary Conditions (ADNBC), for which an arbitrary displacement $\mathbf{u}(\mathbf{x})$ is prescribed on the boundary $\partial \Omega_D$ and a surface force $\sigma(\mathbf{x}) \cdot \mathbf{n}(\mathbf{x})$ is prescribed on the boundary $\partial \Omega_N$.

The cases of KUBC and SUBC will be used to calculate the macroscopic (effective) elasticity tensor using the relation

$$ \langle \sigma \rangle = \bar{C} : \langle \varepsilon \rangle, \qquad (86) $$

where $\langle \sigma \rangle$ is calculated as the volume average of the stress $\sigma$ (in the KUBC case) and $\langle \varepsilon \rangle$ is calculated as the volume average of the strain $\varepsilon$ (in the SUBC case). In the present 2D case, $\bar{C}$ will be expressed using the following Kelvin notation:

$$ \begin{pmatrix} \langle \sigma_{11} \rangle \\ \langle \sigma_{22} \rangle \\ \langle \sigma_{12} \rangle \end{pmatrix} = \begin{pmatrix} \bar{C}_{1111} & \bar{C}_{1122} & \bar{C}_{1112} \\ \bar{C}_{1122} & \bar{C}_{2222} & \bar{C}_{2212} \\ \bar{C}_{1112} & \bar{C}_{2212} & \bar{C}_{1212} \end{pmatrix} \cdot \begin{pmatrix} \langle \varepsilon_{11} \rangle \\ \langle \varepsilon_{22} \rangle \\ 2 \langle \varepsilon_{12} \rangle \end{pmatrix}. \qquad (87) $$

Given the symmetries of the problem considered, the overall elasticity tensor is expected to be cubic (i.e. $\bar{C}_{1112} = \bar{C}_{2212} = 0$ and $\bar{C}_{1111} = \bar{C}_{2222}$).

In all the applications considered, we take the value tol = 10⁻⁸ for the tolerance used in criterion (62).

## 4.2. Local fields for KUBC and SUBC cases

We consider a contrast $\kappa_2/\kappa_1 = 10^2$.

**Kinematic uniform boundary conditions.** We perform a KUBC simulation with the following loading

$$ \langle \varepsilon \rangle = \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}. \qquad (88) $$

Convergence of the iterative scheme is reached for about 2000 iterations. The distribution of the components $u_1$ and $u_2$ of the displacement field is represented in Fig. 1. One can notice that the boundary conditions (84) are correctly applied, for the loading considered (88). The components $\varepsilon_{11}$ and $\varepsilon_{12}$ are represented in Fig. 2. These fields show no numerical artifact at the cell boundaries but classical oscillations are observed near the inclusion corner, which is typical when a continuous Green operator is used. The components $\sigma_{11}$ and $\sigma_{12}$ are finally represented in Fig. 2.

Considering the contrast $\kappa_2/\kappa_1 = 10^2$, the calculations performed using KUBC lead to the macroscopic elasticity tensor (expressed using Eq. (87)):

$$ \bar{C}_{\text{KUBC}} = \begin{pmatrix} 1.93 & 0.51 & 0 \\ 0.51 & 1.93 & 0 \\ 0 & 0 & 0.63 \end{pmatrix}. \qquad (89) $$
<figure><img src="image_6.png" /><figcaption>Fig. 2. Distribution of the components of the strain and stress fields in the KUBC case. (a) Strain component $\epsilon_{11}$, (b) Strain component $\epsilon_{12}$, (c) Stress component $\sigma_{11}$ and (d) Stress component $\sigma_{12}$.</figcaption></figure>

**Static uniform boundary conditions.** We continue with a SUBC simulation considering the following loading

$$ \langle \sigma \rangle = \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}. \qquad (90) $$

In that case, convergence of the iterative scheme is reached for about 6500 iterations. The components $\epsilon_{11}$ and $\epsilon_{12}$ are represented in Fig. 3. As previously, no numerical artifact are observed at the cell boundaries but classical oscillations are observed near the inclusion corner. Finally the components $\sigma_{11}$ and $\sigma_{12}$ are represented in Fig. 3. One can notice that the boundary conditions (85) are correctly applied, for the loading considered (90).

The calculations performed using SUBC lead to the macroscopic elasticity tensor (expressed using the convention of Eq. (87)):

$$ \bar{C}_{\text{SUBC}} = \begin{pmatrix} 1.81 & 0.59 & 0 \\ 0.59 & 1.81 & 0 \\ 0 & 0 & 0.54 \end{pmatrix}. \qquad (91) $$

It is easy to show that $\bar{C}_{\text{KUBC}} > \bar{C}_{\text{SUBC}}$ in the sense of quadratic forms.

### 4.3. Influence of the contrast

We study in this section the influence of the contrast upon (i) the number of iteration to convergence and (ii) the macroscopic elastic properties. We focus on the 2-d plane strain bulk modulus only, which is defined as

$$ \bar{\kappa}^{2d} = \frac{\bar{C}_{1111} + \bar{C}_{1122}}{2}. \qquad (92) $$
<figure><img src="image_7.png" /><figcaption>Fig. 3. Distribution of the components of the strain and stress fields in the SUBC case. (a) Strain component $\varepsilon_{11}$, (b) Strain component $\varepsilon_{12}$, (c) Stress component $\sigma_{11}$ and (d) Stress component $\sigma_{12}$.</figcaption></figure>

In addition to the KUBC and SUBC cases, we also consider the periodic case (PBC) in order to provide a numerical comparison of the number of iterations to convergence between the periodic and the non-periodic approaches (and also the macroscopic 2-d bulk modulus).

First, number of iterations up to convergence is represented in Fig. 4 as a function of the contrast. As expected with a continuous Green operator with a fixed-point iterative procedure, the number of iterations increases with the contrast, similarly to the non-periodic conductivity case [29]. One might remark that, despite the poor spectral radius bound derived in Section 3.5 (which is not related to the contrast) one retrieves the classical linear dependency on the contrast, as with the classical scheme (PBC). Accelerated schemes developed in the periodic case are expected to improve the convergence rate (see e.g. [6]). It must be noted that the SUBC case requires extra iterations due to the one-voxel layer that is required for imposing the surface force; this could be possibly improved by adjusting from one iteration to another (i) the material of the layer and (ii) the prescribed displacement field so as to avoid any strong discontinuities across the interface.

Then, the dependency of the (normalized) macroscopic bulk modulus $\overline{\kappa}^{2d}/\kappa_1^{2d}$ with the contrast is represented in Fig. 4b. As expected it increases with the contrast and it is bounded by the asymptotic values $\overline{\kappa}^{2d}/\kappa_1^{2d} = 0.47$ and $\overline{\kappa}^{2d}/\kappa_1^{2d} = 0.4$ when $\kappa_2^{2d}/\kappa_1^{2d} \to 0$, respectively for the KUBC and SUBC cases, and $\overline{\kappa}^{2d}/\kappa_1^{2d} = 1.53$ and $\overline{\kappa}^{2d}/\kappa_1^{2d} = 1.52$ when $\kappa_2^{2d}/\kappa_1^{2d} \to \infty$, respectively for the KUBC and SUBC cases. The value predicted using PBC lies between the bounds obtained using KUBC and SUBC; this is related to the (standard) property $\overline{C}_{\text{KUBC}} > \overline{C}_{\text{PBC}} > \overline{C}_{\text{SUBC}}$ in the sense of quadratic forms.

### 4.4. Tensile test of square heterogeneous beam

We finally study the case of a square heterogeneous “beam” subjected to a tensile test (see Fig. 5). This case is interesting as it shows the ability of the present approach to be used not only in the context of micromechanics (using KUBC or SUBC) but also for beams (or plates) made of heterogeneous constituents, which can be considered as a first step towards the used of FFT-based for simple structures. The boundary conditions, shown in Fig. 5, can be applied by considering
<figure><img src="image_8.png" /><figcaption>Fig. 4. Influence of the contrast. (a) Number of iterations up to convergence and (b) Normalized macroscopic bulk modulus $\overline{\kappa}^{2d}/\kappa_1^{2d}$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 5. Tensile test of square heterogeneous “beam”.</figcaption></figure>

• For the displacement $u_x$, Dirichlet–Neumann conditions on the x–faces and Neumann–Neumann conditions on the y–faces. This leads to a DST-III for the x–direction and DCT-I for the y–direction;

• For the displacement $u_y$, Neumann–Neumann conditions on the x–faces and Neumann–Neumann conditions on the y–faces. This leads to a DCT-I for the x–direction and DCT-I for the y–direction;

It must be noted that a clamped boundary condition could also be applied by considering, for $u_y$, DST-III instead of DCT-I. We consider the case $v_1 = v_2 = 0.25$ and $E_2/E_1 = 1/10$, which corresponds to a soft inclusion.

The distribution of the displacement field is first represented in Fig. 6. One can remark that the boundary conditions are correctly applied as the component $u_1$ is null on the left boundary. The component $u_1$ is heterogeneous but it is overall increasing with x which is expected from a tensile test. One can remark that $u_1$ is heterogeneous on the right boundary ($x = L_x$), as there is no prescribed displacement on this face. The component $u_2$ is also heterogeneous and with negative and positive values respectively on the top and bottom boundaries, which classically corresponds to a Poisson effect.
<figure><img src="image_10.png" /><figcaption>Fig. 6. Distribution of the components of the displacement field in the “tensile test”. (a) Component $u_1$ and (b) Component $u_2$.</figcaption></figure>

The components of the strain ($\epsilon_{11}$ and $\epsilon_{22}$) and stress ($\sigma_{11}$ and $\sigma_{22}$) fields are then represented in Fig. 7. The component $\sigma_{11}$ is uniform on the right boundary since a uniform surface force has been applied. The presence of the inclusion leads to heterogeneous mechanical fields. In particular, the tensile stress is heterogeneous in the left boundary which is due to the inclusion. Overall, those results enhance the importance of the boundary conditions for this kind of situation, as the presence of the inclusion results in significantly different displacement, strain and stress components on opposite faces, which cannot be predicted by the usual periodic boundary conditions.

## 5. Conclusion

This work was devoted to the development of a new class of spectral FFT-based solvers for materials with heterogeneous mechanical properties and subjected to non-periodic, Dirichlet and/or Neumann, boundary conditions. This works is thus a direct extension of the seminal work of Moulinec and Suquet [1] to non-periodic boundary conditions. The method is based, as in the case of conductivity [27–29], on the decomposition of the displacement field solution of the problem into a known term verifying the boundary conditions and a fluctuation term described by appropriate sine–cosine series inducing no contribution on the boundary. The principle of the method is then very similar to the periodic case, i.e. it is based on an iterative solution of the so-called Lippmann–Schwinger equation [1]. The auxiliary problem is simply modified (through some reference material acting on the displacement gradient and not its symmetrical part) in order not to include cross-derivatives (which were shown to lead to a hard-to-solve linear system). This auxiliary problem is solved using discrete sine–cosine transforms which naturally emerge in the Galerkin method using an approximation space spanned by sine–cosine series. The method therefore relies on the numerical complexity of fast Fourier transforms and allows the consideration of arbitrary boundary conditions on entire faces, including as particular cases kinematic uniform boundary conditions and static uniform boundary conditions. The method has been successfully applied to several cases including the homogenization of a composite (using both KUBC and SUBC) as well as some “structure” (a beam subjected to a tensile loading).

As shown by a study of the influence of the contrast, the present method suffers from the same drawbacks than the initial method of Moulinec and Suquet [1] in the periodic case, i.e. the convergence rate scales linearly with the contrast, making the method not suitable for highly-contrasted materials. This can be improved by the use of a discrete scheme based on finite differences as done by Gélébart [27] in conductivity problems; this will be presented in a future work. This approach based on discrete Green operator is also expected to improve the quality of the solution field, with a notable reduction of spurious oscillations [30]. In addition, accelerated iterative schemes are also expected to improve the convergence of the method [3–7].

This work was restricted to linear elasticity problems for illustrative purposes. The extension of the present method to material non-linearity (e.g. plasticity, damage, etc.) is straightforward [1]: it simply requires to solve an additional *local* problem, as the method itself corresponds to the resolution of the *global step* (i.e. the global equilibrium of the structure).

This work paves the way of using the FFT-based framework of Moulinec and Suquet [1] in applications that are beyond periodic homogenization, including notably inertia effects in elastodynamic problems of heterogeneous materials, higher-order homogenization (of architectured materials) using higher-order boundary type conditions such as quadratic boundary conditions (QBC), and crack propagation in (finite) heterogeneous microstructures (without numerical artifacts in the propagation due to periodicity [40]).
<figure><img src="image_11.png" /><figcaption>Fig. 7. Distribution of the components of the strain and stress fields in the “tensile test”. (a) Strain component $ε_{11}$, (b) Strain component $ε_{22}$, (c) Stress component $σ_{11}$ and (d) Stress component $σ_{22}$.</figcaption></figure>

CRediT authorship contribution statement

**Joseph Paux:** Writing – review & editing, Writing – original draft, Supervision, Software, Methodology, Investigation, Conceptualization.
**Léo Morin:** Writing – review & editing, Writing – original draft, Supervision, Software, Methodology, Investigation, Conceptualization.
**Lionel Gélébart:** Writing – review & editing, Writing – original draft, Supervision, Software, Methodology, Investigation, Conceptualization.
**Abdoul Magid Amadou Sanoko:** Investigation.

Declaration of competing interest

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

Appendix A. Orthogonality of the modes

The weak formulation (38) involves the calculation of elementary integrals which can be of two types, *identical modes* (associated to second order direct partial derivatives) and *cross modes* (associated to cross partial derivatives).
Considering $k_1$ and $k_2$ two integers, the elementary integrals associated to identical mode read

$$ egin{cases} 
\frac{2}{L} \int_{0}^{L} \sin\left(k_{1} \frac{x \pi}{L}\right) \sin\left(k_{2} \frac{x \pi}{L}\right) dx = \begin{cases} 0 & \text{if } k_{1} \neq k_{2} \\ 1 & \text{if } k_{1} = k_{2} \end{cases} \\ \\ 
\frac{2}{L} \int_{0}^{L} \cos\left(k_{1} \frac{x \pi}{L}\right) \cos\left(k_{2} \frac{x \pi}{L}\right) dx = \begin{cases} 0 & \text{if } k_{1} \neq k_{2} \\ 1 & \text{if } k_{1} = k_{2} \neq 0 \\ 2 & \text{if } k_{1} = k_{2} = 0 \end{cases} \\ \\ 
\frac{2}{L} \int_{0}^{L} \sin\left((2k_{1} + 1) \frac{x \pi}{2L}\right) \sin\left((2k_{2} + 1) \frac{x \pi}{2L}\right) dx = \begin{cases} 0 & \text{if } k_{1} \neq k_{2} \\ 1 & \text{if } k_{1} = k_{2} \end{cases} \\ \\ 
\frac{2}{L} \int_{0}^{L} \cos\left((2k_{1} + 1) \frac{x \pi}{2L}\right) \cos\left((2k_{2} + 1) \frac{x \pi}{2L}\right) dx = \begin{cases} 0 & \text{if } k_{1} \neq k_{2} \\ 1 & \text{if } k_{1} = k_{2} \end{cases} 
\end{cases} \quad (A.1) $$

Considering again $k_1$ and $k_2$ two integers, the elementary integrals associated to cross modes read

$$ egin{cases} 
\frac{2}{L} \int_{0}^{L} \cos\left(k_{1} \frac{x \pi}{L}\right) \sin\left(k_{2} \frac{x \pi}{L}\right) dx = \begin{cases} \frac{2k_{2}}{\pi} \frac{(-1)^{k_{1}+k_{2}} - 1}{k_{1}^{2} - k_{2}^{2}} & \text{if } k_{1} \neq k_{2} \\ 0 & \text{if } k_{1} = k_{2} \end{cases} \\ \\ 
\frac{2}{L} \int_{0}^{L} \sin\left((2k_{1} + 1) \frac{x \pi}{2L}\right) \sin\left(k_{2} \frac{x \pi}{L}\right) dx = \frac{8}{\pi} \frac{(-1)^{k_{1}+k_{2}} k_{2}}{(1 + 2k_{1})^{2} - 4k_{2}^{2}} \\ \\ 
\frac{2}{L} \int_{0}^{L} \cos\left((2k_{1} + 1) \frac{x \pi}{2L}\right) \sin\left(k_{2} \frac{x \pi}{L}\right) dx = -\frac{8}{\pi} \frac{k_{2}}{(1 + 2k_{1})^{2} - 4k_{2}^{2}} \\ \\ 
\frac{2}{L} \int_{0}^{L} \cos\left((2k_{1} + 1) \frac{x \pi}{2L}\right) \cos\left(k_{2} \frac{x \pi}{L}\right) dx = \frac{4}{\pi} \frac{(-1)^{k_{1}+k_{2}} (1 + 2k_{1})}{(1 + 2k_{1})^{2} - 4k_{2}^{2}} \\ \\ 
\frac{2}{L} \int_{0}^{L} \sin\left((2k_{1} + 1) \frac{x \pi}{2L}\right) \cos\left(k_{2} \frac{x \pi}{L}\right) dx = \frac{4}{\pi} \frac{1 + 2k_{1}}{(1 + 2k_{1})^{2} - 4k_{2}^{2}} \\ \\ 
\frac{2}{L} \int_{0}^{L} \sin\left((2k_{1} + 1) \frac{x \pi}{2L}\right) \cos\left((2k_{2} + 1) \frac{x \pi}{2L}\right) dx = \begin{cases} \frac{1 + 2k_{1} + (-1)^{k_{1}+k_{2}} (1 + 2k_{2})}{\pi(k_{1} - k_{2})(1 + k_{1} + k_{2})} & \text{if } k_{1} \neq k_{2} \\ \frac{2}{\pi} \frac{1}{1 + 2k_{1}} & \text{if } k_{1} = k_{2} \end{cases} 
\end{cases} \quad (A.2) $$

## Appendix B. Detailed calculation of the elementary integrals

The integrals $I_{lmn}^x$, $K_{lmn}^x$ and $I_{lmn}^{x0}$ involved in the weak formulation of the modified auxiliary problem are calculated as follows.

* The integral $I_{lmn}^x$ reads

$$ I_{lmn}^{x} = \sum_{i=0}^{N_x+1} \sum_{j=0}^{N_y+1} \sum_{k=0}^{N_z+1} \left[ - \left( (\lambda_0 + 2\mu_0)(\xi_l^{xy})^2 + \mu_0 \left( (\xi_j^{xy})^2 + (\xi_k^{xz})^2 \right) \right) \alpha_l^{xx} \alpha_j^{xy} \alpha_k^{xz} U_{ijk}^x \times \int_{x=0}^{x=L_x} g^{xx}(k_l^{xx} x) g^{xx}(k_j^{xx} x) dx \times \int_{y=0}^{y=L_y} g^{xy}(k_j^{xy} y) g^{xy}(k_m^{xy} y) dy \times \int_{z=0}^{z=L_z} g^{xz}(k_k^{xz} z) g^{xz}(k_n^{xz} z) dz \right] \quad (B.1) $$

After calculation of the elementary integrals in Eq. (B.1), it can be shown that

$$ I_{lmn}^{x} = -\frac{L_x L_y L_z}{8} \left( (\lambda_0 + 2\mu_0)(\xi_l^{xy})^2 + \mu_0 \left( (\xi_m^{xy})^2 + (\xi_n^{xz})^2 \right) \right) U_{lmn}^{x}, \quad (B.2) $$

since

$$ \int_{x=0}^{x=L_x} g^{xx}(k_l^{xx} x) g^{xx}(k_j^{xx} x) dx = \frac{\delta_{il} L_x}{2 \alpha_l^{xx}}, \quad \int_{y=0}^{y=L_y} g^{xy}(k_j^{xy} y) g^{xy}(k_m^{xy} y) dy = \frac{\delta_{jm} L_y}{2 \alpha_m^{xy}}, \quad \int_{z=0}^{z=L_z} g^{xz}(k_k^{xz} z) g^{xz}(k_n^{xz} z) dz = \frac{\delta_{kn} L_z}{2 \alpha_n^{xz}}, \quad (B.3) $$

with $\delta$ the Kronecker symbol.

* The integral $I_{lmn}^{x0}$ reads

$$ I_{lmn}^{x0} = \int_{\Omega} \left( (\lambda_0 + 2\mu_0) \frac{\partial^2 u_x^0}{\partial x^2}(\mathbf{x}) + \mu_0 \left( \frac{\partial^2 u_x^0}{\partial y^2}(\mathbf{x}) + \frac{\partial^2 u_x^0}{\partial z^2}(\mathbf{x}) \right) \right) g^{xx}(k_l^{xx} x) g^{xy}(k_m^{xy} y) g^{xz}(k_n^{xz} z) d\Omega. \quad (B.4) $$

Assuming that all second derivatives of the function $u_x^0$ are known analytically, then the integral (B.4) can be simply calculated approximately as

$$ I_{lmn}^{x0} = \frac{L_x L_y L_z}{(N_x + 1)(N_y + 1)(N_z + 1)} (D_{xyz}^{x}(\Delta u_x^0))_{lmn}, \quad (B.5) $$
where $\Delta u_x^0$ is the array of size $(N_x+2) \times (N_y+2) \times (N_z+2)$ containing the values of the function $(\lambda_0+2\mu_0) \frac{\partial^2 u_x^0}{\partial x^2} + \mu_0 \left( \frac{\partial^2 u_x^0}{\partial y^2} + \frac{\partial^2 u_x^0}{\partial z^2} \right)$ at the grid points. It must be noted that this only requires the knowledge of the derivatives of $u_x^0$ at the grid points; consequently, if the analytical expression of the derivatives of the function $u_x^0$ is not known, the integral $I_{lmn}^{x0}$ can still be determined using an approximation of its derivatives (using finite differences for instance).

* The integral $K_{lmn}^x$ reads

$$K_{lmn}^x = K_{lmn}^{xx} + K_{lmn}^{xy} + K_{lmn}^{xz} \quad (B.6)$$

with

$$K_{lmn}^{xx} = \int_{\Omega} \frac{\partial \tau_{xx}}{\partial x}(\mathbf{x}) g^{xx}(k_l^{xx} \mathbf{x}) g^{xy}(k_m^{xy} \mathbf{y}) g^{xz}(k_n^{xz} \mathbf{z}) d\Omega, \quad (B.7)$$

$$K_{lmn}^{xy} = \int_{\Omega} \frac{\partial \tau_{xy}}{\partial y}(\mathbf{x}) g^{yx}(k_l^{yx} \mathbf{x}) g^{yy}(k_m^{yy} \mathbf{y}) g^{yz}(k_n^{yz} \mathbf{z}) d\Omega, \quad (B.8)$$

$$K_{lmn}^{xz} = \int_{\Omega} \frac{\partial \tau_{xz}}{\partial z}(\mathbf{x}) g^{zx}(k_l^{zx} \mathbf{x}) g^{zy}(k_m^{zy} \mathbf{y}) g^{zz}(k_n^{zz} \mathbf{z}) d\Omega. \quad (B.9)$$

Following [29], the integrals defining $K_{lmn}^x$ are calculated using integration by parts since, in the definition of the auxiliary problem, the known fields are $\tau_{xx}$, $\tau_{xy}$ and $\tau_{xz}$, and not their partial derivatives. We only detail hereafter the calculation of $K_{lmn}^{xx}$:

$$K_{lmn}^{xx} = \int_{y=0}^{y=L_y} \int_{z=0}^{z=L_z} \left( g^{xx}(k_l^{xx} L_x) \tau_{xx}(L_x, y, z) - g^{xx}(0) \tau_{xx}(0, y, z) \right) g^{xy}(k_m^{xy} y) g^{xz}(k_n^{xz} z) dy dz \\ - \xi_l^{xx} \int_{\Omega} \tau_{xx}(x, y, z) g^{xx}(k_l^{xx} x) g^{xy}(k_m^{xy} y) g^{xz}(k_n^{xz} z) d\Omega. \quad (B.10)$$

The amount of calculations is reduced by the introduction of the function $\tilde{\tau}_{xx}$, whose expression depends on the type of boundary conditions associated to $u_x$ on the x-faces:

$$ \left\{ \begin{array}{l} \tilde{\tau}_{xx}(x, y, z) = \tau_{xx}(x, y, z) \hfill (\text{DD}) \\ \tilde{\tau}_{xx}(x, y, z) = \tau_{xx}(x, y, z) - \tau_{xx}(0, y, z) - \frac{\tau_{xx}(L_x, y, z) - \tau_{xx}(0, y, z)}{L_x} x \hfill (\text{NN}) \\ \tilde{\tau}_{xx}(x, y, z) = \tau_{xx}(x, y, z) - \tau_{xx}(L_x, y, z) \hfill (\text{DN}) \\ \tilde{\tau}_{xx}(x, y, z) = \tau_{xx}(x, y, z) - \tau_{xx}(0, y, z) \hfill (\text{ND}) \end{array} \right. \quad (B.11) $$

The approximate calculation of $K_{lmn}^{xx}$, using discrete sine-cosine transforms, leads to

$$K_{lmn}^{xx} = \frac{-L_x L_y L_z}{(N_x + 1)(N_y + 1)(N_z + 1)} \left[ \xi_l^{xx} (D_{xyz}^x (\tilde{\tau}_x))_{lmn} + S_{lmn}^{xx} \right], \quad (B.12)$$

where $\tilde{\tau}_{xx}$ denotes the array of size $(N_x + 2) \times (N_y + 2) \times (N_z + 2)$ containing the values of the function $\tilde{\tau}_{xx}$ at the grid points. The term $S_{lmn}^{xx}$ is always null, expect in the case $l=0$ of the NN boundary conditions on the x-component of $u_x$, where it is given by

$$S_{0mn}^{xx} = - \frac{(N_x + 1)}{L_x} (D_{yz}^x (\tau_{xx}^0 - \tau_{xx}^{L_x}))_{mn}. \quad (B.13)$$

The term $K_{lmn}^x$ finally reads

$$K_{lmn}^{x} = \frac{-L_x L_y L_z}{(N_x + 1)(N_y + 1)(N_z + 1)} \left[ \xi_l^{xx} (D_{xyz}^x (\tilde{\tau}_{xx}))_{lmn} + \xi_m^{xy} (D_{xyz}^x (\tilde{\tau}_{xy}))_{lmn} + \xi_n^{xz} (D_{xyz}^x (\tilde{\tau}_{xz}))_{lmn} + S_{lmn}^{x} \right] \quad (B.14)$$

where $\tilde{\tau}_{xx}$, $\tilde{\tau}_{xy}$ and $\tilde{\tau}_{xz}$ respectively denote the arrays of size $(N_x + 2) \times (N_y + 2) \times (N_z + 2)$ containing the values of functions $\tilde{\tau}_{xx}$, $\tilde{\tau}_{xy}$ and $\tilde{\tau}_{xz}$ at the grid points, and $S_{lmn}^x$ is given by

$$S_{lmn}^{x} = 0 \quad \text{if } l \neq 0 \text{ and } m \neq 0 \text{ and } n \neq 0 \quad (B.15)$$

$$S_{0mn}^{x} = - \frac{(N_x + 1)}{L_x} (D_{yz}^{x} (\tau_{xx}^{0} - \tau_{xx}^{L_x}))_{mn} \quad \text{if NN on } x \text{-faces}$$

$$S_{l0n}^{x} = - \frac{(N_y + 1)}{L_y} (D_{xz}^{x} (\tau_{xy}^{0} - \tau_{xy}^{L_y}))_{ln} \quad \text{if NN on } y \text{-faces}$$

$$S_{lm0}^{x} = - \frac{(N_z + 1)}{L_z} (D_{xy}^{x} (\tau_{xz}^{0} - \tau_{xz}^{L_z}))_{lm} \quad \text{if NN on } z \text{-faces,} \quad (B.16)$$

where $\tau_{xx}^0, \tau_{xx}^{L_x}, \tau_{xy}^0, \tau_{xy}^{L_y}, \tau_{xz}^0$ and $\tau_{xz}^{L_z}$ denote the arrays of sizes $(N_y + 2) \times (N_z + 2)$ for $\tau_{xx}^0$ and $\tau_{xx}^{L_x}$, $(N_x + 2) \times (N_z + 2)$ for $\tau_{xy}^0$ and $\tau_{xy}^{L_y}$, and $(N_x + 2) \times (N_y + 2)$ for $\tau_{xz}^0$ and $\tau_{xz}^{L_z}$, containing the values of functions $\tau_{xx}$, $\tau_{xy}$ and $\tau_{xz}$ on the faces.

## Data availability

Data will be made available on request.
References

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

[2] L. Gélébart, Amitex, 2022, https://amitexftfp.github.io/AMITEX/index.html.

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

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

[5] V. Monchiet, G. Bonnet, A polarization-based FFT iterative scheme for computing the effective properties of elastic composites with arbitrary contrast, Internat. J. Numer. Methods Engrg. 89 (2012) 1419-1436.

[6] H. Moulinec, F. Silva, Comparison of three accelerated FFT-based schemes for computing the mechanical response of composite materials, Internat. J. Numer. Methods Engrg. 97 (2014) 960-985.

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

[8] R. Lebensohn, N-site modeling of a 3D viscoplastic polycrystal using Fast Fourier Transform, Acta Mater. 49 (2001) 2723-2737.

[9] R. Brenner, A.J. Beaudoin, P. Suquet, A. Acharya, Numerical implementation of static field dislocation mechanics theory for periodic media, Phil. Mag. 94 (2014) 1764-1787.

[10] N. Bertin, M.V. Upadhyay, C. Pradalier, L. Capolungo, A FFT-based formulation for efficient mechanical fields computation in isotropic and anisotropic periodic discrete dislocation dynamics, Modelling Simul. Mater. Sci. Eng. 23 (2015) 065009.

[11] R. Brenner, Numerical computation of the response of piezoelectric composites using Fourier transform, Phys. Rev. B 79 (2009) 184106, Publisher: American Physical Society.

[12] N. Bilger, F. Auslender, M. Bornert, J.C. Michel, H. Moulinec, P. Suquet, A. Zaoui, Effect of a nonuniform distribution of voids on the plastic response of voided materials: a computational and statistical analysis, Int. J. Solids Struct. 42 (2005) 517-538.

[13] J. Paux, R. Brenner, D. Kondo, Plastic yield criterion and hardening of porous single crystals, Int. J. Solids Struct. 132-133 (2018) 80-95.

[14] M. Schneider, A review of nonlinear FFT-based computational homogenization methods, Acta Mech. 232 (2021) 2051-2100.

[15] S. Lucarini, M.V. Upadhyay, J. Segurado, FFT based approaches in micromechanics: fundamentals, methods and applications, Model. Simul. Mater. Sci. Eng. 30 (2021) 023002, Publisher: IOP Publishing.

[16] L. Gélébart, A simple extension of FFT-based methods to strain gradient loadings-application to the homogenization of beams and plates with linear and non-linear behaviors, J. Theoret. Comput. Appl. Mech. (2022).

[17] T.H. Tran, V. Monchiet, G. Bonnet, A micromechanics-based approach for the derivation of constitutive elastic coefficients of strain-gradient media, Int. J. Solids Struct. 49 (2012) 783-792.

[18] V. Kouznetsova, M.G. Geers, W. Brekelmans, Multi-scale second-order computational homogenization of multi-phase materials: a nested finite element solution strategy, Comput. Methods Appl. Mech. Eng. 193 (2004) 5525-5550.

[19] J. Yvonnet, N. Auffray, V. Monchiet, Computational second-order homogenization of materials with effective anisotropic strain-gradient behavior, Int. J. Solids Struct. 191-192 (2020) 434-448.

[20] R. Sancho, V. Rey-de Pedraza, P. Lafourcade, R.A. Lebensohn, J. Segurado, An implicit FFT-based method for wave propagation in elastic heterogeneous media, Comput. Methods Appl. Mech. Engrg. 404 (2023) 115772.

[21] L. Gélébart, A modified FFT-based solver for the mechanical simulation of heterogeneous materials with Dirichlet boundary conditions, C. R. Mécanique 348 (2020) 693-704.

[22] N.B. Nkoumbou Kaptchouang, L. Gélébart, Multiscale coupling of FFT-based simulations with the LDC approach, Comput. Methods Appl. Mech. Engrg. 394 (2022) 114921.

[23] H. Grimm-Strele, M. Kabel, FFT-based homogenization with mixed uniform boundary conditions, Internat. J. Numer. Methods Engrg. 122 (2021) 7241-7265.

[24] V. Monchiet, G. Bonnet, FFT based iterative schemes for composite conductors with uniform boundary conditions, Eur. J. Mech. A Solids 103 (2024) 105146.

[25] A. Wiegmann, Fast Elliptic Solvers on Rectangular Parallelepipeds, Technical Report, Lawrence Berkeley National Laboratory, 1999.

[26] L. Risthaus, M. Schneider, Imposing different boundary conditions for thermal computational homogenization problems with FFT- and tensor-train-based Green's operator methods, Internat. J. Numer. Methods Engrg. 125 (2024) e7423.

[27] L. Gélébart, FFT-based simulations of heterogeneous conducting materials with combined non-uniform Neumann, periodic and Dirichlet boundary conditions, Eur. J. Mech. A/Solids (2024) 105248.

[28] L. Morin, J. Paux, A fast numerical method for the conductivity of heterogeneous media with Dirichlet boundary conditions based on discrete sine-cosine transforms, Comput. Methods Appl. Mech. Engrg. 421 (2024) 116772.

[29] J. Paux, L. Morin, L. Gélébart, A discrete sine-cosine transforms galerkin method for the conductivity of heterogeneous materials with mixed Dirichlet/Neumann boundary conditions, Int. J. Numer. Methods Eng. (2024) http://dx.doi.org/10.1002/nme.7615.

[30] F. Willot, B. Abdallah, Y.P. Pellegrini, Fourier-based schemes with modified Green operator for computing the electrical response of heterogeneous media with accurate local fields, Internat. J. Numer. Methods Engrg. 98 (2014) 518-533.

[31] L. Risthaus, M. Schneider, Imposing Dirichlet boundary conditions directly for FFT-based computational micromechanics, Comput. Mech. (2024).

[32] V. Fuka, PoisFFT - A free parallel fast Poisson solver, Appl. Math. Comput. 267 (2015) 356-364.

[33] H. Wang, Y. Zhang, X. Ma, J. Qiu, Y. Liang, An efficient implementation of fourth-order compact finite difference scheme for Poisson equation with Dirichlet boundary conditions, Comput. Math. Appl. 71 (2016) 1843-1860.

[34] F. Caforio, S. Imperiale, A high-order spectral element fast Fourier transform for the Poisson equation, SIAM J. Sci. Comput. 41 (2019) A2747-A2771, Publisher: Society for Industrial and Applied Mathematics.

[35] J. Vondrejc, J. Zeman, I. Marek, An FFT-based Galerkin method for homogenization of periodic media, Comput. Math. Appl. 68 (2014) 156-173.

[36] Z. Wang, Fast discrete sine transform algorithms, Signal Process. 19 (1990) 91-102.

[37] G. Strang, The discrete cosine transform, SIAM Rev. 41 (1999) 135-147.

[38] M. Frigo, S. Johnson, FFTW: an adaptive software architecture for the FFT, in: Proceedings of the 1998 IEEE International Conference on Acoustics, Speech and Signal Processing, ICASSP '98 (Cat. No.98CH36181), vol. 3, 1998, pp. 1381-1384, ISSN: 1520-6149.

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

[40] L. Morin, A. Acharya, Analysis of a model of field crack mechanics for brittle materials, Comput. Methods Appl. Mech. Engrg. 386 (2021) 114061.