# Limitations and merits of wavelet discretizations in FFT-based computational micromechanics

Mohit Dave, Elodie Donval, Matti Schneider

► To cite this version:

Mohit Dave, Elodie Donval, Matti Schneider. Limitations and merits of wavelet discretizations in FFT-based computational micromechanics. Computational Mechanics, 2026, 10.1007/s00466-026-02752-3. hal-05523162

HAL Id: hal-05523162

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

Submitted on 23 Feb 2026

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Limitations and merits of wavelet discretizations in FFT-based computational micromechanics

Mohit Dave<fnref n="1" /> · Elodie Donval<fnref n="2" /> ID · Matti Schneider<fnref n="1" />,<fnref n="3" />,<fnref n="4" /> ID

Received: 30 December 2025 / Accepted: 19 January 2026

© The Author(s) 2026

## Abstract

FFT-based computational micromechanics enables to efficiently analyze the solid mechanical behavior of industrial-scale heterogeneous microstructures. This approach relies on the discretization of the microstructure on a regular grid, which makes it particularly convenient for working with 2D or 3D images. However, efficient local refinement of the grid to capture fine-scale phenomena is not possible within the original framework. Among the various strategies addressing this limitation, Kaiser et al. recently introduced a multiresolution wavelet collocation approach for FFT-based computational micromechanics. They showed the benefits of local adaptivity for examples in one and two dimensions. The present contribution focuses on the properties of the wavelet discretizations in three spatial dimensions. More precisely, we analyze the Deslauriers–Dubuc (DD) wavelets of degrees one to four, connecting and comparing them to established FFT-based homogenization frameworks. Numerical studies on regular grids serve to analyze the convergence behavior of the associated Lippmann-Schwinger solvers, the emerging macroscopic properties, and the resulting local fields obtained with these discretizations on industrial-type microstructures. Also, we study the performance for inelastic constituents. Of particular interest is the result that only first-order DD wavelets ensure stable convergence for porous microstructures.

**Keywords** Wavelet collocation · Computational homogenization · Lippmann-Schwinger equation · Eshelby-Green operator · FFT-based solvers · Convergence behavior

# 1 Introduction

## 1.1 State of the art

Since the pioneering work of Moulinec-Suquet [1, 2], computational methods based on the fast Fourier transform (FFT) have become a standard tool for computing the effective solid mechanical behavior of microstructured materials [3] in an efficient way. Designed for cuboid cells discretized by a regular grid, FFT methods come with a natural implementation

based on the Lippmann-Schwinger equation [4–6]. The latter formulation avoids assembling matrices as well as evaluating material tangents and leads to iteration counts which are bounded independently of the grid spacing, at least for well-behaved material models. We refer to the review articles [7–9] for a broader perspective.

In their original works [1, 2], Moulinec-Suquet provided the whole package: Discretization scheme and solution method were intricately linked, and it took some time for researchers to disentangle them. Indeed, there are two types of convergences involved, the convergence of the solver, i.e., the means to approximate the discretized solution, and the convergence under grid refinement, i.e., the increasing fidelity of the discrete approximations for a finer and finer mesh.

Subsequently, W. H. Müller's group [10–12] suggested a similar scheme as Moulinec-Suquet, but based on a finite-difference discretization. However, this approaches did not permeate the community. Rather, alternative discretizations based on a Fourier-Galerkin discretization were introduced into the FFT-based framework independently by different

<footnote>Matti Schneider
matti.schneider@uni-due.de</footnote>

<footnote n="1">¹ Institute of Engineering Mathematics, University of Duisburg-Essen, Essen, Germany</footnote>

<footnote n="2">² CNRS, Univ. Grenoble Alpes, Grenoble INP, 3SR, Grenoble, France</footnote>

<footnote n="3">³ Center for Nanointegration Duisburg-Essen (CENIDE), Duisburg, Germany</footnote>

<footnote n="4">⁴ Fraunhofer Institute for Industrial Mathematics ITWM, Kaiserslautern, Germany</footnote>
groups [13–15]. Interestingly, these schemes did not improve upon the original Moulinec-Suquet discretization in terms of accuracy or computational speed. However, they may be used to derive upper and lower bounds on the effective stiffness for linear elastic materials. An alternative discretization was introduced by Brisard-Dormieux [16, 17] based on the Hashin-Shtrikman variational principle [18–20] and a discretization by piece-wise homogeneous stress polarization fields [21].

All these cited works considered the Moulinec-Suquet discretization as a *perturbation* of a more conventional discretization. However, usually the "perturbation" performed better than the original. Therefore, more effort was invested in understanding the discretization underlying Moulinec-Suquet's construction. Eventually, the discretization can either be understood as a trigonometric collocation discretization [22] or as a non-conforming Fourier-Galerkin discretization [23], where the non-conformity comes from evaluating the weak form of the balance equation via quadrature, i.e., the trapezoidal rule. With this understanding, guarantees for the convergence under grid refinement could be obtained [17, 23, 24].

At that time, based on earlier numerical strategies exploiting the FFT [25, 26], it became more and more apparent [27, 28] that the discretization scheme and the solver could be chosen independently to some degree. Based on these insights, the computational advantages of the Moulinec-Suquet approach could be preserved, while alleviating some of its shortcomings. To be more precise, the primary interest came in uniting the following three qualities:

(i) Accuracy,

(ii) Efficiency,

(iii) Robustness.

The first item concerns both convergence under grid refinement and countermeasures to local solution artifacts, like the Gibbs' oscillations [29] characteristic for Fourier-type discretizations. Efficiency refers to both memory-efficient implementations and rapid convergence of associated solution schemes. The last point, robustness, is more subtle. As FFT-based methods rely on regular grids, they need to mesh pores and rigid inclusions in microstructures, as well. It was observed [30] that the Moulinec-Suquet discretization does not lead to convergent solvers for the latter scenarios due to an inherent ill-conditioning of the associated discrete systems. Robust discretizations avoid this ill-conditioning [31, 32].

The first generation of improved solvers concerned discretizations based on finite differences [30, 33], finite volumes [34, 35] and finite elements [36–38]. These strategies could reduce the solution artifacts and lead to robust convergence for porous or rigidly reinforced materials, at least

in "stable" cases [31, 32]. Which kind of strategy is used, mostly depends on the personal preference. The finite difference schemes [30, 33] appear to be most used for mechanical problems, as they seem to offer the most "bang for the buck", i.e., they resolve the pertinent issues with minimal implementation and runtime overhead.

Still, research continued on alternative discretization schemes compatible with the Lippmann-Schwinger framework conceived by Moulinec-Suquet. Alternative finite difference schemes [39–43] were investigated, in particular for higher order in space. Such discretizations turn out to be favorable for problems with smooth solutions, like phase-field problems, but may lead to artifacts for discontinuous coefficients and potentially divergence for porous materials [9]. Moreover, B-splines splines [44] were integrated, to a similar effect. A completely different approach was proposed by Derraz et al. [45] who introduced a real-space Lippmann-Schwinger scheme based on the Radon transform. However, their discretization scheme lacks convergence under grid refinement unless specific countermeasures [46] are taken. Recently, Finel introduced the tetrahedral (TET) scheme [47], which may be interpreted as a specific but non-classical finite-difference scheme. The discretization satisfies all three requirements (i)-(iii) discussed above and may thus develop into a serious contender for the established finite difference schemes [30, 33]. However, two material-point evaluations are required per voxel, at least away from the interface [48].

Last but not least, we mention the use of composite-voxel methods [49–51], which enhance voxels containing an interface with a material specific and interface-aware constitutive law to increase the accuracy of the effective material response. Composite voxels apply to all discretizations on a regular grid, however, their performance varies depending on the discretization used.

Almost a decade after finally understanding the Moulinec-Suquet discretization, the basic driver of regular-grid methods was identified: superconvergence. In the context of voxel finite-element methods (FEM), Ye-Chung [52] showed that the effective elastic properties converge with twice the rate of the local strain and stress fields (measured in $L^2$). This observation provided the theoretical underpinnings of the previously observed phenomenon: At the interfaces, the solution quality is low. However, the effective properties turn out to be accurate, essential due to fortuitous cancellations upon averaging. The superconvergence phenomenon could also be verified for inelastic scenarios [53] and could also be used for improving the estimated effective properties [54] via postprocessing. Moreover, the convergence rate of the Moulinec-Suquet discretization for elastic problems was shown to be identical to voxel FEM [55], confirming previous observations by FFT practitioners [30, 36]. Moreover, the restriction to the regular grid was established as the main limitation of FFT-based methods, irrespective of the
discretization. In this phase of understanding, the natural emergence of laminate-based composite voxels [56] and the convergence properties for porous/rigidly reinforced materials [31, 32] should be counted, as well.

Recent innovative work was concerned with either improving the accuracy of FFT-based methods or reducing their effort. For instance, a series of works [57–59] considered omitting certain frequencies when evaluating the Lippmann-Schwinger operator in FFT-based schemes. However, no reduction of material evaluations was considered, the primary driver of computational effort.

Building on the extension of FFT-based methods to finite strains [60], mesh-morphing techniques [61–63] were introduced to better represent the geometry of the microstructure. More precisely, an Eulerian-Lagrangian decomposition is considered, where the interface geometry is resolved in the Eulerian setting and the equations are solved in the Lagrangian configuration. Such techniques promise higher accuracy, but come with an additional ill-conditioning due to severe mesh deformations. An alternative approach is based on the extended finite-element method [64, 65], where a robust FFT solver was introduced recently in the context of thermal conductivity [66].

In the regular grid setting, a first attempt at grid refinement within FFT-based schemes was made by Eyre and Milton [67], who proposed a cascading multigrid solver which prolongs a solution obtained on a coarse grid onto a finer grid, providing an improved initial guess for the solver. An alternative idea was developed by Nkoumbou Kaptchouang and Gélébart [68], who "mortar" a fine grid onto a coarse grid in regions of interest. Kabel-Schneider [69] introduced a coarsening strategy based on octrees to reduce the number of material evaluations in nonlinear schemes.

This work is concerned with the innovative approach put forward by Kaiser et al. [70], who proposed to use interpolating wavelets [71] for discretizing the cell problem of periodic homogenization. They derived the associated Lippmann-Schwinger framework exploiting similarities with finite-difference discretizations, and showed the advantages of using Deslauriers-Dubuc wavelets with four compact support points for one- and two-dimensional computational examples.

The use of wavelets was kick-started with Daubechies' landmark construction of $L^2(\mathbb{R})$-orthonormal wavelets with compact support [72]. Such wavelets permit to inductively decompose an image into *coarse* and *fine* features, permitting to leave out these fine features in case they are small. In particular, this kind of adaptivity may be used for image compression, see Daubechies [73]. In the 1990s, wavelets were considered to be strong candidates for constructing high-order adaptive solvers for partial differential equations (PDEs). In a sense, these hopes could be confirmed, as Cohen et al. [74] proved that wavelets provide the optimal trade-

off between effort and accuracy for solving elliptic PDEs. However, only a few years later, a similar optimality statement could be proved for adaptive FEM [75], putting a lot of work on wavelets for PDEs to a halt. In some sense, working with FEM is much simpler than working with wavelets: The latter require special means to operate on general (non-rectangular) domains, lead with *non-sparse* matrices and struggle with non-linearity. To be more precise, most of the wavelet literature is concerned with wavelets in the Hilbert space $L^2$ of square-integrable functions. Wavelets provide convenient multiscale bases, and to identify the coefficient associated to a specific basis vector, an $L^2$-inner product must be evaluated. Even in case the domain and the boundary conditions are no problem, traditional wavelet approaches have shortcomings, as illustrated by Mehraeen-Chen [76], whose wavelet-Galerkin homogenization approach is restricted to one and two dimensions, and lacks behind comparative FEM approaches. We refer to the review article [77] for further background on the application of wavelet methods to solving PDEs.

Just like spectral Galerkin methods may be inferior to collocation methods, illustrated by the advantages of the Moulinec-Suquet discretization over a comparable Fourier-Galerkin discretization discussed above, wavelet collocation methods may be beneficial, see, for instance Vasilyev-Bowman [78]. Here, the weak form of the PDE is not evaluated exactly, but rather "collocated", i.e., evaluated at specific "quadrature" points. Such strategies, however, lead to non-symmetric systems, also in case the non-discretized variational problem was symmetric. A clever remedy for this problem utilizes interpolating wavelets [71] where the coefficients of the multiscale basis vectors may be determined by point evaluation, i.e., via interpolation. Among the different interpolating wavelet families, the Deslauriers-Dubuc wavelets [79] turn out to be particularly favorable, and an application to mechanical problems was reported by Burgos et al. [80].

Kaiser et al. [70] proposed using Deslauriers-Dubuc wavelets with four support points to discretize the cell problem associated to periodic homogenization, discussed the implementation of an associated Eshelby-Green operator on a regular grid and studied computational examples in one and two spatial dimensions. However, the work of Kaiser et al. [70] does neither concern three-dimensional nor inelastic problems. Moreover, a comparison to the extensive list of competing discretizations is lacking.

## 1.2 Contributions

Based on the previously described limitations of the work by Kaiser et al. [70], the primary research question driving this work is:
What are the properties of the Deslauriers-Dubuc wavelet discretization for unit-cell homogenization problems?

Here, we are mostly concerned with the items (i)-(iii) discussed above. To enable a comparison to other discretizations, we restrict to a regular grid. Indeed, upon activating spatial adaptivity, numerical properties are only expected to become worse.

The setting of periodic homogenization of inelastic homogenization problems is discussed in Appendix A, where also the associated Lippmann-Schwinger framework required for the remainder of the article is summarized. In section 2.1, we discuss the construction of the Deslauriers-Dubuc wavelets on the real line for different orders. We cover both lower and higher orders than considered in Kaiser et al. [70]. We claim no novelty in section 2.1 - however, we wish to point out that the development of interpolating wavelets is non-standard. We side-step the conventional (rather heavy) wavelet machinery and provide a streamlined discussion. As a by-product, we notice that the Deslauriers-Dubuc wavelets with two and four support points, respectively, lead to identical stencils as the standard central difference discretizations of order two and four. As a long-shot consequence, the DD wavelet discretization with two support points for the cell problem leads to identical results as second-order central differences, which were considered before [39, 40]. Moreover, the DD wavelets with four support points proposed by Kaiser et al. [70] are identical to the fourth-order central differences put forward by Vidyasagar et al. [41]. To make this clear: We do not want to belittle the innovative character of using DD wavelets - their main advantage is the local adaptivity. Rather, we wish to emphasize the connections between the available literature and the novel schemes.

Section 2.2 introduces the wavelet-collocation discretization in arbitrary dimensions in a general way. The developments are similar to the Moulinec-Suquet discretization [55] - which we consider an advantage - and different from traditional finite-difference and finite-element discretizations. Indeed, for the latter, there are typically two grids: a "nodal grid" and a "strain grid". Both for the Moulinec-Suquet discretization and the Deslauriers-Dubuc discretization, these two grids coincide, and reside in the voxel centers.

Section 3 is the center of the manuscript: It builds on the master's thesis of the first author [81] and reports a number of computational experiments for industrial-scale microstructures and cases involving infinite material contrast. More precisely, we consider a single spherical inclusion in section 3.2, and compare the DD wavelets with different orders to the Moulinec-Suquet discretization [1, 2], the rotated staggered grid [33] and the conventional staggered grid [30], both for finite and infinite material contrast. In section 3.3, a challenging porous sand-core microstructure of industrial scale is investigated. Last but not least, we are concerned with an

inelastic fiber microstructure in section 3.4. We summarize our findings in the conclusion section 4 and provide perspectives.

## 2 Interpolating wavelets for homogenization problems

### 2.1 Deslauriers-Dubuc wavelets on the real line

The starting point for interpolating wavelets is a given function

$$ \phi: \mathbb{R} \to \mathbb{R}, \qquad (2.1) $$

called the fundamental function by Deslauriers-Dubuc [79] and scaling function in the wavelet literature [73]. The function (2.1) is required to have the following elementary properties:

1. The function (2.1) vanishes at all integers except for zero where the function attains the value unity, i.e., the identity

$$ \phi(k) = \delta_{0k}, \quad k \in \mathbb{Z}, \qquad (2.2) $$

holds with Dirac's symbol $\delta_{jk}$.

2. The function (2.1) is continuous.

3. The function (2.1) has compact support, i.e., vanishes identically outside of a bounded interval.

Deslauriers-Dubuc [79] constructed such functions (2.1) via a hierarchical interpolation process based on Lagrange polynomials. More precisely, they considered 2M support points for a positive integer M, and constructed the 2M Lagrange polynomials

$$ \mathcal{L}_{k,M}(x) = \prod_{j=-M+1, j \neq k}^{M} (x-j) / \prod_{j=-M+1, j \neq k}^{M} (k-j) \qquad (2.3) $$

associated to the points

$$ k \in \{-M+1, -M+2, \dots, M-1, M\} \qquad (2.4) $$

on the real line $\mathbb{R}$. These polynomials have order 2M and satisfy the interpolation conditions

$$ \mathcal{L}_{k,M}(j) = \delta_{jk}, \\ j, k \in \{-M+1, -M+2, \dots, M-1, M\}. \qquad (2.5) $$
<figure><img src="image_2.png" /><figcaption>Fig. 1 Scaling functions for DD wavelets of low order M</figcaption></figure>

Then, the function $\phi$ is constructed recursively on dyadic numbers via the rule

$$ \phi\left(\frac{n + \frac{1}{2}}{2^l}\right) = \sum_{k=-M+1}^{M} \mathcal{L}_{k,M}\left(\frac{1}{2}\right) \phi\left(\frac{n + k}{2^l}\right), \quad l = 0, 1, \dots, n \in \mathbb{Z}, \qquad (2.6) $$

starting from the condition (2.2). To be more specific, on the first level, the construction (2.6) means the following

$$ \phi\left(n + \frac{1}{2}\right) = \frac{1}{2}\phi(n) + \frac{1}{2}\phi(n + 1) \qquad (2.7) $$

for $M=1$, i.e., two support points.

In case of $M=2$, i.e., four support points, we obtain

$$ \phi\left(n + \frac{1}{2}\right) = -\frac{1}{16}\phi(n - 1) + \frac{9}{16}\phi(n) + \frac{9}{16}\phi(n + 1) - \frac{1}{16}\phi(n + 2). \qquad (2.8) $$

In turn, for $M=3$, i.e., six support points, it follows

<figure><img src="image_3.png" /><figcaption>Fig. 2 Zoom at the different scaling functions for DD wavelets shown in Fig. 1</figcaption></figure>

$$ \phi\left(n + \frac{1}{2}\right) = \frac{3}{256}\phi(n - 2) - \frac{25}{256}\phi(n - 1) + \frac{75}{128}\phi(n) + \frac{75}{128}\phi(n + 1) - \frac{25}{256}\phi(n + 2) + \frac{3}{256}\phi(n + 3). \qquad (2.9) $$
The symmetry of the interpolation process w.r.t. the point at which the interpolated value is sought turns out to be apparent from these formulas.

The constructed function (2.1) is uniformly continuous, and may thus be extended to a continuous function on the whole real line. Actually, for increasing number of support points $2M$, the regularity of the constructed function is increasing. Indeed, as observed by Saito-Beylkin [82], the representation

$$ \phi(x) = \int_{-\infty}^{\infty} \Phi(y) \Phi(y-x) dy, \quad x \in \mathbb{R}, \qquad (2.10) $$

holds with the scaling function $\Phi$ of the Daubechies wavelets of order $M$. Put differently, the scaling function $\phi$ of the Deslaurier-Dubuc wavelets arises as the autocorrelation function associated to the scaling function $\Phi$ of the corresponding Daubechies wavelets. In particular, the function (2.1) inherits the smoothness properties of the scaling function $\Phi$ characteristic for the Daubechies wavelets. The Deslauriers-Dubuc fundamental functions for low orders $M \in \{1, 2, 3, 4\}$ are shown in Figs. 1 and 2. We observe that these functions do in fact satisfy the interpolation condition (2.2), and show increasing smoothness for increasing order. For the lowest order, $M = 1$, corresponding to two support points, we recover the well-known hat function traditionally used as the shape function for linear finite elements in a single spatial dimension. We also observe the compact support of the interpolating wavelets.

By construction (2.6), continuity and taking into account the definition (2.2), the function (2.1) satisfies the refinement condition

$$ \phi(x) = \sum_{k \in \mathbb{Z}} h_k \phi(2x - k), \quad x \in \mathbb{R}, \qquad (2.11) $$

with the filter coefficients

$$ h_k = \phi\left(\frac{k}{2}\right), \quad k \in \mathbb{Z}, \qquad (2.12) $$

which attain the values [83]

$$ h_{2k} = \delta_{0,k}, \\ h_{2k+1} = \begin{cases} \mathcal{L}_{k,M}(\frac{1}{2}), & k = -M+1, -M+1, \dots, M, \\ 0, & \text{otherwise,} \end{cases} \qquad (2.13) $$

for $k \in \mathbb{Z}$. The filter coefficients for the cases (2.7)–(2.9) are recorded in Tab. 1. These coefficients were used to construct Fig. 1 showing the associated fundamental functions (2.1).

The constructed function (2.1) may be rescaled

$$ \phi_k^\ell : \mathbb{R} \to \mathbb{R}, \quad \phi_k^\ell(x) = \phi\left(2^\ell x - k\right), $$

<table><caption>**Table 1** Potentially non-zero filter coefficients $h_k$ of Deslauriers—Dubuc wavelets for $2M$ support points [70, 83], taking into account the symmetry $h_{-k} = h_k, k \in \mathbb{Z}$</caption><thead><tr><th></th><th>$h_0$</th><th>$h_1$</th><th>$h_2$</th><th>$h_3$</th><th>$h_4$</th><th>$h_5$</th><th>$h_6$</th><th>$h_7$</th></tr></thead><tbody><tr><td>$M = 1$</td><td>1</td><td>$\frac{1}{2}$</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>$M = 2$</td><td>1</td><td>$\frac{9}{16}$</td><td>0</td><td>$-\frac{1}{16}$</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>$M = 3$</td><td>1</td><td>$\frac{75}{128}$</td><td>0</td><td>$-\frac{25}{256}$</td><td>0</td><td>$\frac{3}{256}$</td><td>0</td><td>0</td></tr><tr><td>$M = 4$</td><td>1</td><td>$\frac{1225}{2048}$</td><td>0</td><td>$-\frac{245}{2048}$</td><td>0</td><td>$\frac{49}{2048}$</td><td>0</td><td>$-\frac{5}{2048}$</td></tr></tbody></table>

$$ \ell = 0, 1, \dots, \quad k \in \mathbb{Z}, \qquad (2.14) $$

to set up a number of approximation spaces

$$ V^\ell = \left\{ f : \mathbb{R} \to \mathbb{R} \mid f = \sum_{k \in \mathbb{Z}} s_k^\ell \phi_k^\ell, \ s_k^\ell \in \mathbb{R}, \ k \in \mathbb{Z} \right\}, \\ \ell = 0, 1, \dots, \qquad (2.15) $$

where only a finite number of coefficients $s_k^\ell$ is allowed to be non-zero. Due to the self-similar interpolation (2.6) used to construct the Deslauriers-Dubuc fundamental function (2.1), the approximation spaces (2.15) turn out to be nested, i.e., the subspace relation

$$ V^{\ell+1} \subseteq V^{\ell}, \quad \ell = 0, 1, \dots, \qquad (2.16) $$

holds. Indeed, using the expression (2.12) of the filter coefficients in the fundamental equation (2.11) satisfied by the fundamental function (2.1), we may re-write the scaling function (2.14) in the form

$$ \phi_m^\ell(x) = \phi\left(2^\ell x - m\right) = \sum_{k \in \mathbb{Z}} h_k \phi\left(2\left(2^\ell x - m\right) - k\right) \\ = \sum_{k \in \mathbb{Z}} h_k \phi\left(2^{\ell+1} x - 2m - k\right), \quad x \in \mathbb{R}, \qquad (2.17) $$

i.e., we obtain the equation

$$ \phi_m^\ell = \sum_{k \in \mathbb{Z}} h_k \phi_{2m+k}^{\ell+1}. \qquad (2.18) $$

In particular, the nesting property (2.16) becomes clear.

Any continuous function $f \in C_0(\mathbb{R})$ with compact support may be interpolated

$$ \mathcal{P}^\ell : C_0(\mathbb{R}) \to V^\ell, \quad \mathcal{P}^\ell[f] \\ = \sum_{k \in \mathbb{Z}} s_k^\ell \phi_k^\ell \quad \text{with} \\ s_k^\ell = f\left(\frac{k}{2^\ell}\right), \quad \ell = 0, 1, \dots, \quad k \in \mathbb{Z}, \qquad (2.19) $$
to produce an element of the approximation space (2.15). Apparently, the operator (2.19) is a linear operator which moreover projects onto the space (2.15). The latter property follows directly from the relation

$$ \phi_k^\ell \left( \frac{j}{2^\ell} \right) = \phi \left( 2^\ell \left( \frac{j}{2^\ell} \right) - k \right) = \phi(j-k) = \delta_{jk} \quad (2.20) $$

implied by the definition (2.14) and the property (2.2).

However, the projector (2.19) is not orthogonal in any reasonable sense, which is not surprising because no inner product was used in its construction. Rather, the operator permits to refine the inclusion relation into a (non-orthogonal) direct sum

$$ V^{\ell+1} = V^{\ell} + W^{\ell} \quad (2.21) $$

with the wavelet or detail space

$$ W^{\ell} = \left\{ f \in V^{\ell+1} \mid \mathcal{P}^{\ell}[f] = 0 \right\}, \quad (2.22) $$

which contains all functions that vanish on the coarse grid at level $\ell$. In the interpolating setting, a basis $\{\psi_k^\ell\}_{k \in \mathbb{Z}}$ of the space (2.22) may be readily constructed in the form

$$ \psi_k^\ell(x) = \psi(2^\ell x - k), \quad x \in \mathbb{R}, \ k \in \mathbb{Z}, \quad (2.23) $$

with the mother wavelet

$$ \psi : \mathbb{R} \to \mathbb{R}, \quad \psi(x) = \phi(2x - 1). \quad (2.24) $$

Due to the decomposition (2.21), an arbitrary element $f$ of the space $W^{\ell+1}$ may be written in the form

$$ f = \sum_{k \in \mathbb{Z}} s_k^\ell \phi_k^\ell + \sum_{k \in \mathbb{Z}} d_k^\ell \psi_k^\ell \quad (2.25) $$

with the coarse coefficients $s_k^\ell \in \mathbb{R}$ and the detail coefficients $d_k^\ell \in \mathbb{R}$ in addition to the "standard" representation

$$ f = \sum_{k \in \mathbb{Z}} s_k^{\ell+1} \phi_k^{\ell+1}, \quad (2.26) $$

implied by the definition (2.15). The refinement relation (2.18) together with the identity

$$ \begin{aligned} \psi_k^\ell(x) &= \psi(2^\ell x - k) = \phi(2(2^\ell x - k) - 1) \\ &= \phi(2^{\ell+1} x - 2k - 1) = \phi_{2k+1}^{\ell+1} \end{aligned} \quad (2.27) $$

permits to translate the wavelet representation (2.25) into the scaling-function representation (2.26), called wavelet synthesis [73]. The inverse operation is also possible and is referred to as wavelet analysis [73].

Applying the decomposition (2.25) repeatedly produces a multiscale decomposition of an element $f \in V^{\ell+1}$. By a limit argument, a corresponding multiscale composition applies to an arbitrary continuous function $f \in C_0(\mathbb{R})$ with compact support.

For the work at hand, also derivatives of the investigated functions are required. For a given function

$$ f = \sum_{k \in \mathbb{Z}} s_k^\ell \phi_k^\ell \in V^\ell \quad (2.28) $$

represented in the wavelet basis (2.14), we consider the interpolated derivative

$$ \tilde{\partial} f := \mathcal{P}^{\ell}[f'] \quad (2.29) $$

After some algebra, see Appendix B, we obtain the representation

$$ \mathcal{P}^{\ell}[f'] = \sum_{k \in \mathbb{Z}} \tilde{s}_k^{\ell} \phi_k^{\ell} \in V^{\ell} \quad \text{with} \quad \tilde{s}_k^{\ell} = 2^{\ell} \sum_{j \in \mathbb{Z}} s_{k+j}^{\ell} r_j \quad (2.30) $$

and suitable coefficients $r_k$.

For the considered Deslauriers-Dubuc wavelets of low order, we list the relevant coefficients in Tab. 2.

We close this section by an elementary remark: Up to a prefactor, the finite-difference formulae (B.7) lead to conventional finite-difference formulas for the low orders $M = 1$ and $M = 2$. More precisely, for $M = 1$, we are led – up to a prefactor – to the standard second-order central difference scheme (CDS)

$$ \frac{\tilde{s}_k^\ell}{2^\ell} = \frac{s_{k+1}^\ell - s_{k-1}^\ell}{2}, \quad (2.31) $$

whereas the derivative coefficients for $M = 2$ correspond to the coefficients of the $4^{th}$ order central difference scheme [84]

$$ \frac{\tilde{s}_k^\ell}{2^\ell} = \frac{s_{k+2}^\ell - 8s_{k+1}^\ell + 8s_{k-1}^\ell - s_{k+2}^\ell}{12} \quad (2.32) $$

However, the coefficients for the $6^{th}$ and $8^{th}$-order CDS schemes, for instance given by Fornberg [85] do not match those of the Deslauriers-Dubuc wavelets with three and four compact support points, as provided in Tab. 2. We refer to McCormick-Wells [86] for a more detailed discussion.

## 2.2 Discretizing the cell problem via wavelet collocation

We collected the relevant information on periodic inelastic homogenization of small-strain quasi-static problems in Appendix A. We suppose that the unit cell (A.1) is discretized
<table><caption>**Table 2** Non-zero derivative coefficients (2.30) for positive *k* and Deslauriers-Dubuc wavelets with 2*M* support points. Coefficients for negative *k*-values are obtained from the rule $r_{-k} = -r_k$</caption><thead><tr><th></th><th>r<sub>1</sub></th><th>r<sub>2</sub></th><th>r<sub>3</sub></th><th>r<sub>4</sub></th><th>r<sub>5</sub></th><th>r<sub>6</sub></th></tr></thead><tbody><tr><td>M = 1</td><td><sup>1</sup>&frasl;<sub>2</sub></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>M = 2</td><td>-<sup>2</sup>&frasl;<sub>3</sub></td><td><sup>1</sup>&frasl;<sub>12</sub></td><td></td><td></td><td></td><td></td></tr><tr><td>M = 3</td><td>-<sup>272</sup>&frasl;<sub>365</sub></td><td><sup>53</sup>&frasl;<sub>365</sub></td><td>-<sup>16</sup>&frasl;<sub>1095</sub></td><td>-<sup>1</sup>&frasl;<sub>2920</sub></td><td></td><td></td></tr><tr><td>M = 4</td><td>-<sup>39296</sup>&frasl;<sub>49553</sub></td><td><sup>76113</sup>&frasl;<sub>396424</sub></td><td>-<sup>1664</sup>&frasl;<sub>49553</sub></td><td><sup>2645</sup>&frasl;<sub>1189272</sub></td><td><sup>128</sup>&frasl;<sub>743295</sub></td><td>-<sup>1</sup>&frasl;<sub>1189272</sub></td></tr></tbody></table>

by $N_a$ voxels in each coordinate direction $a \in \{1, 2, \dots, d\}$, and assemble these numbers into an integer-valued vector

$$N = (N_1, N_2, \dots, N_d). \qquad (2.33)$$

The discretized regular grid reads

$$Y_N = \left\{ \mathbf{x} \in Y \mid \begin{array}{l} x_a = \frac{j_a L_a}{N_a} \quad \text{for some } \mathbf{j} \in \mathbb{Z}_N^+ \\ \text{and all } a = 1, 2, \dots, d \end{array} \right\} \qquad (2.34)$$

with the index set

$$\mathbb{Z}_N^+ = \left\{ \mathbf{j} \in \mathbb{Z}^d \mid 0 \le j_a < N_a, \ a = 1, 2, \dots, d \right\} \qquad (2.35)$$

and where we suppress the dependence on the unit-cell dimensions

$$L = (L_1, L_2, \dots, L_d) \qquad (2.36)$$

throughout the manuscript for the sake of readability. We define the ansatz functions

$$\phi_m^N(\mathbf{x}) = \prod_{a=1}^{d} \phi\left(\frac{x_a}{L_a} N_a - m_a \mod N_a\right) \quad \text{for } \mathbf{x} \in Y \\ \text{and } \mathbf{m} \in \mathbb{Z}_N^+ \qquad (2.37)$$

by multiplying rescaled and shifted versions of the fundamental function (2.1) characterizing the Deslauriers-Dubuc wavelets. By construction, i.e., via the use of the modulo operator $\mod$, the rescaled functions (2.37) are periodic. For the lowest order $M = 1$, the fundamental function (2.1) of the Deslauriers-Dubuc wavelets is just a hat function. Therefore, the definition (2.37) leads to the conventional multilinear ansatz functions which are standard for the finite-element method (FEM) on hexahedral grids, i.e., bilinear ansatz functions in two and trilinear ansatz functions in three dimensions. For higher-DD order $M$, the fundamental function (2.1) is even smoother. In particular, the ansatz functions (2.37) define elements of the first-order Sobolev space

$$\phi_m^N \in H_{\text{per}}^1(Y), \quad \mathbf{m} \in \mathbb{Z}_N^+, \qquad (2.38)$$

of periodic functions. Notice that the functions (2.37) do not have vanishing integral mean, i.e., are not contained in the space $H_{\#}^1(Y)$.

With these considerations at hand, it is natural to define the wavelet ansatz space

$$H_N(Y)^d = \left\{ \mathbf{u} = \sum_{m \in \mathbb{Z}_N^+} \phi_m^N \mathbf{u}_m \mid \mathbf{u}_m \in \mathbb{R}^d, \ \mathbf{m} \in \mathbb{Z}_N^+, \right. \\ \left. \text{s.t. } \sum_{m \in \mathbb{Z}_N^+} \mathbf{u}_m = \mathbf{0} \right\} \qquad (2.39)$$

for the displacements and the space

$$L_N^2(Y; \text{Sym}(d)) = \left\{ \tau = \sum_{m \in \mathbb{Z}_N^+} \phi_m^N \tau_m \mid \tau_m \in \text{Sym}(d), \right. \\ \left. \mathbf{m} \in \mathbb{Z}_N^+ \right\} \qquad (2.40)$$

containing discrete strain and stress fields. The discrete symmetrized gradient operator

$$\nabla_N^s : H_N(Y)^d \to L_N^2(Y; \text{Sym}(d)) \qquad (2.41)$$

is realized by interpolating the symmetrized gradient at the points of the grid (2.34), i.e., the explicit representation

$$\nabla_N^s \mathbf{u} = \sum_{m \in \mathbb{Z}_N^+} \phi_m^N \nabla^s \mathbf{u}(\mathbf{x}_m) \\ \text{for all fields } \mathbf{u} \in H_N(Y)^d, \qquad (2.42)$$

emerges at least for DD-order $M$ exceeding unity. By the observations of the previous section on the finite-difference form (B.7) of the interpolated derivative (B.3), we arrive at the explicit component-wise form

$$[\nabla_N^s \mathbf{u}(\mathbf{x}_j)]_{ab} = \sum_{k \in \mathbb{Z}} \frac{r_k}{2} \left( \frac{1}{h_a} \mathbf{e}_b \cdot \mathbf{u}_{j+k} \mathbf{e}_a + \frac{1}{h_b} \mathbf{e}_a \cdot \mathbf{u}_{j+k} \mathbf{e}_b \right), \\ a, b = 1, 2, \dots, d, \qquad (2.43)$$
for a discrete displacement field $\mathbf{u} \in H_N^1(Y)^d$ in the specified form (2.39) and with the mesh spacing

$$h_a = \frac{L_a}{N_a}, \quad a = 1, 2, \dots, d. \qquad (2.44)$$

By the formula (2.43) and the prescription (2.31), we may define the symmetrized gradient operator (2.41) for the DD-order $M = 1$, as well. In this case, the ansatz space (2.37) coincides with $Q_1$-finite elements on a regular grid. However, the wavelet-type derivative operator (2.43) is different. Indeed, the derivative (2.42) is evaluated at the nodes (2.34), where the ansatz functions are not differentiable. For classical $Q_1$ elements with full or reduced integration [36, 87], the derivative is evaluated inside the element.

Similar to the symmetrized gradient, a divergence-type operator may be defined

$$\mathrm{div}_N : L_N^2(Y; \mathrm{Sym}(d)) \to H'_N(Y)^d \equiv H_N(Y)^d, \quad (2.45)$$

which reads explicitly

$$\mathrm{div}_N \tau(\mathbf{x}_j) = \sum_{k \in \mathbb{Z}} \sum_{a=1}^{d} \frac{r_k}{h_a} \tau_{j+k} \mathbf{e}_a \cdot \mathbf{e}_a, \\ \tau \in L_N^2(Y; \mathrm{Sym}(d)). \qquad (2.46)$$

With these notations at hand, and for prescribed macroscopic strain $\bar{\mathbf{e}}$, the DD-wavelet discretization on the regular grid (2.34) seeks fields

$$\mathbf{u}_N \in H_N(Y)^d \quad \text{and} \quad \mathbf{e}_N, \sigma_N \in L_N^2(Y; \mathrm{Sym}(d)), \quad (2.47)$$

s.t. the balance of linear momentum

$$\mathrm{div}_N \sigma_N = \mathbf{0}, \qquad (2.48)$$

the constitute law

$$\sigma_N(\mathbf{x}_j) = \mathbf{S}(\mathbf{x}_j, \mathbf{e}_N(\mathbf{x}_j)), \quad \mathbf{j} \in \mathbb{Z}_N^+, \qquad (2.49)$$

and the kinematic compatibility condition

$$\mathbf{e}_N = \bar{\mathbf{e}} + \nabla_N^s \mathbf{u}_N \qquad (2.50)$$

are satisfied. The problem (2.48)-(2.50) is obtained by wavelet collocation, i.e., neither the balance equation (2.48) nor the compatibility condition (2.50) are satisfied outside of the collocation points. Moreover, the constitutive relation (2.49) only holds at the collocation points (2.34), as well. These characteristics are different than for FEM, where both compatibility (2.50) and the constitutive law (2.49) hold everywhere, whereas only the equilibrium equation (2.48) is

approximated. For the classical Moulinec-Suquet discretization [1, 2, 23], both differential relations (2.48) as well as (2.50) are valid on the entire cell, and the constitutive law holds in the voxel centers only, as in eq. (2.49). We refer to the work [55] for a modern discussion of the Moulinec-Suquet discretization.

At this point we wish to state explicitly the consequences of the findings (2.31)-(2.32): Low-order Deslauriers-Dubuc wavelets on a regular grid coincide with discretization schemes that were already considered in FFT-based computational micromechanics:

* Deslauriers-Dubuc wavelets of order $M = 1$ are equivalent to 2nd order central differences [39, 40] on a regular grid.

* Deslauriers-Dubuc wavelets of order $M = 2$, as considered by Kaiser et al. [70] are equivalent to 4th order central differences [41] on a regular grid.

However, Deslauriers-Dubuc wavelets of orders exceeding two differ from higher-order central difference schemes as, e.g., discussed by Vidyasagar et al. [42].

Once the problem (2.48)-(2.50) is solved, the approximate mean stress is obtained by averaging the stress values at the collocation points

$$\bar{\sigma}_N = \frac{1}{\mathrm{prod}(N)} \sum_{m \in \mathbb{Z}_N^+} \sigma_N(\mathbf{x}_m) \qquad (2.51)$$

with the product

$$\mathrm{prod}(N) = N_1 N_2 \cdots N_d. \qquad (2.52)$$

We will neither discuss the solvability properties of the problem nor the convergence under grid refinement. Rather, we turn our attention to the Lippmann-Schwinger type solvers proposed by Kaiser et al. [70]. Repeating the steps made in the continuous setting for a fixed reference material $\mathbb{C}^0$, see Appendix A.2, we are led to the Lippmann-Schwinger equation

$$\mathbf{e}_N = \bar{\mathbf{e}} - \Gamma_N^0 : \tau_N \qquad (2.53)$$

with the stress polarization $\tau_N \in L_N^2(Y; \mathrm{Sym}(d))$, defined by the interpolation condition

$$\tau_N(\mathbf{x}_j) = \mathbf{S}(\mathbf{x}_j, \mathbf{e}_N(\mathbf{x}_j)) - \mathbb{C}^0 : \mathbf{e}_N(\mathbf{x}_j), \quad \mathbf{j} \in \mathbb{Z}_N^+(2.54)$$

and the wavelet-type Eshelby-Green operator

$$\Gamma_N^0 = \nabla_N^s \mathbf{G}_N^0 \mathrm{div}_N. \qquad (2.55)$$
It remains to specify Green's operator

$$G_N^0 \in L(H'_N(Y)^d, H_N(Y)^d) \qquad (2.56)$$

whose action $u_N = G_N^0 f_N$ is implicitly defined by the solution $u_N \in H_N(Y)^d$ of the reference problem

$$\mathrm{div}_N C^0 : \nabla_N^s u_N = f_N \qquad (2.57)$$

for a given element $f_N \in H'_N(Y)^d$. Due to the similarity to finite-difference discretizations [27, 33, 36], we follow Kaiser et al. [70], and consider the discrete Fourier transform (DFT)

$$\hat{g}[\xi] = \sum_{j \in \mathbb{Z}_N^+} g_j e^{-2\pi i \frac{\xi_j}{N}}, \quad \xi \in \mathbb{Z}_N^+, \qquad (2.58)$$

of a given (possibly vector- or tensor-valued) array

$$\{g_j\}_{j \in \mathbb{Z}_N^+} \qquad (2.59)$$

with the short-hand notation

$$\frac{\xi \cdot j}{N} = \sum_{a=1}^{d} \frac{\xi_a j_a}{N_a}. \qquad (2.60)$$

The transform (2.58) may be readily inverted via the operation

$$g_j = \frac{1}{\mathrm{prod}(N)} \sum_{\xi \in \mathbb{Z}_N^+} \hat{g}[\xi] e^{2\pi i \frac{\xi_j}{N}}, \quad j \in \mathbb{Z}_N^+, \qquad (2.61)$$

the inverse DFT. Due to the finite-difference structure (2.43), the Fourier coefficients of the symmetrized gradient (2.41) may be expressed in the form

$$(\widehat{\nabla}_N^s \mathbf{u})[\xi] = \hat{\mathbf{r}}[\xi] \otimes_s \hat{u}_N[\xi], \quad \xi \in \mathbb{Z}_N^+, \qquad (2.62)$$

involving the vector-valued array $\{\hat{\mathbf{r}}[\xi]\}_{\xi \in \mathbb{Z}_N^+}$ with components

$$\hat{r}_a[\xi] = \sum_{k \in \mathbb{Z}} r_k e^{2\pi i \frac{k \xi_a}{N_a}}, \quad a = 1, 2, \dots, d. \qquad (2.63)$$

As only a finite number of the values $\{r_k\}_{k \in \mathbb{Z}}$ are non-zero, the latter sum (2.63) is actually finite. By a similar reasoning, the discrete divergence operator (2.46) admits the representation

$$(\mathrm{div}_N \boldsymbol{\tau})[\xi] = \hat{\boldsymbol{\tau}}[\xi] \cdot \hat{\boldsymbol{r}}[\xi], \quad \xi \in \mathbb{Z}_N^+. \qquad (2.64)$$

For later reference, we remark that the skew symmetry $r_{-k} = -r_k, k \in \mathbb{Z}$, of the filter coefficients (B.8) implies the alternative representation

$$\hat{r}_a[\xi] = 2i \sum_{k=1}^{\infty} r_k \sin \left( 2\pi \frac{k \xi_a}{N_a} \right), \quad a = 1, 2, \dots, d, \qquad (2.65)$$

of the vector components as a consequence of Euler's identity

$$e^{i\alpha} = \cos \alpha + i \sin \alpha, \quad \alpha \in \mathbb{R}. \qquad (2.66)$$

In particular, the components (2.63) are purely imaginary and satisfy the identity

$$\overline{\hat{r}_a[\xi]} = -\hat{r}_a[\xi], \quad a = 1, 2, \dots, d. \qquad (2.67)$$

With the explicit representations (2.62) and (2.64) at hand, the reference problem (2.57) may be re-written in Fourier space in the form

$$\hat{\mathbf{r}}[\xi] \cdot C^0 : (\hat{\mathbf{r}}[\xi] \otimes_s \hat{u}_N[\xi]) = \hat{f}_N[\xi], \quad \xi \in \mathbb{Z}_N^+. \qquad (2.68)$$

For an isotropic reference material $C^0$ with Lamé constants $\lambda^0$ and $\mu^0$, we may express the formula (2.68) more explicitly

$$(\mu^0 + \lambda^0) (\overline{\hat{\mathbf{r}}[\xi]} \cdot \hat{u}_N[\xi]) \hat{\mathbf{r}}[\xi] + \mu^0 \|\hat{\mathbf{r}}[\xi]\|^2 \hat{u}_N[\xi] = -\hat{f}_N[\xi], \quad \xi \in \mathbb{Z}_N^+. \qquad (2.69)$$

In particular, we may invert eq. (2.68) directly:

$$\hat{u}_N[\xi] = - \frac{1}{\mu^0 \|\hat{\mathbf{r}}[\xi]\|^2} [\mathbf{I} + \frac{\mu^0 + \lambda^0}{2\mu^0 + \lambda^0} \frac{\hat{\mathbf{r}}[\xi]}{\|\hat{\mathbf{r}}[\xi]\|} \otimes \frac{\overline{\hat{\mathbf{r}}[\xi]}}{\|\overline{\hat{\mathbf{r}}[\xi]}\|}] \hat{f}_N[\xi], \quad \xi \in \mathbb{Z}_N^+, \qquad (2.70)$$

for all frequencies $\xi$ in the set

$$A_N = \left\{ \xi \in \mathbb{Z}_N^+ \mid \hat{\mathbf{r}}[\xi] \neq \mathbf{0} \right\} \qquad (2.71)$$

of admissible frequencies. More explicitly, there are two cases where the component (2.63) vanishes, which are evident from the representation (2.65). For trivial reasons, the identity

$$\hat{r}_a[\xi] \Big|_{\xi_a=0} = 2i \sum_{k=1}^{\infty} r_k \underbrace{\sin(0)}_{=0} = 0, \quad a = 1, 2, \dots, d, \qquad (2.72)$$
holds. Moreover, for the middle frequency, we observe

$$ \hat{r}_a[\xi] \Big|_{\xi_a = \frac{N_a}{2}} = 2i \sum_{k=1}^{\infty} r_k \underbrace{\sin(k\pi)}_{=0} = 0, \qquad (2.73) \\ a = 1, 2, \dots, d. $$

Apparently, this scenario can only occur if the voxel count $N_a$ is even.

These are the only cases where the vector $\hat{r}[\xi]$ vanishes, i.e., we find between one and $2^d$ exceptional cases in $d$ spatial dimensions.

With the observation (2.70) at hand, we finally find the explicit representation [70]

$$ \begin{aligned} \widehat{\Gamma}_N^0 : \tau_N[\xi] &= \left\{ \begin{array}{ll} \frac{1}{\mu^0} \frac{\hat{r}}{\|\hat{r}\|} \otimes_s (\hat{\tau}_N[\xi] \frac{\hat{r}}{\|\hat{r}\|}) \frac{\mu^0 + \lambda^0}{\mu^0 (2\mu^0 + \lambda^0)} & \xi \in \mathbb{A}_N, \\ \left( \frac{\hat{r}}{\|\hat{r}\|} \cdot \hat{\tau}_N[\xi] \frac{\hat{r}}{\|\hat{r}\|} \right) \frac{\hat{r}}{\|\hat{r}\|} \otimes \frac{\hat{r}}{\|\hat{r}\|}, & \text{otherwise,} \end{array} \right. \end{aligned} \quad (2.74) $$

of the discrete Eshelby-Green operator (2.55) in Fourier space. In eq. 2.74, we suppressed the dependence on the frequency vector $\xi$ of the vector $\hat{r}$ in the right-hand side.

# 3 Computational experiments

## 3.1 Setup

The Deslauriers-Dubuc wavelet discretizations of order $M$ in the range {1, 2, 3, 4} described in section 2 were implemented into an in-house FFT-based computational micromechanics code. This solver is written in Python with Cython extensions for the time-critical parts, uses OpenMP for shared-memory parallelization and employs the FFTW [88] library for computing the fast Fourier transform (2.58). The computations presented in the following were run on six to twelve cores of a Intel Xenon Platinum 8480+ processor, a single node of the UDE supercomputer amplitUDE.

We compare the wavelet discretizations on regular grids with three discretizations which are presumably the most used nowadays:

(i) The original discretization of Moulinec and Suquet [1, 2] based on trigonometric polynomials and the trapezoid rule [89],

(ii) The discretization on a rotated staggered grid introduced by Willot [33], which is equivalent to $Q_1$-finite elements with reduced integration,

(iii) The staggered grid discretization proposed by Schneider et al. [30].

To solve the discretized systems of equations, we consider three solvers:

1. The basic scheme, i.e., the original solver developed by Moulinec-Suquet [1, 2], see Alg. 1,

2. The Barzilai-Borwein (BB) method [90], which may be interpreted as a basic scheme with adaptive reference material, see Alg. 2,

3. The conjugate gradient (CG) method [16, 91].

These solvers differ in the memory consumption and the convergence speed. CG is the fastest, but is typically implemented on four strain fields. The basic scheme, in contrast, may be implemented on a single strain field only, but is strongly inferior convergence behavior compared to CG. The BB method offers a compromise, i.e., may be implemented on two strain fields, but shows a convergence behavior which is only a little worse than CG.

The interested reader may find a summary of the implementation of these discretizations and solvers in the review articles [7–9]. For the sake of completeness, we use the termination criterion

$$ \sqrt{\alpha_0} \|\Gamma_N^0 : \sigma_N\|_{L^2} \leq tol \|\bar{\sigma}_N\| \quad (3.1) $$

to quantify how well the equilibrium equation (2.48) is satisfied and restrict to reference materials of the form (A.31) with general parameter $\alpha_0$, see Schneider [9, §3.6]. In this case, the formula (2.74) simplifies to

$$ \widehat{\Gamma}_N^0 : \tau_N[\xi] = \left\{ \begin{array}{ll} \frac{2}{\alpha_0} \frac{\hat{r}}{\|\hat{r}\|} \otimes_s (\hat{\tau}_N[\xi] \frac{\hat{r}}{\|\hat{r}\|}) + \frac{1}{\alpha_0} & \xi \in \mathbb{A}_N, \\ \left( \frac{\hat{r}}{\|\hat{r}\|} \cdot \hat{\tau}_N[\xi] \frac{\hat{r}}{\|\hat{r}\|} \right) \frac{\hat{r}}{\|\hat{r}\|} \otimes \frac{\hat{r}}{\|\hat{r}\|}, & \text{otherwise,} \end{array} \right. \quad (3.2) $$

Unless specified otherwise, we use the value $tol = 10^{-5}$ for our studies.

## 3.2 A single-sphere microstructure

We examine a three-dimensional microstructure comprising an aluminium 6061-T6 matrix phase with an embedded spherical silicon carbide inclusion. The centered inclusion occupies 30 % of the total volume, while the remaining space is filled by the matrix phase. Figure 3 shows the microstructure at different resolutions, with a resolution varying from 32 to 256 voxels per axis. The used material properties of aluminium 6061-T6 and silicon carbide (SiC), as reported by Segurado et al. [92], are summarized in table 3.

First, we perform a resolution analysis, for a 1 % uniaxial strain loading in the $x$-direction. Fig. 4a shows the corre-
<figure><img src="image_4.png" /><figcaption>**Fig. 3** 3D ball microstructure with four discretization levels</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 4** 3D ball microstructure subjected to a normal strain: macroscopic stress and solver convergence for various discretizations and solvers</figcaption></figure>

<table><caption>**Table 3** Material properties [92] for the single spherical inclusion shown in Fig. 3</caption><thead><tr><td>Material</td><td><i>E</i> in GPa</td><td><i>ν</i></td></tr></thead><tbody><tr><td>Aluminium 6061-T6</td><td>75</td><td>0.3</td></tr><tr><td>SiC</td><td>400</td><td>0.2</td></tr></tbody></table>

sponding macroscopic stress $\bar{\sigma}_{11}$ obtained with four different discretizations: trigonometric polynomials [1], the rotated staggered grid [33], the staggered grid [30], and the DD2 wavelet discretization [70]. We observe that the trigonometric polynomial and the DD2 wavelet discretizations on the one hand, and the rotated staggered grid and staggered grid discretizations on the other hand, lead to almost identical macroscopic stresses. The rotated staggered grid and classical staggered grid discretizations lead to stresses that are at most 0.02 % higher than the other two discretizations. From

Fig. 4a, it also appears that the macroscopic stress is, on average, reduced by 0.2 % when increasing the resolution from 32 to 128 voxels per axis, but only by 0.01 % when increasing this number from 128 to 256 voxels per axis. Hence, we adopt microstructure with $128^3$ voxels for the subsequent investigations.

Fig. 4b shows the evolution of the residual for three solvers: the basic scheme [1, 2], the Barzilai-Borwein (BB) solver [90], and the conjugate gradient (CG) solver [16, 91]. These solvers are combined with the four discretizations already investigated in Fig. 4a. For all discretizations, the basic scheme converges in 39 iterations, BB in 17 iterations, and CG in 13 iterations. For the present example, the required number of iterations depends on the chosen solver, but the influence of the discretization is rather limited.

For FFT-based homogenization solvers, the necessary number of iterations depends on the material contrast in the microstructure, which was rather low in the previous exam-
<figure><img src="image_6.png" /><figcaption>Fig. 5 Effect of the material contrast on the convergence for the studied discretizations</figcaption></figure>

ple. Fig. 5a reports on the number of required iterations as a function of the material contrast, for various discretizations combined with the CG solver. In this example, the elastic properties of the matrix are identical to those prescribed in Tab. 3. The Poisson ratio of the inclusion is set identical to the one of the matrix, i.e., $\nu = 0.3$. Here, we consider the "material contrast"

$$E_{\text{inclusion}}/E_{\text{matrix}}, \qquad (3.3)$$

which is varied between $10^{-5}$ and $10^5$ and permits us to distinguish soft and hard inclusions. Fig. 5a reveals that the DD2 wavelet discretization shows a similar behavior as the Moulinec-Suquet discretization. To be more precise, the sensitivity to the phase contrast appears to be rather strong for both discretizations: CG terminates after 1029 and 1148 iterations for a contrast of $10^{-5}$ as well as 2268 and 2398 iterations for a contrast of $10^5$. In contrast, both the staggered and the rotated staggered grid discretizations are less sensitive to the phase contrast, especially when approaching the porous case. In particular, for a contrast $E_{\text{inclusion}}/E_{\text{matrix}} = 10^{-5}$, the CG solver converges after 34 iterations for the staggered grid discretization and after 27 iterations for the rotated staggered grid discretization.

To investigate further, we furnish the spherical inclusion with a vanishing Young's modulus while the matrix material stays identical to the previous example. Fig. 5b shows the evolution of the residual of CG for different discretizations.

We observe that for both the Moulinec-Suquet and the DD2 wavelet discretizations the CG solver does not reach the tolerance value $10^{-5}$ within 10000 iterations. The prescribed tolerance is reached in 34 iterations for the staggered grid discretization, and in 27 iterations for the rotated staggered grid discretization.

Fig. 5b also reports on the performance of three other wavelet discretizations: Deslauriers-Dubuc with one, three and four support points, denoted by DD1, DD3 and DD4, see section 2. The convergence behavior of CG for the DD3 and DD4 wavelet discretizations is similar as with the DD2 wavelet: The solver is unable to reach the prescribed tolerance within 10000 iterations. From these results, we conclude that increasing the number of support points of the wavelet discretization does not improve the numerical robustness for porous microstructures. In contrast, reducing the size of the support seems to have a positive effect: For DD1, the residual reaches the tolerance within 33 iterations. This behavior is similar to what we observed for the staggered and the rotated staggered grid discretizations.

Tab. 4 records the runtimes of ten iterations of the basic scheme for different discretizations on a single thread. We use the microstructure with a single inclusion shown in Fig. 3d and the material parameters listed in Tab. 3. We observe that all discretizations lead to a similar runtime, and the differences in the runtime are actually caused by the application of the Eshelby-Green operator (2.74) in Fourier space. The discretization via trigonometric polynomials [1, 2] is
<table><caption>**Table 4** Runtimes of ten iterations of the basic scheme on a single thread for the spherical inclusion shown in Fig. 3 and the material parameters listed in Tab. 3</caption><thead><tr><th>Discretization</th><th>Total runtime in s</th><th>Runtime 0 in s</th></tr></thead><tbody><tr><td>Trigonometric polynomials</td><td>36.1</td><td>3.9</td></tr><tr><td>Staggered grid</td><td>40.0</td><td>7.8</td></tr><tr><td>Rotated staggered grid</td><td>40.1</td><td>7.7</td></tr><tr><td>Wavelets DD1</td><td>37.8</td><td>4.9</td></tr><tr><td>Wavelets DD2</td><td>37.4</td><td>5.0</td></tr><tr><td>Wavelets DD3</td><td>36.8</td><td>4.9</td></tr><tr><td>Wavelets DD4</td><td>37.7</td><td>4.9</td></tr></tbody></table>

<table><caption>**Table 5** Material properties for the sand microstructure [95]</caption><thead><tr><th>Material</th><th><i>E</i> in GPa</th><th><i>ν</i></th></tr></thead><tbody><tr><td>Sand</td><td>66.9</td><td>0.25</td></tr><tr><td>Binder</td><td>43.0</td><td>0.17</td></tr></tbody></table>

the fastest, and only slightly trailed by the different wavelet discretizations. The two staggered grid variants are slightly slower. To sum up, the runtimes per iteration of the wavelet discretizations on a regular grid are comparable to the original discretization proposed by Moulinec-Suquet [1, 2].

### 3.3 A sand grain microstructure

Next, we analyze the porous sand-core microstructure [93, 94] presented in Fig. 6, which comprises 50 % sand and 2 % inorganic binder. The microstructure was generated using the stacked Mechanical Contraction method, proposed by Donval et al. [95], and discretized by 140³ voxels of edge length 4 µm. This example is challenging due to the high porosity of the composite and to the complex geometry of the microstructure. The elastic properties of the constituent materials are given in Tab. 5, and correspond to the ones used in Donval et al. [95].

The microstructure is subjected to a 1 % extension in the z-direction.

Fig. 7 shows the convergence behavior of the basic scheme, the Barzilai-Borwein solver and the Conjugate Gradient method, combined with five different discretizations: trigonometric polynomials integrated with the trapezoid rule, the rotated staggered, the staggered grid and both the DD1 and the DD2 wavelet discretizations. Fig. 7a shows that the basic scheme converges within 10 000 iterations when combined with either the DD1 wavelet discretization or the staggered grid discretization. In that case, the convergence is almost twice as fast for DD1 (3883 iterations) than for DD2 (7549 iterations). After 10 000 iterations, the basic scheme reaches a residual of $7.5 \times 10^{-4}$ when combined with the rotated staggered grid discretization, a residual of $2.8 \times 10^{-3}$ for the considered trigonometric polynomial discretization, and of $3.4 \times 10^{-3}$ for the DD2 wavelet discretization. We

<table><caption>**Table 6** Average macroscopic normal stress $\bar{\sigma}_{zz}$ in z-direction</caption><thead><tr><th>Discretization</th><th>$\bar{\sigma}_{zz}$ in MPa</th></tr></thead><tbody><tr><td>Staggered grid</td><td>49.47</td></tr><tr><td>Rotated staggered grid</td><td>49.56</td></tr><tr><td>Wavelet DD1</td><td>50.94</td></tr></tbody></table>

make similar observations for the BB solver, see Fig. 7b. Convergence is only achieved with the staggered grid and the DD1 wavelet discretizations. The number of iterations to convergence is drastically reduced compared to the basic scheme, however: 192 for the DD1 wavelet discretization and 340 for the staggered grid discretization. Finally, in case the CG solver is used, see Fig. 7d, convergence is reached for three discretizations: the DD1 wavelet, in 172 iterations, the staggered grid discretization, in 210 iterations, and the rotated staggered discretization, in 6378 iterations. After 10 000 iterations, the residual reaches $2.7 \times 10^{-3}$ for the DD2 wavelet discretization, and $3.6 \times 10^{-3}$ for the trigonometric polynomials discretization.

Tab. 6 reports the macroscopic stress component $\bar{\sigma}_{zz}$ for the three discretizations that reach convergence with CG: the staggered and rotated staggered grid discretizations, as well as the DD1 wavelet discretization.

The macroscopic stress obtained with the staggered and rotated staggered grid are only 0.2 % apart. The macroscopic stress obtained using the DD1 wavelet discretization, is 2.7 % higher.

Turning our attention to the local fields, shown in Fig. 8, we compare the local Mises equivalent stress field for all three converging discretizations. The local Mises equivalent stress field obtained using the staggered grid discretization is smooth, even around the binder bridges. For the rotated staggered discretization, we observe a slight checkerboard pattern. This checkerboard pattern is more pronounced for the DD1 wavelet discretization, especially in the vicinity of binder bridges. As a conclusion, for complex and porous microstructures, the DD1 wavelet discretization is as stable as the staggered grid discretization, but leads to a slightly
<figure><img src="image_7.png" /><figcaption>(a) 2D slice</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(b) 3D representation</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(a) Basic scheme</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(b) BB solver with zoom on the first 400 iterations</figcaption></figure>

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

<figure><img src="image_12.png" /><figcaption>(c) CG solver with zoom on the first 250 iterations</figcaption></figure>

<caption>**Fig. 7** Residual vs iterations for the 256³ 3D sand microstructure subjected to a normal strain: convergence behavior for five different discretizations in combination with three solvers</caption>
<figure><img src="image_13.png" /><figcaption>Fig. 8 Local Mises equivalent stress distribution on a slice of sand microstructure using CG solver</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 9 Considered short-fiber microstructure</figcaption></figure>

different value for the macroscopic stress computation, and produces local fields with more artifacts.

## 3.4 A fiber microstructure

For the last example, we investigate the capabilities of the wavelet discretizations in case of material non-linearities.

We consider a short-fiber reinforced polymer, presented in Fig. 9. The microstructure is composed of 10% glass fibers, with a length of 250 µm and a diameter of 10 µm, embedded in a nylon polymer matrix. The fibers were packed with isotropic orientation using the Sequential Addition and Migration approach [96]. The considered volume element, of edge length 512 µm, is discretized on a uniform grid of 256³ voxels. The fibers are considered elastic, while the matrix

follow a $J_2$-elastoplastic model with a linear-exponential hardening law [97]:

$$ \sigma_0(p) = \sigma_y + h_1 p + h_2(1 - \exp(-np)). \quad (3.4) $$

We adopt the material properties of the fiber and matrix phases from Wu et al. [97], as summarized in table 7.

We study the elastoplastic behavior of the microstructure under mixed boundary conditions [98]. A uniaxial strain of 5% is applied in the x-direction, distributed over 50 equidistant load steps, while the macroscopic stress in the transverse plane is set to zero. An affine extrapolation is used for the initial guess $\epsilon^0(t_{n+1})$ of the subsequent load step $t_{n+1}$, as proposed by Moulinec and Suquet [2]:
<table><caption>**Table 7** Material properties for the fiber microstructure [97]</caption><tr><td>Nylon matrix</td><td><i>E</i> = 3 GPa</td><td><i>ν</i> = 0.35</td><td><i>σ</i><sub>Y</sub> = 20 MPa<br><i>h</i><sub>1</sub> = 1 MPa</td><td><i>n</i> = 150<br><i>h</i><sub>2</sub> = 15 MPa</td></tr><tr><td>Glass fibers</td><td><i>E</i> = 72 GPa</td><td><i>ν</i> = 0.2</td><td></td><td></td></tr></table>

$$ \boldsymbol{\varepsilon}^0(t_{n+1}) = \boldsymbol{\varepsilon}^0(t_n) + \frac{t_{n+1} - t_n}{t_n - t_{n-1}} (\boldsymbol{\varepsilon}(t_n) - \boldsymbol{\varepsilon}(t_{n-1})), \quad (3.5) $$

where $\boldsymbol{\varepsilon}(t_n)$ and $\boldsymbol{\varepsilon}(t_{n-1})$ stand for the converged solutions of the two previous steps.

We examine the performance of several discretizations: trigonometric polynomials, rotated staggered grid, staggered grid and DD wavelets with one, two, three, and four support points.

We solve each of these discretizations by the non-linear conjugate gradient method [99], see Alg. 3. Fig. 10a shows the macroscopic strain-stress behavior of the composite. In the elastic phase, that is, until around $\bar{\varepsilon}_{xx} = 0.8\%$, all discretizations lead to the same stress. The difference in macroscopic stress starts to be visible for strains greater than $\bar{\varepsilon}_{xx} = 3\%$. From that strain onwards, the computed macroscopic stress is higher for all wavelet discretizations than for the other three discretizations (trigonometric polynomials, staggered and rotated staggered grid discretizations). The gap is highest for the highest considered strain value, $\bar{\varepsilon}_{xx} = 5\%$. For this value, the macroscopic strain lies between 52.7 MPa and 53.1 MPa for the wavelet discretizations, and between 51.7 MPa and 51.8 MPa for the other three discretizations. The maximum difference between the value obtained with the staggered grid discretization and with the DD2 wavelet discretization is 2.9%.

Figure 10b provides the number of iterations per load step. The staggered grid discretization always leads to the highest iteration count, while the DD1 wavelet discretization leads to the lowest iteration count up to the 33rd load step. After this load step, the DD4 wavelet discretization comes with the least iterations.

To conclude the present example, all the investigated wavelet discretizations are computationally efficient for elastoplastic problems, but slightly overestimate the macroscopic stresses.

## 4 Conclusions and perspectives

The work at hand was concerned with the discretization by wavelet collocation, introduced recently by Kaiser and co-workers [70], for the cell problem of periodic homogenization in linear elasticity. Due to the variety of established options for discretizing the unit-cell problem available for FFT-based methods, our goal was to put this "alien" discretization into perspective, comparing its construction and its properties to available discretization schemes.

In a sense, wavelet analysis has a different paradigm than traditional finite-element discretizations: Wavelet methods aim to construct a full multiscale<fnref n="1" /> basis of the entire function space the solution of the PDE under consideration is sought in. Thus, instead of fixing the discretization a priori, the continuous solution is actually the starting point, and adaptivity is naturally brought in to approximate the typically unknown discontinuities of the solution appropriately. Due to the full multiscale nature, the entire wavelet formalism can be intimidating, and may scare potential users despite potential advantages.

The article at hand provides a stream-lined development of the wavelet method proposed by Kaiser et al. [70], with a few highlights. We discussed the development of the used interpolating wavelets from scratch, dispensing with most of the wavelet machinery. As a byproduct, we could identify similarities to traditional centered finite-difference discretizations. With these developments at hand, we discussed the construction of a wavelet discretization associated to unit-cell problems in micromechanics in arbitrary dimension on a regular grid. The resemblance to the Moulinec-Suquet discretization, as developed from the "modern perspective" [23, 55] is striking: One starts with a discretization of the stored elastic energy, and obtains the balance equation as the first variation in a natural way. Similar to the Moulinec-Suquet discretization, there is only a single grid required to represent both the "nodal" degrees of freedom and their derivatives. This situation contrasts with most other discretizations, e.g., finite elements which distinguish a nodal grid and a quadrature-point grid. This finding highlights the role of the quadrature rule: Whereas the Moulinec-Suquet discretization is intricately linked with the trapezoidal rule – the same is true for the interpolating wavelets discretization on a regular grid – adaptive wavelet discretizations would be naturally adapted to locally refined quadrature rules of trapezoidal type. In particular, such approaches should benefit from the superconvergence properties [53, 54] which were recently found to drive the effectiveness of computational homogenization methods on regular grids. Also, the energetic formulation makes the symmetry of the ensuing discrete problem apparent, and would enable us to establish

<footnote n="1">¹ The term multiscale has a different meaning in this context from multiscale mechanics, for instance: The multiple scales refer to nested levels of discretization fineness.</footnote>
<figure><img src="image_15.png" /><figcaption>**Fig. 10** Macroscopic elastoplastic behavior of an isotropic fiber microstructure.</figcaption></figure>

lish similar contraction estimates for the basic scheme and other Lippmann-Schwinger type solvers that are familiar in FFT-based computational micromechanics. For completeness, we also provided the route to the associated Green's and Eshelby-Green operators. In the process, we discussed the emergence of the "undefined frequencies" for general DD wavelet order in a natural way.

We followed the path established by Kaiser et al. [70] and used the *discrete Fourier transform* (DFT) to construct Green's operator. However, the DFT requires a regular grid to be given. In particular, when realizing spatial adaptivity, a high-resolution background grid needs to be present, which, one the one hand, requires to be treated by computational means, and, on the other hand, may warrant some stabilization to retain computational robustness, as done in Kabel-Schneider [69] in a related scenario. It might be interesting to investigate whether the forward and backward wavelet transforms could be used directly to provide a "Green's operator", as done for $L^2$-wavelets [100, 101]. The authors are not aware of a result stating the required norm equivalence for the case of *interpolating wavelets*. Indeed, the latter were developed for spaces of *continuous functions* [71] for the simple reason that interpolation requires the function to have a well-defined value at a point. The cell problem is set in the first-order Sobolev space $H^1$, whose elements are discontinuous, in general, in spatial dimensions greater than one. Therefore, further research would be required to see if a robust "pure wavelet" Green's operator leads to robust iterative schemes. In a sense, this situation appears to be similar, where trigonometric interpolation is non-standard for discontinuous functions [55].

We reported on a number of computational studies whose purpose is to shed some light on properties of the wavelet discretizations not easily deduced from theoretical considerations. In particular, the three items (i)-(iii) described in the introduction are of particular importance here.

We observed that for moderate (elastic) material contrasts, i.e., between 0.01 and 100, the investigated wavelet discretizations showed an iteration count similar to traditional discretizations. This property was to be expected from the theoretical considerations. For higher contrast, however, things get more interesting. The required iteration counts for the DD2, DD3, and DD4 wavelet discretizations are strongly affected by the material contrast, similar to the original discretization via trigonometric polynomial proposed by Moulinec and Suquet. In contrast, the DD1 discretization converges faster and enables the associated solvers to converge on porous microstructures with a linear rate. In this respect, the behavior of DD1 wavelets is similar to the staggered and rotated staggered-grid discretizations. In the limiting case of infinite contrast, only the DD1 wavelet discretization showed a stable solver behavior. For the specific example considered, the DD1 discretization requires less iterations than both the staggered and the rotated staggered-grid approaches, although the local fields show more pronounced numerical artifacts.

We also conducted computational experiments for an inelastic material. As expected for the associated finite material contrast, the solvers performed similarly for all discretizations considered. Moreover, compared to trigonometric polynomials, staggered-grid, and rotated staggered-grid discretizations, the DD wavelet discretizations for all consid-
ered orders slightly overestimate the macroscopic stresses, both for the elastic and the inelastic case, similar to finite-element discretizations.

To sum up, we could not find convincing arguments for replacing established discretizations used in FFT-based methods by Deslauriers-Dubuc wavelets, at least on a regular grid. In a sense, this is not surprising as the DD wavelets on a regular grid come in two flavors. For low order, DD wavelets coincide with central difference discretizations, and their suitability for being used in FFT methods was investigated earlier. Deslauriers-Dubuc wavelets of high order, on the other hand, do not improve both the solver properties nor the convergence properties. Whether or not wavelet discretizations offer a real advantage over alternative approaches hinges on whether suitable preconditioning strategies based on the fast wavelet transform can be established for adaptive discretizations, dispensing with the FFT altogether. In fact, the competitors [61–63, 66] available promise higher convergence rates than adaptive wavelets. In three spatial dimensions, to achieve a specific fidelity in the effective stresses, the effort of traditional FFT methods scales as $N^3$, where N denotes the number of voxels per axis, whereas X-FFT leads to an $N^{3/2}$ scaling - similar to interface conforming finite elements. For adaptive methods on a regular grid, a scaling of $N^2$ is expected, proportional to the surface area. Thus, the scaling appears inferior. However, the salient scaling $N^{3/2}$ requires the interface to be sufficiently smooth, so adaptive methods may shine in the presence of geometric singularities, e.g., the tips of cylindrical fibers or the edges of polycrystalline grains/fillers.

# Appendix A Periodic homogenization in a nutshell

## Appendix A.1 The effective stiffness and the equilibrium equation

We consider a rectangular unit cell

$$Y = [0, L_1] \times [0, L_2] \times \dots \times [0, L_d] \qquad (A.1)$$

in d spatial dimensions with positive edge lengths $L_i$ ($i = 1, 2, \dots, d$). Suppose a free energy density

$$W : Y \times \mathrm{Sym}(d) \to \mathbb{R} \qquad (A.2)$$

is given which we assume to be measurable in the first slot and continuously differentiable in the second slot [31]. We assume that there is a constant $\alpha_+$, s.t. the estimate

$$\| \mathcal{S}(x, \boldsymbol{\varepsilon}_1) - \mathcal{S}(x, \boldsymbol{\varepsilon}_2) \| \le \alpha_+ \| \boldsymbol{\varepsilon}_1 - \boldsymbol{\varepsilon}_2 \| \qquad (A.3)$$

holds for (almost) all $x \in Y$ as well as $\boldsymbol{\varepsilon}_1, \boldsymbol{\varepsilon}_2 \in \mathrm{Sym}(d)$ and where we introduced the stress operator

$$\mathcal{S}(x, \boldsymbol{\varepsilon}) = \frac{\partial W}{\partial \boldsymbol{\varepsilon}}(x, \boldsymbol{\varepsilon}), \quad x \in Y, \quad \boldsymbol{\varepsilon} \in \mathrm{Sym}(d). \qquad (A.4)$$

We use the Frobenius norm

$$\| \boldsymbol{\varepsilon} \| = \sqrt{\boldsymbol{\varepsilon} : \boldsymbol{\varepsilon}}, \quad \boldsymbol{\varepsilon} \in \mathrm{Sym}(d), \qquad (A.5)$$

where the colon refers to a double contraction

$$\boldsymbol{\varepsilon} : \boldsymbol{\tau} = \sum_{i,j=1}^{d} \varepsilon_{ij} \tau_{ij} \quad \text{with} \quad \boldsymbol{\varepsilon} = \sum_{i,j=1}^{d} \varepsilon_{ij} \boldsymbol{e}_i \otimes \boldsymbol{e}_j \\ \text{and} \quad \boldsymbol{\tau} = \sum_{i,j=1}^{d} \tau_{ij} \boldsymbol{e}_i \otimes \boldsymbol{e}_j \qquad (A.6)$$

of indices with the canonical basis $(\boldsymbol{e}_1, \boldsymbol{e}_2, \dots, \boldsymbol{e}_d)$ of the Cartesian space $\mathbb{R}^d$.

An important special case arises in case a (measurable) microscopic distribution of stiffness tensors

$$C : Y \to L(\mathrm{Sym}(d)), \quad x \mapsto C(x), \qquad (A.7)$$

is given, which we consider pointwise as a linear mapping on the vector space $\mathrm{Sym}(d)$ of symmetric $d \times d$-tensors. We assume that the stiffness distribution (A.7) is symmetric

$$[C(x) : \boldsymbol{\varepsilon}_1] : \boldsymbol{\varepsilon}_2 = [C(x) : \boldsymbol{\varepsilon}_2] : \boldsymbol{\varepsilon}_1 \qquad (A.8)$$

for all strains $\boldsymbol{\varepsilon}_1, \boldsymbol{\varepsilon}_2 \in \mathrm{Sym}(d)$ and almost every continuum point $x \in Y$ in the cell (A.1). Then, we may consider the Hookean free energy

$$W(x, \boldsymbol{\varepsilon}) = \frac{1}{2} \boldsymbol{\varepsilon} : C(x) : \boldsymbol{\varepsilon}, \quad x \in Y, \quad \boldsymbol{\varepsilon} \in \mathrm{Sym}(d), \qquad (A.9)$$

with the framework (A.2). General inelastic problems are treated via progressing time steps, and the free energy (A.2) is taken as an incremental energy [102, 103].

In periodic homogenization [104–106], we are interested in the periodic displacement fluctuation

$$\boldsymbol{u} \in H_{\#}^{1}(Y)^{d} \qquad (A.10)$$

with vanishing mean in the first-order Sobolev space which solves the quasistatic balance of linear momentum

$$\mathrm{div} \mathcal{S}(\cdot, \bar{\boldsymbol{\varepsilon}} + \nabla^s \boldsymbol{u}) = \mathbf{0} \qquad (A.11)$$

on the microscopic scale for a prescribed macroscopic strain $\bar{\boldsymbol{\varepsilon}} \in \mathrm{Sym}(d)$. The symbol $\nabla^s$ stands for the symmetrized gradient operator. Under suitable conditions on the free energy
(A.2), existence and uniqueness of a solution (A.10) to the equilibrium equation (A.11) may be established [31]. The simplest such condition is (uniform) strong monotonicity of the stress operator (A.4), i.e., the existence of a positive constant, s.t., the inequality

$$[\mathcal{S}(\mathbf{x}, \boldsymbol{\epsilon}_1) - \mathcal{S}(\mathbf{x}, \boldsymbol{\epsilon}_2)] : (\boldsymbol{\epsilon}_1 - \boldsymbol{\epsilon}_2) \geq \alpha_- \|\boldsymbol{\epsilon}_1 - \boldsymbol{\epsilon}_2\|^2 \quad (A.12)$$

holds for all strains $\boldsymbol{\epsilon}_1, \boldsymbol{\epsilon}_2 \in \text{Sym}(d)$ and almost all points $\mathbf{x} \in Y$. Physically speaking the condition (A.12) rules out degeneracies like voids/pores and rigid inclusions. Such cases can be treated under suitable conditions on the arrangement of the phases, and we refer to the pertinent literature [31, 32] for a discussion. Moreover, the case of negative eigenvalues in the material tangent leads to a violation of condition (A.12).

For prescribed macroscopic strain $\bar{\boldsymbol{\epsilon}}$, we denote by

$$\sigma_{\bar{\boldsymbol{\epsilon}}} = \mathcal{S}(\cdot, \bar{\boldsymbol{\epsilon}} + \nabla^s \mathbf{u}_{\bar{\boldsymbol{\epsilon}}}) \quad (A.13)$$

the associated local stress field, where $\mathbf{u}_{\bar{\boldsymbol{\epsilon}}} \in H_{\#}^1(Y)^d$ stands for the unique solution to the balance equation (A.11). Of particular significance is the average stress response

$$\bar{\boldsymbol{\sigma}}_{\bar{\boldsymbol{\epsilon}}} = \langle \boldsymbol{\sigma}_{\bar{\boldsymbol{\epsilon}}} \rangle_Y \quad (A.14)$$

to the applied strain $\bar{\boldsymbol{\epsilon}}$ which arises by spatially averaging the local stress field, i.e., we use the definition

$$\langle \phi \rangle_Y = \frac{1}{\text{vol}(Y)} \int_Y \phi(\mathbf{x}) \, \text{dV}(\mathbf{x}) \quad (A.15)$$

for the mean of any (integrable) field $\phi$ on the cell (A.1).

In the special case (A.9) of linear elasticity, the effective stress $\bar{\boldsymbol{\sigma}}_{\bar{\boldsymbol{\epsilon}}}$ depends linearly on the applied strain $\bar{\boldsymbol{\epsilon}}$, and we may define the effective stiffness tensor

$$C^{\text{eff}} \in L(\text{Sym}(d)) \quad (A.16)$$

implicitly via the relation

$$C^{\text{eff}} : \bar{\boldsymbol{\epsilon}} = \bar{\boldsymbol{\sigma}}_{\bar{\boldsymbol{\epsilon}}}. \quad (A.17)$$

In practice, the effective stiffness (A.16) is evaluated by selecting a linear basis of the vector space $\text{Sym}(d)$ and solving the balance equation (A.11) once for every basis tensor. Thus, three solutions are required in two spatial dimensions, whereas six cases must be considered in three spatial dimensions, in the general scenario. In case the microstructure (A.7) has additional symmetries, sometimes fewer solutions are required [107].

The effective stiffness (A.16) is symmetric in the sense (A.8) and inherits the bounds (A.3) and (A.12) from the local

stiffness. These properties are most easily seen by considering the variational principle

$$\mathcal{E}_{\bar{\boldsymbol{\epsilon}}}(\mathbf{u}) := \langle W(\cdot, \bar{\boldsymbol{\epsilon}} + \nabla^s \mathbf{u}) \rangle_Y \rightarrow \min_{\mathbf{u} \in H_{\#}^1(Y)^d} \quad (A.18)$$

which seeks periodic displacement fields minimizing the average free energy (A.2). The associated critical points coincide with (weak) solutions of the balance equation, even in the general case (A.2). For linear elasticity, the effective stiffness (A.17) emerges naturally

$$\frac{1}{2} \bar{\boldsymbol{\epsilon}} : C^{\text{eff}} : \bar{\boldsymbol{\epsilon}} = \min_{\mathbf{u} \in H_{\#}^1(Y)^d} \mathcal{E}_{\bar{\boldsymbol{\epsilon}}}(\mathbf{u}) \quad (A.19)$$

via energy equivalence [54].

Except for a few special cases, see Milton's treatise [108], no analytical solutions for the equilibrium equation (A.11) are available. Therefore, it is imperative to resort to computational approaches to resolve the balance equation (A.11). There are two essential steps for realizing this task: First, the balance equation (A.11) needs to be suitably discretized, i.e., turned into a finite-dimensional problem. Subsequently, the discretized problems need to be solved, or, to be more precise, the solution needs to be approximated to the desired accuracy.

Discretizing the problem (A.11) accurately and solving it efficiently afterwards are contradictory goals. Indeed, it is not difficult to produce false solutions quickly. To offer a compromise and to inherit the positive properties of the continuum formulation, it is convenient to follow the pioneering ideas of Moulinec-Suquet [1, 2] and to study potential solvers in a continuum setting first, see the next section. In such a way one may differentiate between the characteristics of the continuum problem and the properties of the discretization scheme.

## Appendix A.2 The Lippmann-Schwinger equation and associated solvers

To treat the balance equation (A.11) numerically, it is convenient to introduce a homogeneous reference stiffness $C^0$, i.e., a symmetric and positive definite element of the space $L(\text{Sym}(d))$ which will later serve as a numerical parameter. With its help, we may re-write the balance equation (A.11) in the form

$$\text{div}C^0 : \nabla^s \mathbf{u} = - \div \boldsymbol{\tau} \quad (A.20)$$

with the stress polarization

$$\boldsymbol{\tau} = \mathcal{S}(\cdot, \bar{\boldsymbol{\epsilon}} + \nabla^s \mathbf{u}) - C^0 : (\bar{\boldsymbol{\epsilon}} + \nabla^s \mathbf{u}). \quad (A.21)$$
Considering the left-hand side in the eq. (A.20) as an operator with inverse $G^0$, we may consider the equivalent form of eq. (A.20)

$$ u = -G^0 \div \tau \equiv -G^0 \operatorname{div} [\mathcal{S}(\cdot, \bar{\varepsilon} + \nabla^s u) - C^0 : (\bar{\varepsilon} + \nabla^s u)], \quad (A.22) $$

the displacement form of the celebrated Lippmann-Schwinger equation [4–6]. More often than not, eq. (A.22) is differentiated, adding the macroscopic strain to both sides in the process, resulting in the equation

$$ \bar{\varepsilon} + \nabla^s u = \bar{\varepsilon} - \nabla^s G^0 \operatorname{div} [\mathcal{S}(\cdot, \bar{\varepsilon} + \nabla^s u) - C^0 : (\bar{\varepsilon} + \nabla^s u)], \quad (A.23) $$

which may be re-written in the form

$$ \varepsilon = \bar{\varepsilon} - \Gamma^0 : [\mathcal{S}(\cdot, \varepsilon) - C^0 : \varepsilon] \quad (A.24) $$

with the (total) strain field

$$ \varepsilon = \bar{\varepsilon} + \nabla^s u \in L^2(Y; \operatorname{Sym}(d)) \quad (A.25) $$

and the Eshelby-Green operator

$$ \Gamma^0 = \nabla^s G^0 \operatorname{div} \in L(L^2(Y; \operatorname{Sym}(d))). \quad (A.26) $$

Equation (A.24) is known as the Lippmann-Schwinger equation in strain form. Both Lippmann-Schwinger equations (A.22) and (A.24) are formally equivalent, i.e., every displacement solution of the eq. (A.22) gives rise to a solution to eq. (A.24) via the prescription (A.25). In turn, every strain field $\varepsilon$ solving eq. (A.24) has the form (A.25) with a displacement field $u \in H_\#^1(Y)^d$ solving the displacement-form equation (A.22). Moreover, every displacement solution to the Lippmann-Schwinger equation (A.22) solves the original balance equation (A.11) and vice versa. This property, albeit elementary, is critical, as it emphasizes the role of the reference stiffness as a numerical parameter. Indeed, usually a parameter involved in an equation leads to a solution depending on the parameter. In our case, no such dependence is observed.

Either of the Lippmann-Schwinger equations has a suggestive fixed-point form, i.e., leading to the natural fixed-point scheme

$$ \varepsilon^{k+1} = \bar{\varepsilon} - \Gamma^0 : [\mathcal{S}(\cdot, \varepsilon^k) - C^0 : \varepsilon^k], \quad k = 0, 1, 2, \dots, \quad (A.27) $$

where $k$ is the iteration index and some initial guess $\varepsilon^0$ must be provided to start the iterative scheme. We restrict to the strain form (A.24), as this formulation typically leads to faster schemes, whereas the displacement formulation (A.22) requires less memory, in general. We refer to the works [36–38, 87] for more details.

The approach (A.27) was introduced by Moulinec-Suquet [1, 2] as their basic scheme. To turn this strategy into an iterative scheme, two questions must be answered first. For a start, the operator $\Gamma^0$ – which we introduced rather formally just yet via a concatenation of

operators (A.26) – must be evaluated. Moreover, a recommendation for the reference material should be given.

To evaluate Green’s operator $G^0$, it is convenient to exploit the periodicity contained in the definition (A.10) of the displacement fluctuation and to express the field in terms of the Fourier series

$$ u(x) = \sum_{\xi \in \mathbb{Z}^d} \hat{u}(\xi) e^{i \tilde{\xi} \cdot x}, \quad \tilde{\xi}_k = \frac{2\pi}{L_k} \xi_k, \quad k = 1, 2, \dots, d, \quad (A.28) $$

with suitable Fourier coefficients $(\hat{u}(\xi))_{\xi \in \mathbb{Z}^d}$ and the imaginary unit $i$. Then, the left-hand side of eq. (A.20) takes the form

$$ \operatorname{div} C^0 : \nabla^s u(x) = - \sum_{\xi \in \mathbb{Z}^d} \tilde{\xi} \cdot [C^0 : (\tilde{\xi} \otimes_s \hat{u}(\xi))] e^{i \tilde{\xi} \cdot x}. \quad (A.29) $$

Thus, one is led to solve a linear equation

$$ - \tilde{\xi} \cdot [C^0 : (\tilde{\xi} \otimes_s \hat{u}(\xi))] = \hat{f}(\xi) \quad (A.30) $$

for every non-zero frequency vector $\xi \in \mathbb{Z}^d \setminus \{0\}$ and where $(\hat{f}(\xi))$ refer to the Fourier coefficients of the right-hand side $f = -\operatorname{div}\tau$ in eq. (A.20). Solving the linear system (A.30) provides the action of Green’s operator $G^0$ in Fourier space. Both the symmetrized gradient operator $\nabla^s$ and the divergence $\operatorname{div}$ may be expressed in Fourier space. Therefore, the entire Eshelby-Green operator (A.26) may be evaluated efficiently in Fourier space.

Concerning the choice of the reference material, a simple choice in the scenario (A.12) and (A.3) selects the reference material in the form

$$ C^0 = \alpha_0 I, \quad \alpha_0 = \frac{\alpha_- + \alpha_+}{2} \quad (A.31) $$

with the identity tensor $I$ on symmetric $d \times d$-tensors. Then, for prescribed fidelity, the iteration count required for the basic scheme (A.27) scales linearly with the material contrast

$$ \kappa = \frac{\alpha_+}{\alpha_-}. \quad (A.32) $$

Interestingly, as these considerations were made for the continuous formulation, dedicated discretizations preserving the outlined structure lead to the same iteration bound. In particular, a resolution-independent convergence behavior may be expected here. This property contrasts with traditional approaches for solving the balance equation (A.11), where an increase in the condition number of the associated linear system is observed for finer and finer meshes. In this sense, working with the Lippmann-Schwinger formulation (A.26) leads to a natural preconditioning [31, 38] which provides resolution-independent convergence guarantees at the expense of a material-contrast dependence of the iteration-count bounds.

Under more restrictive hypotheses on the symmetries of the local stiffnesses, more clever choices for the reference material are possible [109, 110]. However, these alterations offer only a minor improvement due to intrinsic limitations of the basic scheme (A.27). In fact, the basic scheme (A.27) may be interpreted as
a gradient-descent method for the average elastic energy (A.18), and the limitations of gradient methods are well-known and well-reported in the optimization literature [111]. As a remedy, more efficient solution schemes may be considered, which also provide resolution-independent convergence guarantees but require fewer iterations than the basic scheme (A.27). For linear problems, the conjugate gradient (CG) method, introduced by Brisard-Dormieux [16] and Zeman et al. [91] for this type of problems, serves as the gold standard for other solvers. The iteration count scales with the square root $\sqrt{\kappa}$ of the material contrast (A.32), which is known to be optimal [111], in general. There are also other strategies which are not as fast, but may be implemented with less memory [90, 99, 112].

## Appendix B Interpolated derivatives of interpolating wavelets

The purpose of this appendix is to provide the arguments for the validity of the representation (2.30) and the values recorded in Tab. 2. We consider a function

$$f = \sum_{k \in \mathbb{Z}} s_k^\ell \phi_k^\ell \in V^\ell \quad (B.1)$$

and compute the exact derivative

$$f' = \sum_{k \in \mathbb{Z}} s_k^\ell \frac{d\phi_k^\ell}{dx}. \quad (B.2)$$

To retain the space $V^\ell$, we "project" the derivative (B.2) back onto the space $V^\ell$ via the interpolation operator (2.19)

$$\tilde{\partial} f := \mathcal{P}^\ell [f'] . \quad (B.3)$$

A special role here play the values of the derivative of the fundamental function (2.1) at the integers

$$r_k := \phi'(k), \quad k \in \mathbb{Z}, \quad (B.4)$$

in case the fundamental function (2.1) is differentiable. Due to the compact support, only a finite number of the coefficients $\{r_k\}$ are actually non-zero.

By the linearity of the projection (2.19) and the scaling property (2.14), we obtain the explicit description

$$\mathcal{P}^\ell [f'] = 2^\ell \sum_{k \in \mathbb{Z}} \sum_{m \in \mathbb{Z}} s_m^\ell r_{k-m} \phi_k^\ell = 2^\ell \sum_{k \in \mathbb{Z}} \sum_{j \in \mathbb{Z}} s_{k+j}^\ell r_j \phi_k^\ell \quad (B.5)$$

for the derivative (B.3). The last two representations are related by the index change $j = k - m$. In particular, if we express the derivative

$$\mathcal{P}^\ell [f'] = \sum_{k \in \mathbb{Z}} \tilde{s}_k^\ell \phi_k^\ell \in V^\ell \quad (B.6)$$

in the natural basis with coefficients $\{\tilde{s}_k^\ell\}$, the latter may be written in the explicit form

$$\tilde{s}_k^\ell = 2^\ell \sum_{j \in \mathbb{Z}} s_{k+j}^\ell r_j. \quad (B.7)$$

The latter formula shows that the derivative operator (B.3) on the approximation space $V^\ell$ may be regarded as a special finite-difference operator acting on the coefficients $\{s_k^\ell\}$ used to represent the function (B.1).

It remains to identify the values of the coefficients (B.4). For this purpose, the representation (2.10) is helpful, i.e., following the literature [70, 83], we observe

$$\begin{align}
r_k &= \phi'(k) = - \int_{-\infty}^{\infty} \Phi(y) \Phi'(y-k) \\
dy &= - \int_{-\infty}^{\infty} \Phi(z+k) \Phi'(z) \\
dz &= \int_{-\infty}^{\infty} \Phi(z-k) \Phi'(z) dz, \quad k \in \mathbb{Z}, \tag{B.8}
\end{align}$$

where we used the substitution $z = y - k$ and the symmetry $r_k = -r_{-k}$, $k \in \mathbb{Z}$, implied by the symmetry $\phi(-x) = \phi(x)$, $x \in \mathbb{R}$, of the fundamental function. The last expression in eq. (B.8) was computed by Beylkin [113] for the Daubechies wavelets. For the considered Deslauriers-Dubuc wavelets of low order, we list the following coefficients in Tab. 2. Beylkin [113] also provided a formal computation of the coefficients for the case $M = 1$, corresponding to DD1 wavelets, which are not classically differentiable. We included the corresponding coefficients in the table, as they are useful for the hat-shaped Deslauriers-Dubuc wavelets of order $M = 1$.
## Appendix C Pseudo-codes

For the convenience of the reader, this appendix comprises pseudo-codes for the relevant algorithms: the basic scheme, Alg. 1, the Barzilai-Borwein method, Alg. 2, and the non-linear conjugate gradient algorithm, Alg. 3.

**Algorithm 1 Basic scheme [1, 2] (maxit, tol)**

1: Determine initial guess $\epsilon$ (with the correct mean) and initial reference material $C^0$

2: $\text{res} \leftarrow \text{tol} + 1$

3: $k \leftarrow 0$

4: **while** $k < \text{maxit}$ **and** $\text{res} > \text{tol}$ **do**

5: $k \leftarrow k + 1$

6: $\boldsymbol{\epsilon}^{\text{old}} \leftarrow \boldsymbol{\epsilon}$

7: $\boldsymbol{\epsilon} \leftarrow \mathcal{S}(\boldsymbol{\epsilon}) - C^0 : \boldsymbol{\epsilon}$ ▷ Compute stress polarization (A.21), apply constitutive law (A.4)

8: $\hat{\boldsymbol{\epsilon}} \leftarrow \text{FFT}(\boldsymbol{\epsilon})$

9: $\bar{\boldsymbol{\sigma}} \leftarrow \hat{\boldsymbol{\epsilon}}(0) + C^0 : \bar{\boldsymbol{\epsilon}}$

10: $\hat{\boldsymbol{\epsilon}}(\xi) \leftarrow \begin{cases} \bar{\boldsymbol{\epsilon}}, \quad \xi = 0 \\ -\hat{\Gamma}^0(\xi) : \hat{\boldsymbol{\epsilon}}(\xi), \quad \text{otherwise} \end{cases}$ ▷ See eq. (2.74)

11: $\boldsymbol{\epsilon} \leftarrow \text{FFT}^{-1}(\hat{\boldsymbol{\epsilon}})$

12: $r \leftarrow \langle \|\boldsymbol{\epsilon} - \boldsymbol{\epsilon}^{\text{old}}\|^2 \rangle_Y$

13: $\text{res} \leftarrow \sqrt{r} / \|\bar{\boldsymbol{\sigma}}\|$

14: **end while**

15: **return** $\boldsymbol{\epsilon}$, $\text{res}$

**Algorithm 2 Barzilai-Borwein scheme [90] (maxit, tol)**

1: Determine initial guess $\epsilon$ (with the correct mean) and initial reference material $C^0$

2: $\text{res} \leftarrow 1$

3: $k \leftarrow 0$

4: $\rho \leftarrow 1$

5: $\mathbf{R} \leftarrow 0$

6: **while** $n < \text{maxit}$ **and** $\text{res} > \text{tol}$ **do**

7: $k \leftarrow k + 1$

8: $\mathbf{R} \leftarrow \boldsymbol{\epsilon}$

9: $\begin{bmatrix} \mathbf{R} \\ q \end{bmatrix} \leftarrow \begin{bmatrix} \mathcal{S}(\mathbf{R}) \\ \langle \mathcal{S}(\mathbf{R}) : \mathbf{R} \rangle_Y \end{bmatrix}$ ▷ Apply constitutive law (A.4)

10: $\bar{\boldsymbol{\sigma}} \leftarrow \langle \mathbf{R} \rangle_Y$

11: $\hat{\mathbf{R}} \leftarrow \text{FFT}(\mathbf{R})$

12: $\hat{\mathbf{R}}(\xi) \leftarrow \begin{cases} 0, & \xi = 0 \\ \hat{\Gamma}^0(\xi) : \hat{\mathbf{R}}(\xi), & \text{otherwise} \end{cases}$ ▷ See eq. (2.74)

13: $\mathbf{R} \leftarrow \text{FFT}^{-1}(\hat{\mathbf{R}})$

14: $\rho \leftarrow \rho/(1 - q/r)$

15: $\boldsymbol{\epsilon} \leftarrow \boldsymbol{\epsilon} - \rho \mathbf{R}$

16: $r \leftarrow \langle \|\mathbf{R}\|^2 \rangle_Y$

17: $\text{res} \leftarrow \sqrt{r} / \|\bar{\boldsymbol{\sigma}}\|$

18: **end while**

19: **return** $\boldsymbol{\epsilon}$, $\text{res}$

**Algorithm 3 Non-linear CG method [99] (maxit, tol)**

1: Determine initial guess $\epsilon$ (with the correct mean) and initial reference material $C^0$

2: $\text{res} \leftarrow \text{tol} + 1$

3: $k \leftarrow 0$

4: $\gamma \leftarrow 1$

5: $\mathbf{D} \leftarrow 0$

6: **while** $k < \text{maxit}$ **and** $\text{res} > \text{tol}$ **do**

7: $k \leftarrow k + 1$

8: $\mathbf{G} \leftarrow \mathcal{S}(\boldsymbol{\epsilon})$ ▷ Apply constitutive law (A.4)

9: $\bar{\boldsymbol{\sigma}} \leftarrow \langle \mathbf{G} \rangle_Y$

10: $\hat{\mathbf{G}} \leftarrow \text{FFT}(\mathbf{G})$

11: $\hat{\mathbf{G}}(\xi) \leftarrow \begin{cases} \mathbf{0}, & \xi = 0 \\ \hat{\Gamma}^0(\xi) : \hat{\mathbf{G}}(\xi), & \text{otherwise} \end{cases}$ ▷ See eq. (2.74)

12: $\mathbf{G} \leftarrow \text{FFT}^{-1}(\hat{\mathbf{G}})$

13: $\gamma^{\text{old}} \leftarrow \gamma$

14: $\gamma \leftarrow \|\mathbf{G}\|^2$

15: $\beta \leftarrow \gamma / \gamma^{\text{old}}$

16: $\mathbf{D} \leftarrow -\mathbf{G} + \beta \mathbf{D}$

17: $\boldsymbol{\epsilon} \leftarrow \boldsymbol{\epsilon} + \alpha \mathbf{D}$

18: $\text{res} \leftarrow \sqrt{\gamma / \|\bar{\boldsymbol{\sigma}}\|}$

19: **end while**

20: **return** $\boldsymbol{\epsilon}$, $\text{res}$

**Funding** Open Access funding enabled and organized by Projekt DEAL.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

**Acknowledgements** Funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – 507778349 – is gratefully acknowledged by ED. MS acknowledges support from the European Research Council within the Horizon Europe program - project 101040238. We thank the anonymous reviewers for their constructive comments and suggestions.
References

1. Moulinec H, Suquet P (1994) "A fast numerical method for computing the linear and nonlinear mechanical properties of composites," *Comptes Rendus de l'Académie des Sciences. Série II* 318(11):1417-1423

2. Moulinec H, Suquet P (1998) A numerical method for computing the overall response of nonlinear composites with complex microstructure. *Comput Methods Appl Mech Eng* 157:69-94

3. Matouš K, Geers MGD, Kouznetsova VG, Gillman A (2017) A review of predictive nonlinear theories for multiscale modeling of heterogeneous materials. *J Comput Phys* 330:192-220

4. Zeller R, Dederichs PH (1973) "Elastic constants of polycrystals," *physica status solidi*, 55(2), 831-842

5. Kröner E (1977) Bounds for effective elastic moduli of disordered materials. *J Mech Phys Solids* 25(2):137-155

6. Mura T (1987) *Micromechanics of Defects in Solids*. Martinus Nijhoff, Dordrecht

7. Lebensohn RA, Rollett AD (2020) Spectral methods for full-field micromechanical modelling of polycrystalline material. *Comput Mater Sci* 173:109336

8. Segurado J, Lebensohn RA, LLorca J (2018) "Chapter One - Computational Homogenization of Polycrystals," *Advances in Applied Mechanics*, vol. 51, pp. 1-114

9. Schneider M (2021) A review of non-linear FFT-based computational homogenization methods. *Acta Mech* 232:2051-2100

10. Müller WH (1996) "Mathematical vs. Experimental Stress Analysis of Inhomogeneities in Solids," *Journal de Physique IV*, vol. 6, pp. C1.139-C1.148

11. Müller WH (1998) "Fourier transforms and their application to the formation of textures and changes of morphology in solids," in *IUTAM Symposium on Transformation Problems in Composite and Active Materials* (Y. A. Bahei-El-Din and G. J. Dvorak, eds.), vol. 60, (Dordrecht), pp. 61-72, Kluwer

12. Dreyer W, Müller WH, Olschewski J (1999) An approximate analytical 2D-solution for the stresses and strains in eigenstrained cubic materials. *Acta Mech* 136(3-4):171-192

13. Bonnet G (2007) Effective properties of elastic periodic composite media with fibers. *J Mech Phys Solids* 55:881-899

14. Monchiet V (2015) Combining FFT methods and standard variational principles to compute bounds and estimates for the properties of elastic composites. *Comput Methods Appl Mech Eng* 283:454-473

15. Vondřejc J (2014) Improved guaranteed computable bounds on homogenized properties of periodic media by Fourier-Galerkin method with exact integration. *Int J Numer Meth Eng* 107:1106-1135

16. Brisard S, Dormieux L (2010) FFT-based methods for the mechanics of composites: A general variational framework. *Comput Mater Sci* 49(3):663-671

17. Brisard S, Dormieux L (2012) Combining Galerkin approximation techniques with the principle of Hashin and Shtrikman to derive a new FFT-based numerical method for the homogenization of composites. *Comput Methods Appl Mech Eng* 217-220:197-212

18. Hashin Z, Shtrikman S (1962) On some variational principles in anisotropic and nonhomogeneous elasticity. *J Mech Phys Solids* 10(4):335-342

19. Hashin Z, Shtrikman S (1962) A variational approach to the theory of the elastic behaviour of polycrystals. *J Mech Phys Solids* 10(4):343-352

20. Hashin Z, Shtrikman S (1963) A variational approach to the theory of the elastic behaviour of multiphase materials. *J Mech Phys Solids* 11(2):127-140

21. Schneider M (2019) On the mathematical foundations of the self-consistent clustering analysis for non-linear materials at small strains. *Comput Methods Appl Mech Eng* 354:783-801

22. Nguyen N-T, Licht C, Kweon J-H (2011) An efficient homogenization method using the trigonometric interpolation and the fast Fourier transform. *Vietnam Journal of Mechanics* 33(4):215-223

23. Vondřejc J, Zeman J, Marek I (2014) An FFT-based Galerkin method for homogenization of periodic media. *Comput Math Appl* 68:156-173

24. Schneider M (2015) Convergence of FFT-based homogenization for strongly heterogeneous media. *Mathematical Methods in the Applied Sciences* 38(13):2761-2778

25. Wiegmann A, Zemitis A (2006) *EJ-HEAT: A fast explicit jump harmonic averaging solver for the effective heat conductivity of composite materials*. *Berichte des Fraunhofer ITWM* 94:1-21

26. Willot F, Pellegrini YP (2008) "Fast Fourier transform computations and build-up of plastic deformation in 2D, elastic-perfectly plastic, pixelwise disordered porous media," in *Continuum Models and Discrete Systems* (D. Jeulin and S. Forest, eds.), vol. CMDS11, (Paris), pp. 443-449, Presse Ecole des Mines Paris

27. Willot F, Abdallah B, Pellegrini Y-P (2014) *Fourier-based schemes with modified Green operator for computing the electrical response of heterogeneous media with accurate local fields*. *Int J Numer Meth Eng* 98:518-533

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

29. Gibbs JW (1898) *Fourier's Series*. Nature 59(1522):200

30. Schneider M, Ospald F, Kabel M (2016) Computational homogenization of elasticity on a staggered grid. *Int J Numer Meth Eng* 105(9):693-720

31. Schneider M (2020) *Lippmann-Schwinger solvers for the computational homogenization of materials with pores*. *Int J Numer Meth Eng* 121(22):5017-5041

32. Sab K, Bleyer J, Brisard S, Dolbeau M (2024) An FFT-based adaptive polarization method for infinitely contrasted media with guaranteed convergence. *Comput Methods Appl Mech Eng* 427:117012

33. Willot F (2015) *Fourier-based schemes for computing the mechanical response of composites with accurate local fields*. *Comptes Rendus Mécanique* 343:232-245

34. Dorn C, Schneider M (2019) *Lippmann-Schwinger solvers for the explicit jump discretization for thermal computational homogenization problems*. *Int J Numer Meth Eng* 118(11):631-653

35. Ernesti F, Schneider M (2021) A fast Fourier transform based method for computing the effective crack energy of a heterogeneous material on a combinatorially consistent grid. *Int J Numer Meth Eng* 122:6283-6307

36. Schneider M, Merkert D, Kabel M (2017) FFT-based homogenization for microstructures discretized by linear hexahedral elements. *Int J Numer Meth Eng* 109:1461-1489

37. Leuschner M, Fritzen F (2018) *Fourier-Accelerated Nodal Solvers (FANS)* for homogenization problems. *Comput Mech* 62:359-392

38. Ladecký M, Leute RJ, Falsafi A, Pultarová I, Pastewka L, Junge T, Zeman J (2023) An optimal preconditioned FFT-accelerated finite element solver for homogenization. *Appl Math Comput* 446:127835

39. Berbenni S, Taupin V, Djaka KS, Fressengeas C (2018) A numerical spectral approach for solving elasto-static Field Dislocation and G-Dislocation Mechanics. *International Journal for Solids and Structures* 51:4157-4175

40. Lebensohn RA, Needleman A (2016) Numerical implementation of non-local polycrystal plasticity using fast Fourier transforms. *J Mech Phys Solids* 97:333-351
41. Vidyasagar A, Tan WL, Kochmann DM (2017) Predicting the effective response of bulk polycrystalline ferroelectric ceramics via improved spectral phase field methods. Mechanics and Physics of Solids 106:133–151

42. Vidyasagar A, Tutcuoglu AD, Kochmann DM (2018) Deformation patterning in finite-strain crystal plasticity by spectral homogenization with application to magnesium. Comput Methods Appl Mech Eng 335:584–609

43. Eloh KS, Jacques A, Berbenni S (2019) Development of a new consistent discrete Green operator for FFT-based methods to solve heterogeneous problems with eigenstrains. Int J Plast 116:1–23

44. Tu F, Jiao Y, Zhou X, Cheng Y, Tan F (2020) The implementation of B-splines to Hashin and Shtrikman variational principle based FFT method for the homogenization of composite. Int J Solids Struct 191–192:133–145

45. Derraz MR, Boukour M, El Mamouni A, El Omri A (2022) "Three-Dimensional Finite Radon Transform and Linear Homogenization," in Advanced Intelligent Systems for Sustainable Development (AI2SD'2020) (J. Kycprzyk, V. E. Balas, and M. Ezziyyani, eds.), vol. 1418 of Advances in Intelligent Systems and Computing, pp. 164–182, Springer International Publishing

46. Jabs L, Schneider M (2024) A consistent discretization via the finite radon transform for FFT-based computational micromechanics. Comput Mech 75:1451–1470

47. Finel A (2025) A tetrahedron-based discretization for FFT-based computational homogenization with smooth solution fields. Comput Methods Appl Mech Eng 436:117703

48. Gehrig F, Schneider M (2025) "Element-based internal variable formulations for finite element discretizations in FFT-based homogenization methods," International Journal for Numerical Methods in Engineering, vol. Online, pp. 1–34,

49. Gélébart L, Ouaki F (2015) Filtering material properties to improve FFT-based methods for numerical homogenization. J Comput Phys 294:90–95

50. Kabel M, Merkert D, Schneider M (2015) Use of composite voxels in FFT-based homogenization. Comput Methods Appl Mech Eng 294:168–188

51. Keshav S, Fritzen F, Kabel M (2022) FFT-based homogenization at finite strains using composite boxels (ComBo). Comput Mech 71:191–212

52. Ye C, Chung ET (2023) Convergence of trigonometric and finite-difference discretization schemes for FFT-based computational micromechanics. BIT Numer Math 63:11

53. Schneider M, Wicht D (2022) Superconvergence of the effective Cauchy stress in computational homogenization of inelastic materials. Int J Numer Methods Eng. 124(4): 959–978

54. Schneider M (2022) Superaccurate effective elastic moduli via postprocessing in computational homogenization. Int J Numer Meth Eng 123(17):4119–4135

55. Schneider M (2023) On the effectiveness of the Moulinec-Suquet discretization for composite materials. Int J Numer Meth Eng 124(14):3191–3218

56. Lendvai J, Schneider M (2024) Assumed strain methods in micromechanics, laminate composite voxels and level sets. Int J Numer Meth Eng 125(11):e7459

57. Kochmann J, Manjunatha K, Gierden C, Wulfinghoff S, Svendsen B, Reese S (2019) A simple and flexible model order reduction method for FFT-based homogenization problems using a sparse sampling technique. Comput Methods Appl Mech Eng 347:622–638

58. Gierden C, Waimann J, Svendsen B, Reese S (2021) FFT-based simulation using a reduced set of frequencies adapted to the underlying microstructure. Computer Methods in Materials Science 21(1):51–58

59. Waimann J, Gierden C, Schmidt A, Svendsen B, Reese S (2022) Reduced FFT-Based Simulation of a Mechanically Loaded Clus-

tered Microstructure Using an Adaptive Set of Fourier Modes. Key Eng Mater 926:2285–2292

60. Lahellec N, Michel JC, Moulinec H, Suquet P (2003) "Analysis of Inhomogeneous Materials at Large Strains using Fast Fourier Transforms," in IUTAM Symposium on Computational mechanics of Solid Materials at Large Strains (C. Miehe, ed.), vol. 108 of Solid Mechanics and Its Applications, pp. 247–258, Springer Netherlands

61. Zecevic M, Lebensohn RA, Capolungo L (2022) New large-strain FFT-based formulation and its application to model strain localization in nano-metallic laminates and other strongly anisotropic crystalline materials. Mech Mater 166:104208

62. Bellis C, Ferrier R (2024) Numerical homogenization by an adaptive Fourier spectral method on non-uniform grids using optimal transport. Comput Methods Appl Mech Eng 419:116658

63. Zecevic M, Lebensohn RA, Capolungo L (2026) Achieving geometric accuracy in FFT-based micromechanical models using conformal grid. Mech Mater 212:105512

64. Moës N, Cloirec M, Cartraud P, Remacle J-F (2003) A computational approach to handle complex microstructure geometries. Comput Methods Appl Mech Eng 192(28-30):3163–3177

65. Lian W-D, Legrain G, Cartraud P (2013) Image-based computational homogenization and localization: comparison between X-FEM/levelset and voxel-based approaches. Comput Mech 51(3):279–293

66. Gehrig F, Schneider M (2025) An X-FFT Solver for Two-Dimensional Thermal Homogenization Problems. Int J Numer Meth Eng 126(7):e70022

67. Eyre DJ, Milton GW (1999) A fast numerical scheme for computing the response of composites using grid refinement. The European Physical Journal - Applied Physics 6(1):41–47

68. Nkoumbou Kaptchouang NB, Gélébart L (2022) "Multiscale coupling of FFT-based simulations with the LDC approach," Computer Methods in Applied Mechanics and Engineering, vol. 394, p. 114921

69. Kabel M, Schneider M (2024) Adaptive material evaluation by stabilized octree and sandwich coarsening in FFT-based computational micromechanics. Int J Numer Meth Eng 125:e7399

70. Kaiser T, Raasch T, Remmers JJC, Geers MGD (2023) A wavelet-enhanced adaptive hierarchical FFT-based approach for the efficient solution of microscale boundary value problems. Comput Methods Appl Mech Eng 409:115959

71. Donoho DL (1992) "Interpolating wavelet transforms," Preprint, Department of Statistics, Stanford University, vol. 2, no. 3, pp. 1–54

72. Daubechies I (1988) Orthonormal bases of compactly supported wavelets. Communications in Pure and Applied Mathematics 41:909–996

73. Daubechies I (1992) Ten Lectures On Wavelets. Society for Industrial and Applied Mathematics

74. Cohen A, Dahmen W, DeVore R (2001) Adaptive wavelet methods for elliptic operator equations: convergence rates. Math Comput 70(233):27–75

75. Cohen A, Dahmen W, DeVore R (2001) Adaptive finite element methods with convergence rates. Numer Math 97(2):219–268

76. Mehraeen S, Chen J-S (2005) Wavelet Galerkin method in multi-scale homogenization of heterogeneous media. Int J Numer Meth Eng 66(3):381–403

77. Li B, Chen X (2014) Wavelet-based numerical analysis: A review and classification. Finite Elem Anal Des 81:14–31

78. Vasilyev OV, Bowman C (2000) Second-generation wavelet collocation method for the solution of partial differential equations. J Comput Phys 165(2):660–693

79. Deslauriers G, Dubuc S (1989) "Symmetric iterative interpolation processes," Constructive Approximation: Special Issue: Fractal Approximation, pp. 49–68
80. Burgos RB, Santos MAC, Rosas e Silva R (2013) "Deslauriers-Dubuc interpolating wavelet beam finite element," *Finite Elements in Analysis and Design*, vol. 75, pp. 71–77

81. Dave M (2025) "Wavelet-Enhanced FFT-Based Homogenization," Master's thesis, University of Duisburg-Essen, Essen, Germany

82. Saito N, Beylkin G (1993) Multiresolution Representations Using the Auto-Correlation Functions of Compactly Supported Wavelets. IEEE Trans Signal Process 41(12):3584–3590

83. Fukuda N (2013) On the wavelet-Galerkin method with Deslauriers-Dubuc interpolating scaling functions. Tsukuba Journal of Mathematics 37(2):321–337

84. Leer BV (1977) Towards the ultimate conservative difference scheme III. Upstream-centered finite-difference schemes for ideal compressible flow. J Comput Phys 23(3):263–275

85. Fornberg B (1988) Generation of finite difference formulas on arbitrarily spaced grids. Math Comput 51(184):699–706

86. McCormick K, Wells RO (1994) Wavelet calculus and finite difference operators. Math Comput 63(207):155–173

87. Schneider M (2022) Voxel-based finite elements with hourglass control in FFT-based computational homogenization. Int J Numer Meth Eng 123(24):6286–6313

88. Frigo M, Johnson SG (2005) The Design and Implementation of FFTW3. Proc IEEE 93:216–231

89. Trefethen LN, Weideman JAC (2014) The Exponentially Convergent Trapezoidal Rule. SIAM Rev 56(3):385–458

90. Schneider M (2019) On the Barzilai-Borwein basic scheme in FFT-based computational homogenization. Int J Numer Meth Eng 118(8):482–494

91. Zeman J, Vondřejc J, Novák J, Marek I (2010) Accelerating a FFT-based solver for numerical homogenization of periodic media by conjugate gradients. J Comput Phys 229(21):8065–8071

92. Segurado J, Llorca J, González C (2002) On the accuracy of mean-field approaches to simulate the plastic deformation of composites. Scripta Mater 46(7):525–529

93. Schneider M, Hofmann T, Andrä H, Lechner P, Ettemeyer F, Volk W, Steeb H (2018) Modeling the microstructure and computing effective elastic properties of sand core materials. Int J Solids Struct 143:1–17

94. Ettemeyer F, Lechner P, Hofmann T, Andrä H, Schneider M, Grund D, Volk W, Günther D (2020) Digital Sand Core Physics: Predicting physical properties of sand cores by simulations on digital microstructures. Int J Solids Struct 188–189:155–168

95. Donval E, Schneider M, Grimm-Strele H, Godehardt M, Burger R, Lechner P, Günther D, Andrä H (2025) A directional contraction method to model sand-based binder jet 3D printed materials. Int J Solids Struct 312:113260

96. Schneider M (2017) The Sequential Addition and Migration method to generate representative volume elements for the homogenization of short fiber reinforced plastics. Comput Mech 59:247–263

97. Wu L, Adam L, Doghri I, Noels L (2017) An incremental-secant mean-field homogenization method with second statistical moments for elasto-visco-plastic composite materials. Mech Mater 114:180–200

98. Kabel M, Fliegener S, Schneider M (2016) Mixed boundary conditions for FFT-based homogenization at finite strains. Comput Mech 57(2):193–210

99. Schneider M (2020) A dynamical view of nonlinear conjugate gradient methods with applications to FFT-based computational micromechanics. Comput Mech 66:239–257

100. Jaffard S (1992) Wavelet Methods for Fast Resolution of Elliptic Problems. SIAM J Numer Anal 29(4):965–986

101. Dahmen W, Kunoth A (1992) Multilevel preconditioning. Numer Math 63:315–344

102. Ortiz M, Stainier L (1999) The variational formulation of viscoplastic constitutive updates. Comput Methods Appl Mech Eng 171(3–4):419–444

103. Miehe C (2002) Strain-driven homogenization of inelastic microstructures and composites based on an incremental variational formulation. Int J Numer Meth Eng 55:1285–1322

104. Sanchez-Palencia E (1980) *Non-homogeneous media and vibration theory*. No. 127 in Lecture Notes in Physics, New York: Springer

105. Bakhvalov N, Panasenko G (1989) Homogenisation: Averaging Processes in Periodic Media. Kluwer, Mathematics and its Applications, Dordrecht

106. Bensoussan A, Lion J-L, Papanicolaou G (1978) *Asymptotic Analysis for Periodic Structures*. Studies in Mathematics and its Applications, Amsterdam: North-Holland

107. Nguyen BH, Schneider M (2025) "Symmetries in stochastic homogenization and adjustments for the RVE method," International Journal of Solids and Structures, vol. Accepted, pp. 1–48

108. Milton GW (2002) The Theory of Composites. Cambridge University Press, Cambridge

109. Michel JC, Moulinec H, Suquet P (2001) A computational scheme for linear and non-linear composites with arbitrary phase contrast. Int J Numer Meth Eng 52:139–160

110. Moulinec H, Silva F (2014) Comparison of three accelerated FFT-based schemes for computing the mechanical response of composite materials. Int J Numer Meth Eng 97:960–985

111. Nesterov Y (2004) Introductory Lectures on Convex Optimization: A Basic Course. Kluwer Academic Publishers, Mathematics and its applications, Doordrecht

112. Schneider M, Wicht D, Böhlke T (2019) On polarization-based schemes for the FFT-based computational homogenization of inelastic materials. Comput Mech 64(4):1073–1095

113. Beylkin G (1992) On the representation of operators in bases of compactly supported wavelets. SIAM J Numer Anal 29(6):1716–1740

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.