FFT-based homogenization on periodic anisotropic translation invariant spaces

TL;DR

This paper develops a unified mathematical framework for homogenization in elasticity using anisotropic translation invariant spaces, combining Fourier and finite element methods, and introduces a periodised Green operator with proven properties.

Contribution

It unifies Fourier series and finite element approaches in a common anisotropic translation invariant space framework for homogenization.

Findings

The periodised Green operator is bounded.

The operator is a projection if generated by a Dirichlet kernel.

Numerical examples demonstrate framework flexibility.

Abstract

In this paper we derive a discretisation of the equation of quasi-static elasticity in homogenization in form of a variational formulation and the so-called Lippmann-Schwinger equation, in anisotropic spaces of translates of periodic functions. We unify and extend the truncated Fourier series approach, the constant finite element ansatz and the anisotropic lattice derivation. The resulting formulation of the Lippmann-Schwinger equation in anisotropic translation invariant spaces unifies and analyses for the first time both the Fourier methods and finite element approaches in a common mathematical framework. We further define and characterize the resulting periodised Green operator. This operator coincides in case of a Dirichlet kernel corresponding to a diagonal matrix with the operator derived for the Galerkin projection stemming from the truncated Fourier series approach and to the anisotropic lattice derivation for all other Dirichlet kernels. Additionally, we proof the boundedness of the periodised Green operator. The operator further constitutes a projection if and only if the space of translates is generated by a Dirichlet kernel. Numerical examples for both the de la Vall\'ee Poussin means and Box splines illustrate the flexibility of this framework.

Equations

$$⟨f,g⟩=\frac{1}{(2\pi {)}^d}{\int }_{{\mathbb{T}}^d}f(\mathbf{x})\overline{g(\mathbf{x})}\mathrm{d}\mathbf{x},f,g\in L^2({\mathbb{T}}^d)\text{.}$$

$$f(\mathbf{x})=\sum _{\mathbf{k}\in {\mathbb{Z}}^d}{c}_{\mathbf{k}}(f){\mathrm{e}}^{\mathrm{i}{\mathbf{k}}^{\mathrm{T}}\mathbf{x}},$$

$$⟨\mathbf{c},\mathbf{d}⟩=\sum _{\mathbf{k}\in {\mathbb{Z}}^d}{c}_{\mathbf{k}}\overline{d_{\mathbf{k}}},\mathbf{c},\mathbf{d}\in {\ell }^2({\mathbb{Z}}^d)\text{.}$$

$$⟨f,g⟩=⟨\mathbf{c}(f),\mathbf{c}(g)⟩=\sum _{\mathbf{k}\in {\mathbb{Z}}^d}{c}_{\mathbf{k}}(f)\overline{c_{\mathbf{k}}(g)}\text{.}$$

$$\mathbf{h}\equiv \mathbf{k}\mathrm{mod}\mathbf{M}\iff \exists \mathbf{z}\in {\mathbb{Z}}^d:\mathbf{k}=\mathbf{h}+\mathbf{Mz}\text{.}$$

$$\Lambda (\mathbf{M}):={\mathbf{M}}^{-1}{\mathbb{Z}}^d=\{\mathbf{y}\in {\mathbb{R}}^d:\mathbf{My}\in {\mathbb{Z}}^d\},$$

$$\bigl{[}\mathbf{a}\bigr{]}_{\mathbf{k}}^{\mathbf{M}}\coloneqq\sum_{\mathbf{z}\in\mathbb{Z}^{d}}a_{\mathbf{k}+\mathbf{M}^{\mathrm{T}}\mathbf{z}},\qquad\mathbf{k}\in\mathbb{Z}^{d}.$$

$$\mathcal{F}(\mathbf{M})\coloneqq\frac{1}{\sqrt{m}}\Bigl{(}{\,\mathrm{e}}^{-2\pi{\mathrm{i}}\mathbf{h}^{\mathrm{T}}\mathbf{y}}\Bigr{)}_{\mathbf{h}\in\operatorname{\mathcal{G}}(\mathbf{M}^{\mathrm{T}}),\,\mathbf{y}\in\operatorname{\mathcal{P}}(\mathbf{M})},$$

$$\hat{\mathbf{a}}=({\hat{a}}_{\mathbf{h}}{)}_{\mathbf{h}\in G({\mathbf{M}}^{\mathrm{T}})}=\mathcal{F}(\mathbf{M})\mathbf{a},$$

$$V_{\mathbf{M}}^{f}\coloneqq\operatorname{span}\bigl{\{}\operatorname{\mathcal{T}}(\mathbf{y})f\,:\,\mathbf{y}\in\operatorname{\mathcal{P}}(\mathbf{M})\bigr{\}}$$

$$g=\sum _{\mathbf{y}\in P(\mathbf{M})}{a}_{\mathbf{y}}T(\mathbf{y})f$$

$$c_{\mathbf{h}+{\mathbf{M}}^{\mathrm{T}}\mathbf{z}}(g)={\hat{a}}_{\mathbf{h}}{c}_{\mathbf{h}+{\mathbf{M}}^{\mathrm{T}}\mathbf{z}}(f)\text{for all }\mathbf{h}\in G({\mathbf{M}}^{\mathrm{T}}),\mathbf{z}\in {\mathbb{Z}}^d\text{,}$$

$$\operatorname{\mathcal{T}}_{\mathbf{M}}\coloneqq\Bigl{\{}f\,:\,f=\sum_{\mathbf{h}\in\operatorname{\mathcal{G}}(\mathbf{M}^{\mathrm{T}})}a_{\mathbf{h}}{\,\mathrm{e}}^{{\mathrm{i}}\mathbf{h}^{\mathrm{T}}\circ},\ a_{\mathbf{h}}\in\mathbb{C}\Bigr{\}}\text{,}$$

$$S_{\mathbf{M}}f:=\sum _{\mathbf{h}\in G({\mathbf{M}}^{\mathrm{T}})}{c}_{\mathbf{h}}(f){\mathrm{e}}^{\mathrm{i}{\mathbf{h}}^{\mathrm{T}}\circ }.$$

$$c_{\mathbf{h}}^{\mathbf{M}}(f):=\frac{1}{m}\sum _{\mathbf{y}\in P(\mathbf{M})}f(2\pi \mathbf{y}){\mathrm{e}}^{-2\pi \mathrm{i}{\mathbf{h}}^{\mathrm{T}}\mathbf{y}},\mathbf{h}\in G({\mathbf{M}}^{\mathrm{T}})\text{.}$$

$$c_{\mathbf{h}}^{\mathbf{M}}(f)=\sum_{\mathbf{z}\in\mathbb{Z}^{d}}c_{\mathbf{h}+\mathbf{M}^{\mathrm{T}}\mathbf{z}}(f)=\bigl{[}\mathbf{c}(f)\bigr{]}_{\mathbf{h}}^{\mathbf{M}},\quad\mathbf{h}\in\operatorname{\mathcal{G}}(\mathbf{M)^{\mathrm{T}}}\text{.}$$

$$I_{\mathbf{M}}(2\pi \mathbf{y}):={\delta }_{0,\mathbf{y}}^{{\mathbf{E}}_d},\mathbf{y}\in P(\mathbf{M}),\text{where }{\delta }_{\mathbf{x},\mathbf{y}}^{\mathbf{M}}:=\{\begin{array}{ll}1,& \text{ if }\mathbf{y}\equiv \mathbf{x}\mathrm{mod}\mathbf{M},\\ 0,& \text{ else.}\end{array}$$

$$\bigl{[}\mathbf{c}(f)\bigr{]}_{\mathbf{h}}^{\mathbf{M}}\neq 0\quad\text{ for all }\mathbf{h}\in\operatorname{\mathcal{G}}(\mathbf{M}^{\mathrm{T}}).$$

$$\sum _{\mathbf{z}\in {\mathbb{Z}}^d}|c_{\mathbf{h}+{\mathbf{M}}^{\mathrm{T}}\mathbf{z}}(f){|}^2>0\text{holds for all }\mathbf{h}\in G({\mathbf{M}}^{\mathrm{T}}).$$

$$\sum _{\mathbf{z}\in {\mathbb{Z}}^d}|c_{\mathbf{h}+{\mathbf{M}}^{\mathrm{T}}\mathbf{z}}(f){|}^2=\frac{1}{m}\text{holds for all }\mathbf{h}\in G({\mathbf{M}}^{\mathrm{T}}).$$

$$\tilde{g}(2\pi \mathbf{y})=g(2\pi \mathbf{y}),\mathbf{y}\in P(\mathbf{M}),$$

$$\hat{a}_{\mathbf{h}}=\frac{\bigl{[}\mathbf{c}(\tilde{g})\bigr{]}_{\mathbf{h}}^{\mathbf{M}}}{\bigl{[}\mathbf{c}(f)\bigr{]}_{\mathbf{h}}^{\mathbf{M}}},\qquad\mathbf{h}\in\operatorname{\mathcal{G}}(\mathbf{M}^{\mathrm{T}}),$$

$$\widehat{B}^{c}_{\mathbf{\Xi}}(\mathbf{y})=\prod_{\mathbf{\xi}\in\mathbf{\Xi}}\operatorname{sinc}\bigl{(}\frac{1}{2}\mathbf{\xi}^{\mathrm{T}}\mathbf{y}\bigr{)},\qquad\operatorname{sinc}(t)\coloneqq\frac{\sin(t)}{t},$$

$$g_{\mathrm{p}}(\mathbf{x})\coloneqq\sum_{\mathbf{z}\in\mathbb{Z}^{d}}g\bigl{(}\frac{\mathbf{x}}{2\pi}-\mathbf{z}\bigr{)}.$$

$$c_{\mathbf{k}}(g_{\mathrm{p}})=\hat{g}(2\pi \mathbf{k}),\mathbf{k}\in {\mathbb{Z}}^d$$

$$c_{\mathbf{k}}(B_{\mathbf{\Xi}})\coloneqq\widehat{B}^{c}_{\mathbf{\Xi}}(2\pi\mathbf{k})=\prod_{\xi\in\Xi}\operatorname{sinc}\bigl{(}\pi\xi^{\mathrm{T}}\mathbf{k}\bigr{)},\qquad\mathbf{k}\in\mathbb{Z}^{d}.$$

$$f_{\mathbf{M},\mathbf{\Xi }}(\mathbf{x}):=B_{\mathbf{\Xi }}({\mathbf{M}}^{-1}\mathbf{x}).$$

$$g_{\alpha }(\mathbf{x}):=B_{\mathbf{\Xi }}(\mathbf{x}),\mathbf{\Xi }=\left(\begin{array}{c}\mathrm{diag}(\alpha ){\mathbf{I}}_d\end{array}\right)\in {\mathbb{R}}^{d\times 2d},\alpha \in [0,1{]}^d$$

$$c_{\mathbf{k}}(f_{\mathbf{M},g}):=\frac{1}{\sqrt{m}}g({\mathbf{M}}^{-\mathrm{T}}\mathbf{k}),\mathbf{k}\in {\mathbb{Z}}^d$$

$$g(\mathbf{x})\coloneqq\begin{cases}1,&\mathbf{x}\in\bigl{[}-\frac{1}{2},\frac{1}{2}\bigr{)}^{d}\\ 0,&\text{otherwise}.\end{cases}$$

Full text

Keywords: homogenization, anisotropic lattices, translation invariant spaces, Lippmann-Schwinger equation

MSC 2000: 42B35, 42B37, 65T40, 74B05, 74E30

1 Introduction

Many modern tools and products use composites, i.e. mixtures of two or more materials with distinct elastic properties to obtain certain flexible behaviour, dampening effects, or longevity. Homogenization aims to simplify simulations by replacing the microscopically composed material by a homogeneous one which behaves the same on the macroscopic scale. Mathematically one assumes is a periodic microstructure, i.e. a structure that can be represented by a certain unit cell with periodic boundary conditions.

For the simulation of such elastic composite structures Moulinec and Suquet [18, 18] derive an algorithm based on the fast Fourier transform. This algorithm, called the Basic Scheme, inspired many similar numerical approaches based on using discretised differential operators [28, 23, 22] and extensions to porous media [15, 22]. Information on sub-structures of the geometry is incorporated into the solution method in [12]. The Basic Scheme is generalized to problems of higher order, i.e. derivatives of strain and stiffness [26], and the solution of the arising linear system by Krylov subspace methods is analysed in [29].

Vondřejc et.al. [27] show that the method of Moulinec and Suquet can also be understood as a Galerkin projection using truncated Fourier series. This idea is generalized in [3] to anisotropic lattices thus allowing to take directional information on the geometrical structure or the orientation of interfaces between materials into account. Brisard and Dormieux [8, 9] use constant finite elements to arrive at the Basic Scheme with a modified linear operator, based on an energy based formulation.

In this paper we unify and extend the approaches of Vondřejc et.al., Brisard and Dormieux, and the anisotropic lattice ansatz obtaining a discretisation of the equations for quasi-static elasticity in homogenization in anisotropic spaces of periodic translates. Vondřejc et.al. show that a variational equation and a formulation with the strain as a fixed-point, the so-called Lippmann-Schwinger equation, are equivalent by means of a projection operator derived from the Green operator. This paper introduces a periodised Green operator on the space of translates for the Lippmann-Schwinger equation. Furthermore we classify the properties of this operator in case of spaces of translates and prove that it induces a projection operator if and only if the space of translates is generated by the Dirichlet kernel. Hence the (anisotropic) truncated Fourier series emerges as a case with special properties of the general setting introduced here. This introduces further insight into the equivalence of the variational and the Lippmann-Schwinger equation formulation of Vondřejc et.al. The mathematical framework this paper introduces unifies and analyses the approaches of Fourier methods and finite elements for the first time. Especially using translates of Box splines [6] as ansatz functions incorporates the constant finite elements of Brisard and Dormieux allowing also for anisotropic finite elements of arbitrary smoothness. A different approach to solve the equation using linear finite elements with full quadrature is shown in [22] and is based on replacing the continuous differential operator by a discrete one.

Spaces of translates can for example be generated by de la Vallée Poussin means which provide a generalization of the Dirichlet and Fejér kernel. They combine a finite support in frequency domain with good localization in space [11]. These functions introduce a trade off between damping of the Gibbs phenomenon and reproduction of multivariate trigonometric monomials. They allow for better predictions of the elastic macroscopic properties of the composite material and result in smoother — and thus better — solutions. Further, different Box splines and their influence on the solution are demonstrated.

The remainder of the paper is structured as follows. After reviewing important properties of anisotropic spaces of translates in Section 2 the partial differential equation of quasi-static elasticity in homogenization is introduced in Section 3. Then, the periodised Green operator on spaces of translates is introduced. This operator is subsequently analysed regarding its equivalence to a projection operator and then used to discretise the Lippmann-Schwinger equation. Numerical examples are then provided in Section 4 and a conclusion is drawn in Section 5.

2 Preliminaries

Throughout this paper we will employ the following notation: the symbols $a\in\mathbb{C}$, $\mathbf{a}\in\mathbb{C}^{d}$ and $\mathbf{A}\in\mathbb{C}^{d\times d}$ denote scalars, vectors, and matrices, respectively. The only exception from this are $f,g,h$ which are reserved for functions. We denote the inner product of two vectors by $\mathbf{a}^{\mathrm{T}}\mathbf{b}\coloneqq\sum_{i}a_{i}b_{i}$ and reserve the symbol $\langle\cdot,\cdot\rangle$ for inner products of two functions or two generalized sequences, respectively. For a complex number $a=b+{\mathrm{i}}c$, $b,c\in\mathbb{R}$, we denote the complex conjugate by $\overline{a}\coloneqq b-{\mathrm{i}}c$. Constants like Euler’s number ${\,\mathrm{e}}$ or the imaginary unit ${\mathrm{i}}$, i.e. ${\mathrm{i}}^{2}=-1$, are set upright.

Usually, we are concerned with $d$-dimensional data, where $d=2,3$, but the theory is written in arbitrary dimensions. Sets are denoted by capital case calligraphic letters, e.g. $\mathcal{P}$ or $\mathcal{G}$ and the same for the Fourier transform $\mathcal{F}$ which all might depend on parameters given in round brackets. We denote second-order tensors by small Greek letters as $\lambda,\varepsilon$ with entries $\lambda_{ij}$ are indexed again by scalars $i,j$ and similarly we denote fourth-order tensors by capital calligraphic letters, where $\mathcal{C}$ is the most prominent one.

2.1 Arbitrary patterns and the Fourier transform

The space of functions we are concerned with is the Hilbert space $L^{2}(\mathbb{T}^{d})$ of (equivalence classes of) square integrable functions on the $d$-dimensional torus $\mathbb{T}\cong[-\pi,\pi)^{d}$ with inner product

[TABLE]

In several cases, the functions of interest are tensor-valued. For these functions, we take the tensor product of the Hilbert space, e.g. $L^{2}(\mathbb{T}^{d})^{n\times n}$ for the space of functions $f\colon\mathbb{T}^{d}\to\mathbb{C}^{n\times n}$ that have values being $n\times n$-dimensional matrices. The following preliminaries can be generalized to these tensor product spaces by performing the operations element wise. We restrict the following of this subsection therefore to the case of $L^{2}(\mathbb{T}^{d})$.

Every function $f\in L^{2}(\mathbb{T}^{d})$ can be written in its Fourier series representation

[TABLE]

introducing the multivariate Fourier coefficients $c_{\mathbf{k}}(f)=\langle f,{\,\mathrm{e}}^{{\mathrm{i}}\mathbf{k}^{\mathrm{T}}\circ}\rangle$, $\mathbf{k}\in\mathbb{Z}^{d}$. The equality in (1) is meant in $L^{2}(\mathbb{T}^{d})$ sense. We denote by \mathbf{c}(f)=\bigl{\{}c_{\mathbf{k}}(f)\bigr{\}}_{\mathbf{k}\in\mathbb{Z}^{d}}\in\ell^{2}(\mathbb{Z}^{d}) generalized sequences which form a Hilbert space with the inner product

[TABLE]

The Parseval equation reads

[TABLE]

The pattern and the generating set.

For any regular matrix $\mathbf{M}\in\mathbb{Z}^{d\times d}$ we define the congruence relation for $\mathbf{h},\mathbf{k}\in\mathbb{Z}^{d}$ with respect to $\mathbf{M}$ by

[TABLE]

We define the lattice

[TABLE]

and the pattern $\operatorname{\mathcal{P}}(\mathbf{M})$ as any set of congruence representants of the lattice with respect to $\bmod\ \mathbf{1}$, e.g. $\Lambda(\mathbf{M})\cap[0,1)^{d}$ or \Lambda(\mathbf{M})\cap\bigl{[}-\tfrac{1}{2},\tfrac{1}{2}\bigr{)}^{d}. For the rest of the paper we will refer to the set of congruence class representants in the symmetric unit cube \bigl{[}-\tfrac{1}{2},\tfrac{1}{2}\bigr{)}^{d}. The generating set $\operatorname{\mathcal{G}}(\mathbf{M})$ is defined by $\operatorname{\mathcal{G}}(\mathbf{M})\coloneqq\mathbf{M}\operatorname{\mathcal{P}}(\mathbf{M})$ for any pattern $\operatorname{\mathcal{P}}(\mathbf{M})$. For both, the number of elements is given by $\lvert\operatorname{\mathcal{P}}(\mathbf{M})\rvert=\lvert\operatorname{\mathcal{G}}(\mathbf{M})\rvert=\lvert\det{\mathbf{M}}\rvert\eqqcolon m,$ which follows directly from [6, Lemma II.7].

For a regular integer matrix $\mathbf{M}\in\mathbb{Z}^{d\times d}$ and an absolutely summable generalized sequence $\mathbf{a}=\{a_{\mathbf{k}}\}_{\mathbf{k}\in\mathbb{Z}^{d}}$ we further define the bracket sum

[TABLE]

The bracket sum is periodic with respect to $\mathbf{M}^{\mathrm{T}}$, i.e., \bigl{[}\mathbf{a}\bigr{]}_{\mathbf{k}}^{\mathbf{M}}=\bigl{[}\mathbf{a}\bigr{]}_{\mathbf{k}+\mathbf{M}^{\mathrm{T}}\mathbf{z}}^{\mathbf{M}} holds for any $\mathbf{k},\mathbf{z}\in\mathbb{Z}^{d}$.

A fast Fourier transform on patterns.

The discrete Fourier transform on the pattern $\operatorname{\mathcal{P}}(\mathbf{M})$ is defined [10] by

[TABLE]

where $\mathbf{h}\in\operatorname{\mathcal{G}}(\mathbf{M}^{\mathrm{T}})$ indicate the rows and $\mathbf{y}\in\operatorname{\mathcal{P}}(\mathbf{M})$ indicate the columns of the Fourier matrix $\mathcal{F}(\mathbf{M})$. The discrete Fourier transform on $\operatorname{\mathcal{P}}(\mathbf{M})$ is defined for a vector $\mathbf{a}=(a_{\mathbf{y}})_{\mathbf{y}\in\operatorname{\mathcal{P}}(\mathbf{M})}\in\mathbb{C}^{m}$ arranged in the same ordering as the columns in (4) by

[TABLE]

where the resulting vector $\mathbf{\hat{a}}$ is ordered as the columns of $\mathcal{F}(\mathbf{M})$ in (4). Its implementation yields complexity of $\mathcal{O}(m\log m)$ similar to the classical Fourier transform, when the ordering is fixed as described in [1, Theorem 2]. Note that for the so-called rank-1-lattices, the Fourier transform on the pattern even reduces to a one-dimensional FFT for patterns in arbitrary dimensions [13].

2.2 Translation invariant spaces of periodic functions

Spaces of translates and interpolation.

A space of functions $V\subset L^{2}(\mathbb{T}^{d})$ is called $\mathbf{M}$-invariant, if for all $\mathbf{y}\in\operatorname{\mathcal{P}}(\mathbf{M})$ and all functions $f\in V$ the translates $\operatorname{\mathcal{T}}({\mathbf{y}})f\coloneqq f(\cdot-2\pi\mathbf{y})\in V$. Especially the space

[TABLE]

of translates of $f$ is $\mathbf{M}$-invariant. A function $g\in V_{\mathbf{M}}^{f}$ is of the form

[TABLE]

For $f\in L^{2}(\mathbb{T}^{d})$ an easy calculation on the Fourier coefficients using the unique decomposition of $\mathbf{k}\in\mathbb{Z}^{d}$ into $\mathbf{k}=\mathbf{h}+\mathbf{M}^{\mathrm{T}}\mathbf{z}$, $\mathbf{h}\in\operatorname{\mathcal{G}}(\mathbf{M}^{\mathrm{T}}),\mathbf{z}\in\mathbb{Z}^{d}$, yields, that $g\in V_{\mathbf{M}}^{f}$ holds if and only if [14, Theorem 3.3]

[TABLE]

holds, where \mathbf{\hat{a}}=\bigl{(}\hat{a}_{\mathbf{h}}\bigr{)}_{\mathbf{h}\in\operatorname{\mathcal{G}}(\mathbf{M}^{\mathrm{T}})}=\sqrt{m}\mathcal{F}(\mathbf{M})\mathbf{a} denotes the discrete Fourier transform of \mathbf{a}=\bigl{(}a_{\mathbf{y}}\bigr{)}_{\mathbf{y}\in\operatorname{\mathcal{P}}(\mathbf{M})}\in\mathbb{C}^{m}, see [14]. Using the space of trigonometric polynomials on the generating set $\operatorname{\mathcal{G}}(\mathbf{M}^{\mathrm{T}})$, which is denoted by

[TABLE]

we define for a function $f\in L^{2}(\mathbb{T}^{d})$ the Fourier partial sum $\operatorname{S}_{\mathbf{M}}f\in\operatorname{\mathcal{T}}_{\mathbf{M}}$ by

[TABLE]

The discrete Fourier coefficients $c_{\mathbf{k}}^{\mathbf{M}}(f)$ of a function $f$ that is evaluated pointwise on the pattern $\operatorname{\mathcal{P}}(\mathbf{M})$ are defined by

[TABLE]

The discrete Fourier coefficients $c_{\mathbf{h}}^{\mathbf{M}}$ are related to the Fourier coefficients for a function $f\in A(\mathbb{T}^{d})$, where $A(\mathbb{T}^{d})$ denotes the Wiener Algebra, i.e. the space of functions with absolutely convergent Fourier series. This relation is given in the following Lemma, also known as the aliasing formula, see e.g. [4, Lemma 2].

Lemma 2.1**.**

Let $f\in A(\mathbb{T}^{d})$ and the regular matrix $\mathbf{M}\in\mathbb{Z}^{d\times d}$ be given. Then the discrete Fourier coefficients $c_{\mathbf{h}}^{\mathbf{M}}(f)$ are given by

[TABLE]

When looking at the space $V_{\mathbf{M}}^{f}$ of translates, the following definition is crucial in order to approximate a function $g$ by using these translates.

Definition 2.2**.**

Let $\mathbf{M}\in\mathbb{Z}^{d\times d}$ be a regular matrix. A function $\operatorname{I}_{\mathbf{M}}\in V_{\mathbf{M}}^{\varphi}$ is called fundamental interpolant or Lagrange function of $V_{\mathbf{M}}^{\varphi}$ if

[TABLE]

The following lemma characterizes the existence of such a fundamental interpolant in a space $V_{\mathbf{M}}^{f}$ of translates and collects some properties of the translates themselves, see [2, Lemma 1.23] and [4, Lemma 2].

Lemma 2.3**.**

Given a regular matrix $\mathbf{M}\in\mathbb{Z}^{d\times d}$ and a function $f\in A(\mathbb{T}^{d})$, then the following holds.

1. a)

The fundamental interpolant $\operatorname{I}_{\mathbf{M}}\in V_{\mathbf{M}}^{f}$ exists if and only if

[TABLE]

If the fundamental interpolant $\operatorname{I}_{\mathbf{M}}\in V_{\mathbf{M}}^{f}$ exists, it is uniquely determined. 2. b)

The set of translates \bigl{\{}\operatorname{\mathcal{T}}(\mathbf{y})f:\ \mathbf{y}\in\operatorname{\mathcal{P}}(\mathbf{M})\bigr{\}} is linear independent if and only if

[TABLE] 3. c)

The set of translates \bigl{\{}\operatorname{\mathcal{T}}(\mathbf{y})f:\ \mathbf{y}\in\operatorname{\mathcal{P}}(\mathbf{M})\bigr{\}} is an orthonormal basis of $V_{\mathbf{M}}^{f}$ if and only if

[TABLE] 4. d)

Given a function $\tilde{g}\in A(\mathbb{T}^{d})$ we can obtain a function $g\in V_{\mathbf{M}}^{f}$ fulfilling

[TABLE]

provided that the fundamental interpolant exists (which also implies linear independence of the translates on $f$) as

[TABLE]

where the coefficients $\hat{a}_{\mathbf{h}}$ yield $g$ in Fourier coefficients by (6).

By using Lemma 2.3 changing from sampling values, i.e. the coefficients on the pattern $\operatorname{\mathcal{P}}(\mathbf{M})$ of the fundamental interpolant, to coefficients with respect to $f$ in the corresponding space $V_{\mathbf{M}}^{f}$ of translates can be done by using the Fourier transform (5) and the Fourier coefficients $c_{\mathbf{k}}(f)$ of $f$.

For the remainder of this paper, two special spaces of translates are of interest, the periodised Box splines and the de la Vallée Poussin means.

Periodised Box splines.

Let $\mathbf{\Xi}=(\mathbf{\xi}_{1},\ldots,\mathbf{\xi}_{s})\in\mathbb{R}^{d\times s}$ denote a set of column vectors $\mathbf{\xi}_{i}$, $i=1,\ldots,s$, where we assume that these vectors span the $\mathbb{R}^{d}$, i.e. especially we have $s\geq d$. Then the centred Box spline $B^{c}_{\mathbf{\Xi}}$ can be defined via its Fourier transform as

[TABLE]

cf. [6, p. 11]. A Box spline has compact support. For a function $g\colon\mathbb{R}^{d}\to\mathbb{C}$ we can introduce its periodisation

[TABLE]

Its Fourier coefficients can be directly computed from the continuous Fourier transform $\hat{g}$ of $g$, cf. e.g. [4, p. 41], as

[TABLE]

We combine these two to introduce the periodised Box Spline $B_{\mathbf{\Xi}}\colon\mathbb{T}^{d}\to\mathbb{R}$ via its Fourier coefficients as

[TABLE]

Finally, we obtain by scaling the periodised pattern Box Spline $f_{\mathbf{M},\mathbf{\Xi}}$

[TABLE]

Note that its translates might not be linearly independent for an arbitrary set of vectors in $\mathbf{\Xi}$, see also [20]. However, by [4], see also [7, Sect. 4], the matrices $\mathbf{\Xi}\in\mathbb{R}^{d\times(p+q+r)}$ of the form $\xi_{1}=\ldots,\xi_{p}=(1,0)^{\mathrm{T}}$, $\xi_{p+1}=\ldots=\xi_{p+q}=(0,1)^{\mathrm{T}}$, and $\xi_{p+q+1}=\ldots=\xi_{p+q+r}=(1,1)^{\mathrm{T}}$, where at least two of the values $p,q,r$ are larger than [math], induce a periodised pattern Box spline $f_{\mathbf{M},\Xi}$ with linear independent translates.

De la Vallée Poussin means.

A special case of $\mathbf{M}$-invariant spaces are the ones defined via de la Vallée Poussin means, following the construction of [5]. We call a function $g\colon\mathbb{R}^{d}\to\mathbb{R}$ admissible if the function fulfils

1. a)

$g(\mathbf{x})\geq 0$ for all $\mathbf{x}\in\mathbb{R}^{d}$, 2. b)

$g(\mathbf{x})>0$ for \mathbf{x}\in\bigl{[}-\tfrac{1}{2},\frac{1}{2})^{d}, 3. c)

$\displaystyle\sum_{\mathbf{z}\in\mathbb{Z}^{d}}g(\mathbf{x}+\mathbf{z})=1$ for all $\mathbf{x}\in\mathbb{R}^{d}$.

This can be for example the Box splines of the form

[TABLE]

where $\mathbf{I}_{d}$ is the $d$-dimensional unit matrix. In the following we define the de la Vallée Poussin means as follows, which is a special case of [5, Definition 4.2] setting $l=0$ therein.

Definition 2.4**.**

Let $\mathbf{M}\in\mathbb{Z}^{d\times d}$ be a regular matrix and $g$ be an admissible function. The function $f_{\mathbf{M},g}$, which is defined by their Fourier coefficients as

[TABLE]

is called de la Vallée Poussin mean.

In case of Box splines $g_{\mathbf{\alpha}}$ this generalizes the one-dimensional de la Vallée Poussin means from [24, 21] to arbitrary patterns including the tensor product case for diagonal matrices $\mathbf{M}$, which where for example used in [25]. We will use the short hand notation $f_{\mathbf{M},\mathbf{\alpha}}\coloneqq f_{\mathbf{M},g_{\mathbf{\alpha}}}$ and omit $\alpha$ whenever its clear from the context. It is easy to see, that by admissibility of $g$ the fundamental interpolant exists for any de la Vallée Poussin mean $f_{\mathbf{M},g}$. The functions $f_{\mathbf{M},\alpha}$ generalize the classical de la Vallée Poussin means to higher dimensions and anisotropic patterns, for which examples are shown in Figure 1 and explained in the following.

Finally, the Dirichlet kernel $D_{\mathbf{M}}$ is defined by using the function $g(\mathbf{x})$ with

[TABLE]

This kernel is comprised in the definition of the generalized de la Vallée Poussin mean as well. Furthermore we obtain the modified Dirichlet Kernel $f_{\mathbf{M},\mathbf{0}}$ as a limiting case of the de la Vallée Poussin case.

As an example we choose \mathbf{M}_{1}=\bigl{(}\begin{smallmatrix}8&0\\ 0&8\end{smallmatrix}\bigr{)} and \mathbf{M}_{2}=\bigl{(}\begin{smallmatrix}4&-2\\ 4&14\end{smallmatrix}\bigr{)}. For $\mathbf{M}_{1}$ we obtain the usual rectangular (pixel grid) pattern $\operatorname{\mathcal{P}}(\mathbf{M}_{1})$ while $\operatorname{\mathcal{P}}(\mathbf{M}_{2})$ models a certain anisotropy, cf. [3, Fig. 2.1]. By further setting $\alpha=\frac{1}{10}\begin{pmatrix}1&1\end{pmatrix}^{\mathrm{T}}$ we obtain the de la Vallée Poussin means $f_{\mathbf{M}_{1},\alpha}$ and $f_{\mathbf{M}_{2},\alpha}$. Their Fourier coefficients $c_{\mathbf{k}}(f_{\mathbf{M}_{1},\alpha})$ and $c_{\mathbf{k}}(f_{\mathbf{M}_{2},\alpha})$ after orthonormalising the translates, cf. Lemma 2.3 c), are shown in Figs. 1 1 (a) and 1 (b), respectively. Note that the first results in $64$ translates, while the second determinant is smaller and results in $58$ translates. The functions in time domain are plotted in Figs. 1 1 (c) and 1 (d), respectively. While the first can also be obtained by a tensor product of one-dimensional de la Vallée Poussin means, cf., e.g. [24], the second one prefers in time domain certain directions due to its anisotropic form.

3 Homogenization on spaces of translates

In the following steps we use anisotropic spaces of translates to discretise the quasi-static equation of linear elasticity in homogenization. First we introduce the necessary spaces and differential operators. With these we can state the partial differential equation we are interested in and two equivalent formulations, a variational equation and the so-called Lippmann-Schwinger equation. These formulations make use of the Green operator $\Gamma^{0}$. Based on this operator we introduce the periodised Green operator $\Gamma^{\mathrm{p}}$ and subsequently analyse its properties and special cases. Next we use this operator to discretise the partial differential equation while splitting the derivation into two steps.

3.1 The elasticity problem in periodic homogenization

FFT-based methods for the equations of linear elasticity in homogenization based on the Lippmann-Schwinger equation are first introduced by Moulinec and Suquet [18, 19]. Based on their method, Vondřejc et. al. [27] interpret the resulting discretisation as a Galerkin projection using trigonometric sums.

In the following we generalize the interpretation using trigonometric sums to spaces of translates of a periodic function. Therein the trigonometric sums will appear as a special case, namely when choosing the Dirichlet kernel’s translates and a diagonal matrix $\mathbf{M}$.

Let $S$ be an arbitrary set, then we introduce the notations

[TABLE]

The space $\operatorname{Sym}_{d}(S)$ corresponds to symmetric matrices built of elements of $S$ and $\operatorname{SSym}_{d}(S)$ corresponds to fourth-order tensors \mathcal{C}=\bigl{(}\mathcal{C}_{ijkl}\bigr{)}_{ijkl}\in S^{d\times d\times d\times d} with minor and major symmetries, i.e. $\mathcal{C}_{ijkl}=\mathcal{C}_{jikl}=\mathcal{C}_{ijlk}=\mathcal{C}_{klij}$.

We endow the space of symmetric matrices \operatorname{Sym}_{d}\bigl{(}\mathbb{R}\bigr{)} with the Frobenius inner product

[TABLE]

We call a function \mathcal{A}\in\operatorname{SSym}_{d}\bigl{(}L^{2}\bigl{(}\mathbb{T}^{d}\bigr{)}\bigr{)} uniformly elliptic if there exist constants $0<l_{\mathcal{A}^{0}}\leq u_{\mathcal{A}^{0}}<\infty$ such that for almost all functions \gamma\in\operatorname{Sym}_{d}\bigl{(}L^{2}\bigl{(}\mathbb{T}^{d}\bigr{)}\bigr{)} it holds true that

[TABLE]

A uniformly elliptic and constant function \mathcal{A}^{0}\in\operatorname{SSym}_{d}\bigl{(}L^{2}\bigl{(}\mathbb{T}^{d}\bigr{)}\bigr{)} can be identified with an element from \operatorname{SSym}_{d}\bigl{(}\mathbb{R}\bigr{)} and is called elliptic in the following.

The Sobolev space $H^{1}(\mathbb{T}^{d})$ is defined via Fourier series, i.e.

[TABLE]

To simplify the equations of linear elasticity and the differential operators occurring therein we additionally introduce for $\mathbf{u}\in H^{1}(\mathbb{T}^{d})^{d}$ the symmetrized gradient operator

[TABLE]

and the divergence operator $\operatorname{div}(\mathbf{u})$ as the formal $L^{2}$-adjoint of $\operatorname{\operatorname{\nabla}_{\mathrm{Sym}}}(\mathbf{u})$. The action of the symmetrized gradient operator in Fourier space is for $\mathbf{k}\in\mathbb{Z}^{d}$ given by

[TABLE]

The basic solution space for the PDE we want to analyse is given by symmetric gradient fields with zero mean

[TABLE]

With these preparations we can now state the partial differential equation describing quasi-static linear elasticity in homogenization.

Definition 3.1**.**

Let \mathcal{C}\in\operatorname{SSym}_{d}\bigl{(}L^{\infty}\bigl{(}\mathbb{T}^{d}\bigr{)}\bigr{)} be uniformly elliptic and let \varepsilon^{0}\in\operatorname{Sym}_{d}\bigl{(}\mathbb{R}\bigr{)} be given. We want to find the strain \varepsilon\in\mathcal{E}\bigl{(}\mathbb{T}^{d}\bigr{)} such that

[TABLE]

holds for all \tilde{\gamma}\in\mathcal{E}\bigl{(}\mathbb{T}^{d}\bigr{)}. With $\mathcal{C}:\varepsilon$ we hereby denote the product of the fourth-order stiffness tensor $\mathcal{C}$ and the second-order strain $\varepsilon$, the symmetry making the order of multiplication irrelevant.

The stiffness distribution $\mathcal{C}$ describes the material behaviour and is in practical applications usually a piece-wise constant function. The role of the macroscopic strain $\varepsilon^{0}$ is that of an overall strain that is applied to the composed material and corresponds to pulling (or compressing) the composite in a certain direction.

For applications one is interested either in the strain field $\varepsilon$ —or derivates like stress or displacements, respectively— or in the so-called effective stiffness of the medium. The overall stiffness behaviour of the composite given by the action of \mathcal{C}^{\mathrm{eff}}\in\operatorname{SSym}_{d}\bigl{(}\mathbb{R}\bigr{)} on \varepsilon^{0}\in\operatorname{Sym}_{d}\bigl{(}\mathbb{R}\bigr{)} is defined by

[TABLE]

where $\varepsilon$ is the solution of (8) corresponding to $\varepsilon^{0}$.

The space of suitable ansatz functions \mathcal{E}\bigl{(}\mathbb{T}^{d}\bigr{)} is in practice difficult to work with, especially regarding the discretisation steps that follow. A method to deal with that is introduced by Vondřejc et. al. [27]. They derive a projection operator $\Gamma^{0}\mathcal{C}^{0}$ with $L^{2}$-adjoint $\mathcal{C}^{0}\Gamma^{0}$ that maps \operatorname{Sym}_{d}\bigl{(}L^{2}\bigl{(}\mathbb{T}^{d}\bigr{)}\bigr{)} onto \mathcal{E}\bigl{(}\mathbb{T}^{d}\bigr{)} and thus replaces the structurally complicated space by a simpler one necessitating a more involved PDE. The operator $\Gamma^{0}$ acts as a second order derivative of a preconditioner that solves the constant coefficient PDE $\operatorname{div}\mathcal{C}^{0}\operatorname{\operatorname{\nabla}_{\mathrm{Sym}}}\mathbf{u}=\mathbf{f}$. For a derivation of the formula see for example [18, 19, 22].

Definition 3.2**.**

For a constant elliptic reference stiffness \mathcal{C}^{0}\in\operatorname{SSym}_{d}\bigl{(}\mathbb{R}\bigr{)} define the Green operator \Gamma^{0}\colon\operatorname{Sym}_{d}\bigl{(}L^{2}\bigl{(}\mathbb{T}^{d}\bigr{)}\bigr{)}\rightarrow\mathcal{E}(\mathbb{T}^{d}) which acts as a Fourier multiplier. Its action on a field \varepsilon\in\operatorname{Sym}_{d}\bigl{(}L^{2}\bigl{(}\mathbb{T}^{d}\bigr{)}\bigr{)} is given by the equation

[TABLE]

This operator allows to reformulate (8) with test functions in \operatorname{Sym}_{d}\bigl{(}L^{2}\bigl{(}\mathbb{T}^{d}\bigr{)}\bigr{)} by projecting them onto the required function space. Further, this additional operator is easy do apply as it acts as a convolution operator, i.e. a Fourier multiplier.

Proposition 3.3**.**

Let \varepsilon^{0}\in\operatorname{Sym}_{d}\bigl{(}\mathbb{R}\bigr{)}, then \varepsilon\in\mathcal{E}\bigl{(}\mathbb{T}^{d}\bigr{)} fulfils

[TABLE]

for all \tilde{\gamma}\in\mathcal{E}\bigl{(}\mathbb{T}^{d}\bigr{)} if and only if

[TABLE]

holds true for all \gamma\in\operatorname{Sym}_{d}\bigl{(}L^{2}\bigl{(}\mathbb{T}^{d}\bigr{)}\bigr{)}.

Proof.

For a proof see [27, Proposition 3]. ∎

By the properties of the operator $\Gamma^{0}\mathcal{C}^{0}$ this equation is equivalent to the so-called Lippmann-Schwinger equation with fixed point $\varepsilon$, cf.[27, Lemma 2]:

Proposition 3.4**.**

Let \varepsilon^{0}\in\operatorname{Sym}_{d}\bigl{(}\mathbb{R}\bigr{)}, then \varepsilon\in\mathcal{E}\bigl{(}\mathbb{T}^{d}\bigr{)} fulfils

[TABLE]

for all \tilde{\gamma}\in\mathcal{E}\bigl{(}\mathbb{T}^{d}\bigr{)} if and only if

[TABLE]

is fulfilled for all \gamma\in\operatorname{Sym}_{d}\bigl{(}L^{2}\bigl{(}\mathbb{T}^{d}\bigr{)}\bigr{)}.

Remark 3.5**.**

The proof makes use of the adjoint operator $\mathcal{C}^{0}\Gamma^{0}$ and the identity

[TABLE]

which holds in weak sense. The above identity is equivalent to $\Gamma^{0}\mathcal{C}^{0}$ being a projection operator that maps constants to zero. These properties are proven in [27, Lemma 2].

With this at hand we can now proceed to discretising the Lippmann-Schwinger equation in a space of translates $V_{\mathbf{M}}^{f}$, $f\in A(\mathbb{T}^{d})$.

3.2 The periodised Green operator

The Lippmann-Schwinger equation is first discretised by Moulinec and Suquet [18, 19] using a Fourier collocation scheme, cf. [29]. The resulting fixed point algorithm inspired many publications, for an overview see for example [17].

In contrast, Vondřejc et. al. [27] employ a Galerkin projection of (10) using trigonometric polynomials as ansatz functions and obtain the same discretisation. In the following we want to generalize this approach to spaces of translates on anisotropic lattices.

Throughout this section $\mathbf{M}\in\mathbb{Z}^{d\times d}$ denotes a regular pattern matrix that defines a translation invariant space $V_{\mathbf{M}}^{f}$ spanned by the translates $\operatorname{\mathcal{T}}(\mathbf{y})f$, $\mathbf{y}\in\operatorname{\mathcal{P}}(\mathbf{M})$, such that a fundamental interpolant exists, see Lemma 2.3 a). Especially, this implies that the translates $\operatorname{\mathcal{T}}(\mathbf{y})f$ of $f$ are linearly independent.

This fundamental interpolant is denoted by $\operatorname{I}_{\mathbf{M}}\in V_{\mathbf{M}}^{f}$, i.e. there exist coefficients $\hat{a}_{\mathbf{h}}$, $\mathbf{h}\in\operatorname{\mathcal{G}}(\mathbf{M}^{\mathrm{T}})$ such that

[TABLE]

holds true for all $\mathbf{h}\in\operatorname{\mathcal{G}}(\mathbf{M}^{\mathrm{T}})$ and $\mathbf{z}\in\mathbb{Z}^{d}$.

From now on we assume for functions \gamma\in\operatorname{Sym}_{d}\bigl{(}V^{f}_{\mathbf{M}}\bigr{)} that

[TABLE]

holds true. We further denote the discrete Fourier transform of $\mathbf{G}=(G_{\mathbf{y}})_{\mathbf{y}\in\operatorname{\mathcal{P}}(\mathbf{M})}$ by $\mathbf{\hat{G}}=(\hat{G}_{\mathbf{h}})_{\mathbf{h}\in\operatorname{\mathcal{G}}(\mathbf{M}^{\mathrm{T}})}=\operatorname{\mathcal{F}}(\mathbf{M})\mathbf{G}$

A Galerkin projection of (10) onto the space of translates requires the definition of a Green operator similar to Definition 3.2. To account for the finite dimensional space $V_{\mathbf{M}}^{f}$ the operator $\Gamma^{0}$ has to be periodised in frequency domain.

Definition 3.6**.**

We call the Fourier multiplier $\Gamma^{\mathrm{p}}$ the periodised Green operator on $V^{f}_{\mathbf{M}}$ and define its action onto a field \gamma\in\operatorname{Sym}_{d}\bigl{(}V^{f}_{\mathbf{M}}\bigr{)} by

[TABLE]

Properties of the periodised Green operator.

In the trigonometric collocation case of Vondřejc et.al. [27] the Green operator keeps the same form after discretisation, i.e. a restriction of the Fourier series to a bounded cube. This also holds true for the generalization to anisotropic patterns [3] where the cube is replaced by a parallelotope, i.e. to the set $\operatorname{\mathcal{G}}(\mathbf{M}^{\mathrm{T}})$. The properties of the Green operator and the projection operator $\Gamma^{0}\mathcal{C}^{0}$ are shown via properties of its Fourier coefficients. Hence, the proofs in the continuous and the discretised case can be done analogously. This is no longer the case for the approach using translation invariant spaces.

Theorem 3.7**.**

The operator $\Gamma^{\mathrm{p}}\mathcal{C}^{0}$ has the $L^{2}$-adjoint $\mathcal{C}^{0}\Gamma^{p}$.

Proof.

The proof for $\Gamma^{0}\mathcal{C}^{0}$ in [27, Lemma 2 (ii)] relies purely on the symmetry of the operator. This symmetry is preserved by the periodised Green operator $\Gamma^{\mathrm{p}}\mathcal{C}^{0}$ and therefore the proof is analogous. ∎

Vondřejc et. al. introduce the operator $\Gamma^{0}\mathcal{C}^{0}$ to project functions in \operatorname{Sym}_{d}\bigl{(}L^{2}\bigl{(}\mathbb{T}^{d}\bigr{)}\bigr{)} onto \mathcal{E}\bigl{(}\mathbb{T}^{d}\bigr{)}. The operator $\Gamma^{\mathrm{p}}\mathcal{C}^{0}$ has similar properties and maps onto the respective discretised versions of these spaces.

Theorem 3.8**.**

*For all \gamma\in\operatorname{Sym}_{d}\bigl{(}V_{\mathbf{M}}^{f}\bigr{)} is holds true that \Gamma^{\mathrm{p}}\mathcal{C}^{0}:\gamma\in\mathcal{E}\bigl{(}\mathbb{T}^{d}\bigr{)}\cap\operatorname{Sym}_{d}\bigl{(}V_{\mathbf{M}}^{f}\bigr{)}. *

Proof.

For $\mathbf{y}\in\operatorname{\mathcal{P}}(\mathbf{M})$ we have that

[TABLE]

First, observe this can be rewritten with

[TABLE]

for $\mathbf{h}\in\operatorname{\mathcal{G}}(\mathbf{M}^{\mathrm{T}})$ as

[TABLE]

Using (6) the result is a function in \operatorname{Sym}_{d}\bigl{(}V_{\mathbf{M}}^{f}\bigr{)}.

Expanding the bracket sum and the Green operator $\Gamma^{0}$ and yields

[TABLE]

We define new Fourier coefficients

[TABLE]

With the decomposition $\mathbf{k}=\mathbf{h}+\mathbf{M}^{\mathrm{T}}\mathbf{z}$ for $\mathbf{h}\in\operatorname{\mathcal{G}}(\mathbf{M}^{\mathrm{T}})$ and $\mathbf{z}\in\mathbb{Z}^{d}$ the Fourier coefficients from the formula above can be collected with respect to congruence classes of the generating set $\operatorname{\mathcal{G}}(\mathbf{M}^{\mathrm{T}}$ and this yields

[TABLE]

for $\mathbf{y}\in\operatorname{\mathcal{P}}(\mathbf{M})$.

We finally take a closer look at the Fourier series

[TABLE]

and analyse its convergence. With $f\in A(\mathbb{T}^{d})\subset L^{2}(\mathbb{T}^{d})$ and because the Fourier coefficients $c_{\mathbf{h}}(\gamma)$ and $c_{\mathbf{h}}^{\mathbf{M}}(f)$ depend only on $\mathbf{h}\in\operatorname{\mathcal{G}}(\mathbf{M}^{\mathrm{T}})$ and not on $\mathbf{z}\in\mathbb{Z}^{d}$ the Fourier series

[TABLE]

converges and the result is at least in \operatorname{Sym}_{d}\bigl{(}L^{2}\bigl{(}\mathbb{T}^{d}\bigr{)}\bigr{)}. When we apply the differential operators $\overline{\operatorname{\operatorname{\nabla}_{\mathrm{Sym}}}_{\mathbf{k}}}^{\mathrm{T}}$ and \bigl{(}\overline{\operatorname{\operatorname{\nabla}_{\mathrm{Sym}}}_{\mathbf{k}}}^{\mathrm{T}}\mathcal{C}^{0}\operatorname{\operatorname{\nabla}_{\mathrm{Sym}}}_{\mathbf{k}}\bigr{)}^{-1} we differentiate the $L^{2}$ function once and integrate twice. The resulting function of interest is at least once weakly differentiable, i.e. in $H^{1}(\mathbb{T}^{d})^{d}$. Thus, another application of $\operatorname{\operatorname{\nabla}_{\mathrm{Sym}}}$ is admissible and the Fourier series (15) converges. Further (15) is a gradient field with mean zero, cf. (7), and the proof is concluded. ∎

A special choice for the space $V^{f}_{\mathbf{M}}$ comes from using the Dirichlet kernel $D_{\mathbf{M}}$, where $\Gamma^{\mathrm{p}}$ coincides with $\Gamma^{0}$, which also occurs in the derivation in [19].

Theorem 3.9**.**

For the Dirichlet kernel $D_{\mathbf{M}}$ the periodised Green operator $\Gamma^{\mathrm{p}}$ of (14) on $V^{D_{\mathbf{M}}}_{\mathbf{M}}$ coincides with the Green operator $\Gamma^{0}$.

Proof.

When we insert the formula for the Dirichlet kernel $D_{\mathbf{M}}$ into (14), the sum reduces to one single term which is exactly $\hat{\Gamma}^{0}_{\mathbf{k}}$ and thus the proof is completed. ∎

In contrast to the operator $\Gamma^{0}\mathcal{C}^{0}$, the periodised Green operator $\Gamma^{\mathrm{p}}\mathcal{C}^{0}$ corresponding to $V_{\mathbf{M}}^{f}$ is in general no longer a projection.

Theorem 3.10**.**

The periodised Green operator $\Gamma^{\mathrm{p}}\mathcal{C}^{0}$ corresponding to $V_{\mathbf{M}}^{f}$ is a projection operator if and only if $V_{\mathbf{M}}^{f}=V_{\mathbf{M}}^{D_{\mathbf{M}}}$, i.e. iff either $f$ or (one of) its orthonormalised translates is the Dirichlet kernel $D_{\mathbf{M}}$.

Proof.

Consider for a field \gamma\in\operatorname{Sym}_{d}\bigl{(}V_{\mathbf{M}}^{f}\bigr{)} the Fourier series

[TABLE]

for $\mathbf{y}\in\operatorname{\mathcal{P}}(\mathbf{M})$ and insert the definition of the periodised Green operator (14) to get

[TABLE]

From [27, Lemma 2 (iii)] we know that for $\mathbf{k}\in\mathbb{Z}^{d}$ we have that $\hat{\Gamma}^{0}_{\mathbf{k}}\mathcal{C}^{0}\hat{\Gamma}^{0}_{\mathbf{k}}\mathcal{C}^{0}=\hat{\Gamma}^{0}_{\mathbf{k}}\mathcal{C}^{0}$, i.e. that $\mathcal{C}^{0}\Gamma^{0}$ is a projection. This does not hold true for the mixed terms in (16), i.e. summands with $\mathbf{z}\neq\mathbf{z}^{\prime}$. These only vanish if $c_{\mathbf{k}}(f)=0$ for $\mathbf{k}\not\in\operatorname{\mathcal{G}}(\mathbf{M}^{\mathrm{T}})$, i.e. $f\in V_{\mathbf{M}}^{D_{\mathbf{M}}}$, cf. Theorem 3.9 ∎

Hence $\Gamma^{\mathrm{p}}$ is a projection if and only if $\Gamma^{\mathrm{p}}=\Gamma^{0}$ with respect to the corresponding generating set $\operatorname{\mathcal{G}}(\mathbf{M}^{\mathrm{T}})$. In addition the periodised Green operator $\Gamma^{\mathrm{p}}\mathcal{C}^{0}$ is bounded with the same bound as $\Gamma^{0}\mathcal{C}^{0}$.

Theorem 3.11**.**

Let the translates of $f$ be orthonormal and let \mathcal{C}^{0}\in\operatorname{SSym}_{d}\bigl{(}\mathbb{R}\bigr{)} be elliptic with constants $0<l_{\mathcal{C}^{0}}\leq u_{\mathcal{C}^{0}}<\infty$. Then the periodised Green operator $\Gamma^{\mathrm{p}}\mathcal{C}^{0}$ corresponding to $V^{f}_{\mathbf{M}}$, is bounded by

[TABLE]

for all \gamma\in\operatorname{Sym}_{d}\bigl{(}V^{f}_{\mathbf{M}}\bigr{)}.

Proof.

The Parseval equation together with the splitting $\mathbf{k}=\mathbf{h}+\mathbf{M}^{\mathrm{T}}\mathbf{z}$ with $\mathbf{h}\in\operatorname{\mathcal{G}}(\mathbf{M}^{\mathrm{T}})$ and $\mathbf{z}\in\mathbb{Z}^{d}$ yields

[TABLE]

The Cauchy-Schwarz theorem together with inserting the formula for the bracket sums (3) bounds this expression from above by

[TABLE]

A standard estimate for $\hat{\Gamma}^{0}_{\mathbf{k}}$ with $\mathbf{k}\in\mathbb{Z}^{d}$ is $\lVert\hat{\Gamma}^{0}_{\mathbf{k}}\rVert^{2}\leq\frac{1}{l_{\mathcal{C}^{0}}}$, see e.g. [27]. This leads to

[TABLE]

and Jensen’s inequality together with Lemma 2.3 c) results in

[TABLE]

With another application of the Parseval equation the desired estimate

[TABLE]

is obtained. ∎

The computation of the Green operator $\Gamma^{p}$ on $V^{f}_{\mathbf{M}}$ involves computing the value of the series

[TABLE]

for all $\mathbf{h}\in\operatorname{\mathcal{G}}(\mathbf{M}^{\mathrm{T}})$. This evaluation simplifies for functions $f$ having compact support in the frequency domain and where the series reduces to a sum over finitely many terms. In this case the operator can be evaluated exactly, i.e. without introducing any additional numerical error. This is the case for example for de la Vallée Poussin means, where each sum only consists of up to $4$ terms.

Functions that have compact support in space no longer allow for an exact evaluation of $\Gamma^{p}$. An example are Box splines in space domain which can be interpreted as finite elements integrated by only one quadrature point. In addition, the Box splines allow for finite elements which have different degrees of differentiability in directions other than the grid.

Brisard and Dormieux [8] derive a Green operator from an energy based formulation using a discretisation with element-wise constant finite elements. Their Green operator corresponds to using a Box spline of order zero in the approach here and is thus contained in the framework.

3.3 Discretisation of the Lippmann-Schwinger equation

With the definition of the periodised Green operator $\Gamma^{p}$ we can now proceed to derive a corresponding discretisation of the Lippmann-Schwinger equation. This derivation is split into two theorems. The discretised version of the space \mathcal{E}\bigl{(}\mathbb{T}^{d}\bigr{)} is given by

[TABLE]

This space allows to state the discretised version of the PDE correctly. For the following theorems we assume that \bigl{(}\mathcal{C}-\mathcal{C}^{0}\bigr{)}:(\varepsilon_{\mathbf{M}}+\varepsilon^{0})\in\operatorname{Sym}_{d}\bigl{(}A(\mathbb{T}^{d})\bigr{)} for \varepsilon_{\mathbf{M}}\in\mathcal{E}_{\mathbf{M}}^{f}\bigl{(}\mathbb{T}^{d}\bigr{)} so an interpolation on \operatorname{Sym}_{d}\bigl{(}V_{\mathbf{M}}^{f}\bigr{)} is possible with

[TABLE]

for all $\mathbf{x}\in\operatorname{\mathcal{P}}(\mathbf{M})$.

Additionally, a test functions $\gamma$ can be written as

[TABLE]

with $\mathbf{G}=(G_{\mathbf{y}})_{\mathbf{y}\in\operatorname{\mathcal{P}}(\mathbf{M})}\in\mathbb{C}^{m}$ and its discrete Fourier transform by $\mathbf{\hat{G}}=(\hat{G}_{\mathbf{h}})_{\mathbf{h}\in\operatorname{\mathcal{G}}(\mathbf{M}^{\mathrm{T}})}=\operatorname{\mathcal{F}}(\mathbf{M})\mathbf{G}$.

Let the strain \varepsilon_{\mathbf{M}}\in\mathcal{E}_{\mathbf{M}}^{f}\bigl{(}\mathbb{T}^{d}\bigr{)} be written in terms of translates of the fundamental interpolant as

[TABLE]

and again the discrete Fourier transform of the coefficient vector as $\mathbf{\hat{E}}=\operatorname{\mathcal{F}}(\mathbf{M})\mathbf{E}$.

Theorem 3.12**.**

Let the translates of $f$ be orthonormal, let \bigl{(}\mathcal{C}-\mathcal{C}^{0}\bigr{)}:(\varepsilon_{\mathbf{M}}+\varepsilon^{0})\in\operatorname{Sym}_{d}\bigl{(}A(\mathbb{T}^{d})\bigr{)}, let \gamma\in\operatorname{Sym}_{d}\bigl{(}V^{f}_{\mathbf{M}}\bigr{)}, and \varepsilon_{\mathbf{M}}\in\mathcal{E}_{\mathbf{M}}^{f}\bigl{(}\mathbb{T}^{d}\bigr{)}. Then it holds

[TABLE]

Proof.

Starting with the left-hand side of (18) applying the Parseval equation (2) to transform it to Fourier space yields the equal form

[TABLE]

where we make use of (9). Equation (6) together with the formula for the translate coefficients of the fundamental interpolant (13) and the splitting $\mathbf{k}=\mathbf{h}+\mathbf{M}^{\mathrm{T}}\mathbf{z}$ with $\mathbf{k},\mathbf{z}\in\mathbb{Z}^{d}$ and $\mathbf{h}\in\operatorname{\mathcal{G}}(\mathbf{M}^{\mathrm{T}})$ results in the expressions

[TABLE]

Inserting these equations into the expression above yields

[TABLE]

Collecting the terms depending on $\mathbf{z}$ and employing the bracket sums (3) to simplify the expression, one obtains

[TABLE]

Since by assumption translates of $f$ are orthonormal they fulfil \bigl{[}\bigl{\{}\lvert c_{\mathbf{k}}(f)\rvert^{2}\bigr{\}}\bigr{]}_{\mathbf{h}}^{\mathbf{M}}=\frac{1}{m} by Lemma 2.3 c) and we get the desired result. ∎

The following theorem states the result of a Galerkin projection of (10) onto the space of translates.

Theorem 3.13**.**

Let the translates of $f$ be orthonormal, let \mathcal{C}\in\operatorname{SSym}_{d}\bigl{(}A\bigl{(}\mathbb{T}^{d}\bigr{)}\bigr{)}, and let \varepsilon_{\mathbf{M}}\in\mathcal{E}_{\mathbf{M}}^{f}\bigl{(}\mathbb{T}^{d}\bigr{)}. Then $\varepsilon_{\mathbf{M}}$ fulfils the weak form

[TABLE]

for all \gamma\in\operatorname{Sym}_{d}\bigl{(}V_{\mathbf{M}}^{f}\bigr{)} if and only if

[TABLE]

for all $\mathbf{x}\in\operatorname{\mathcal{P}}(\mathbf{M})$, where $\Gamma^{\mathrm{p}}_{\mathbf{y}}$ is defined in Definition 3.6.

Proof.

With \varepsilon_{\mathbf{M}}\in\mathcal{E}_{\mathbf{M}}^{f}\bigl{(}\mathbb{T}^{d}\bigr{)}\subset\operatorname{Sym}_{d}\bigl{(}V_{\mathbf{M}}^{f}\bigr{)}\subset\operatorname{Sym}_{d}\bigl{(}A(\mathbb{T}^{d})\bigr{)} and \mathcal{C}\in\operatorname{SSym}_{d}\bigl{(}A\bigl{(}\mathbb{T}^{d}\bigr{)}\bigr{)} it follows that \bigl{(}\mathcal{C}-\mathcal{C}^{0}\bigr{)}:(\varepsilon_{\mathbf{M}}+\varepsilon^{0})\in\operatorname{Sym}_{d}\bigl{(}A(\mathbb{T}^{d})\bigr{)} and hence the assumptions of Theorem 3.12 are fulfilled. Therefore (19) is equivalent to

[TABLE]

with the notation from (17). A necessary and sufficient condition for (21) to hold true is that it is fulfilled for all

[TABLE]

for all $\mathbf{h}\in\operatorname{\mathcal{G}}(\mathbf{M}^{\mathrm{T}})$ and $\mathbf{y}\in\operatorname{\mathcal{P}}(\mathbf{M})$ and $p,q\in\{1,\dots,d\}$. The vector $\alpha_{p}\in\mathbb{R}^{d}$ denotes the $p$-th unit vector and $\beta_{pq}\coloneqq 1-\frac{1}{2}\delta_{pq}$ normalizes the resulting matrix. This parametrization is the trigonometric basis of \operatorname{Sym}_{d}\bigl{(}V_{\mathbf{M}}^{f}\bigr{)} on the pattern $\operatorname{\mathcal{P}}(\mathbf{M})$.

Hence an equivalent condition stems from looking at (21) component-wise, i.e.

[TABLE]

bearing in mind the necessary complex conjugate, for all $\mathbf{y}\in\operatorname{\mathcal{P}}(\mathbf{M})$. This, however, is an inverse discrete Fourier transform on the pattern $\operatorname{\mathcal{P}}(\mathbf{M})$ and together with Definition 3.6 and setting

[TABLE]

yields

[TABLE]

The coefficients \bigl{(}\tilde{B}_{\mathbf{y}}\bigr{)}_{\mathbf{y}\in\operatorname{\mathcal{P}}(\mathbf{M})}\coloneqq\frac{1}{\sqrt{m}}\overline{\operatorname{\mathcal{F}}}(\mathbf{M})^{\mathrm{T}}\bigl{(}\hat{\tilde{B}}_{\mathbf{h}}\bigr{)}_{\mathbf{h}\in\operatorname{\mathcal{G}}(\mathbf{M}^{\mathrm{T}})} can now be interpreted as coefficients of translates of the fundamental interpolant, i.e. it holds

[TABLE]

for all $\mathbf{x}\in\mathbb{T}^{d}$. By Definition 3.6 the operator $\Gamma^{\mathrm{p}}$ acts as a Fourier multiplier with Fourier coefficients (14). This transforms the above equation to

[TABLE]

The coefficients $B_{\mathbf{y}}$ were chosen such that they coincide with the function values of \bigl{(}\mathcal{C}(\mathbf{y})-\mathcal{C}^{0}\bigr{)}:\bigl{(}\varepsilon_{\mathbf{M}}(\mathbf{y})+\varepsilon^{0}\bigr{)} at points $\mathbf{y}\in\operatorname{\mathcal{P}}(\mathbf{M})$. Likewise, $\varepsilon_{\mathbf{M}}$ coincides in the points $\mathbf{y}\in\operatorname{\mathcal{P}}(\mathbf{M})$ with the coefficients $E_{\mathbf{y}}$. Inserting these relations one obtains

[TABLE]

which yields the desired result. ∎

When discretising the PDE (10) in a similar way, one arrives at the following discretised form.

Theorem 3.14**.**

Let the translates of $f$ be orthonormal and let \mathcal{C}\in\operatorname{SSym}_{d}\bigl{(}A\bigl{(}\mathbb{T}^{d}\bigr{)}\bigr{)} and let \varepsilon_{\mathbf{M}}\in\mathcal{E}_{\mathbf{M}}^{f}\bigl{(}\mathbb{T}^{d}\bigr{)}. Then $\varepsilon_{\mathbf{M}}$ fulfils the weak form (10)

[TABLE]

for all \gamma\in\operatorname{Sym}_{d}\bigl{(}V_{\mathbf{M}}^{f}\bigr{)} if and only if

[TABLE]

for all $\mathbf{x}\in\operatorname{\mathcal{P}}(\mathbf{M})$.

Proof.

The proof follows the same steps as for the Theorems 3.12 and 3.13 and we omit it therefore. ∎

In Remark 3.5 we already mentioned that in the continuous case the Lippmann-Schwinger equation (11) and the variational equation (10) coincide. This is, as shown in [27], also the case when using trigonometric collocation for the discretisation. With the equations (20) and (22) using spaces of translates this is in general no longer the case. When looking at the identity (12) one can see this rather quickly.

Remark 3.15**.**

For \varepsilon_{\mathbf{M}}\in\mathcal{E}_{\mathbf{M}}^{f}\bigl{(}\mathbb{T}^{d}\bigr{)} it holds true that

[TABLE]

almost everywhere if and only if $V_{\mathbf{M}}^{f}=V_{\mathbf{M}}^{D_{\mathbf{M}}}$, i.e. if $f$ can be written as a sum of translates of the Dirichlet kernel $D_{\mathbf{M}}$.

Proof.

The proof in [27, Proposition 3] uses that $\Gamma^{0}\mathcal{C}^{0}$ projects a constant function $g$ onto the function that is [math] almost everywhere, i.e. it is only characterized by the Fourier coefficient $c_{\mathbf{0}}(g)$. When interpolating in $V^{f}_{\mathbf{M}}$ this is in general no longer the case and an application of $\Gamma^{\mathrm{p}}\mathcal{C}^{0}$ does not result in the zero function. ∎

The connections between the variational formulation (10) and the Lippmann-Schwinger equation (11) and their discretisations on spaces of translates are summarized in Figure 2. The continuous equations are shown to be equivalent in [27, Proposition 3]. The same holds also true for the equations discretised on $V_{\mathbf{M}}^{f}$ if and only if $V_{\mathbf{M}}^{f}=V_{\mathbf{M}}^{D_{\mathbf{M}}}$ For the special case $f=D_{\mathbf{M}}$ and $\mathbf{M}$ a diagonal matrix, i.e. for a tensor product grid, this equivalence is already proven in [27, Proposition 12]. Box splines allow to solve the equations in terms of (simplified) finite elements, which generalizes the constant finite element approach of [8]. This emerges when discretising the Lippmann-Schwinger equation with $f=B_{\operatorname{Id}}$.

4 Numerics

In this section we study the effect of choosing different functions for the translation invariant spaces to illustrate the capabilities of the generalization presented in this paper. In the publication [3] the authors study the influence of different patterns on the solution quality and their numerical effects.

A prototypical structure that is introduced in [3] is the generalized Hashin structure, a geometry that is based on publications of Milton, see [16]. It consists of two confocal ellipses embedded in a surrounding material, see Figure 3. The centre ($\Omega_{c}$) and coating ($\Omega_{e}$) ellipses have isotropic behaviour and the matrix material ($\Omega_{m}$) is built in such a way that it is unaffected by the inclusion for a chosen macroscopic strain $\varepsilon^{0}$. For this special kind of structure analytic expressions for the strain field $\varepsilon$ and the action of the effective matrix $\mathcal{C}^{\mathrm{eff}}:\varepsilon^{0}$ are known.

In the following we take exactly the same parameters as in [3], i.e. for the ellipses with parametrization

[TABLE]

and $(x_{1},x_{2})^{\mathrm{T}}\in\mathbb{T}^{2}$ we choose $c_{1}=0.05$, $c_{2}=0.35$ and for the inner and outer ellipsis $\rho=0$ and $\rho=0.09$, respectively. The structure is then rotated by $60^{\circ}$. For the inner and outer ellipsis we choose isotropic material laws with Poisson’s ratio $\nu=0.3$ in both ellipses and the matrix material and Young’s moduli $E=1$ and $E=10$ for the inner and outer ellipsis, respectively. The material law for the surrounding matrix material can then be determined by the formulae in [3, Section 4.2].

A solution for the first component of the strain field $\varepsilon_{11}$ is depicted in Figure 3 and taken for instructional purposes directly from [3].

All numerical results in this section were obtained by solving (20) with a fixed-point iteration on $E_{\mathbf{y}}$ as described in [19] up to an relative error of $10^{-10}$ using a Cauchy criterion.

4.1 De la Vallée Poussin means

In [3] the authors study the influence of the pattern matrix $\mathbf{M}$ on the solution field and the quality of the effective matrix $\mathcal{C}^{\mathrm{eff}}$. We take the following pattern matrices:

[TABLE]

They correspond to a tensor product grid ($\mathbf{M}_{1}$), a tensor product grid rotated by $45^{\circ}$ ($\mathbf{M}_{2}$) and the minimal $\ell^{2}$-error ($\mathbf{M}_{3}$) achieved. The matrices $\mathbf{M}_{1}$ and $\mathbf{M}_{3}$ have a determinant of $2^{14}$ whereas the so-called quincux pattern $\operatorname{\mathcal{P}}(\mathbf{M}_{2})$ has $2^{13}$ sampling points.

First we study how using de la Vallée Poussin means changes the quality of the effective stiffness matrix and the strain field. The Box splines used in frequency domain to define the means allow to use different slopes in each direction of the de la Vallée Poussin mean and thus one can expect to reduce the Gibbs phenomenon in different directions. For the following study the functions $f_{\mathbf{M}_{i},\alpha}$ are parametrized with $\alpha=(\alpha_{1},\alpha_{2})^{\mathrm{T}}$ and $\alpha_{1},\alpha_{2}\in[0,0.5]$ for $i=1,\dots,3$.

The parameters $\alpha_{1}$ and $\alpha_{2}$ correspond to damping the Fourier coefficients of the de la Vallée Poussin mean along the directions $\mathbf{M}^{\mathrm{T}}(1,0)^{\mathrm{T}}$ and $\mathbf{M}^{\mathrm{T}}(0,1)^{\mathrm{T}}$, respectively. In the space domain they introduce a better localization [11] along $\mathbf{M}^{-1}(0,1)^{\mathrm{T}}$ and $\mathbf{M}^{-1}(1,0)^{\mathrm{T}}$, respectively.

The result of these experiments regarding the effective stiffness matrix is shown in Figure 4. The relative effective error

[TABLE]

is depicted with $E_{\mathbf{y}}$ from (20). Parameters $\alpha_{1}=\alpha_{2}=0$ correspond to the modified Dirichlet kernel $f_{\mathbf{M},\mathbf{0}}$ and $\alpha_{1}=\alpha_{2}=0.5$ correspond to the Fejér kernel $f_{\mathbf{M},\frac{1}{2}\mathbf{1}}$.

For $\mathbf{M}_{1}$, the error in the effective stiffness matrix changes from $0.0037$ for the modified Dirichlet kernel to the optimum at $\alpha_{1}=0.45$ and $\alpha_{2}=0$ with $e_{\mathrm{eff}}=0.0017$. For $\mathbf{M}_{2}$ it changes from $0.0034$ to $e_{\mathrm{eff}}=0.0013$ for $\alpha_{1}=0.25$ and $\alpha_{2}=0$. The solution for $\mathbf{M}_{3}$ can be improved upon by de la Vallée Poussin means with $\alpha_{1}=0.4$ and $\alpha_{2}=0$ which results in $e_{\mathrm{eff}}=0.0036$ of the Fourier approach on patterns being reduced to $e_{\mathrm{eff}}=0.0024$. When approaching either $\alpha_{1}=0.5$ or $\alpha_{2}=0.5$ the value of $e_{\mathrm{eff}}$ drastically increases. The errors for the solution using the Dirichlet kernel and the modified Dirichlet kernel, i.e. for $\alpha_{1}=\alpha_{2}=0$, are almost the same and thus not marked in the plots above.

The study in [3] suggests that by changing the pattern one can improve the quality of the solution and the effective stiffness matrix. The experiments from above show that one can improve these results even further by extending the theory to translation invariant spaces. Especially for the tensor product grid, i.e. for data given on a regular voxel grid like for example from a computer tomography image, a suitable choice of the space of translates could here reduce the error in the effective stiffness by $55\%$.

Figure 5 shows the logarithmic error of the analytic solution $\tilde{\varepsilon}$ to solution $\varepsilon$, i.e. $e_{\text{log}}\coloneqq\log(1+\lvert\tilde{\varepsilon}_{11}-\varepsilon_{11}\rvert)$. The top row shows the error corresponding to the de la Vallée Poussin means with parameters $\alpha_{1}$ and $\alpha_{2}$ which give the smallest $e_{\mathrm{eff}}$. The middle row shows the error using the Dirichlet kernel. For illustration purposes each pixel has the form of the corresponding unit cell $\mathbf{M}^{-1}[-\frac{1}{2},\frac{1}{2}]^{2}$ centred at each pattern point $\mathbf{y}\in\operatorname{\mathcal{P}}(\mathbf{M})$.

The value of the relative $\ell_{2}$-error is given by the formula

[TABLE]

While the error in the effective stiffness matrix can be drastically improved using de la Vallée Poussin means, the $\ell^{2}$-error gets larger, however, only be a few percent. The strain field stemming from the numerical computation using de la Vallée Poussin means shows less Gibbs phenomena as can be seen most prominently in the solution for pattern matrix $\mathbf{M}_{2}$.

The decrease of $e_{\mathrm{eff}}$ is caused by the smoothing of the Gibbs phenomena around discontinuities of the solution with higher $\alpha_{1}$. With smaller de la Vallée Poussin parameters $\alpha_{1}$ and $\alpha_{2}$ the polynomial reproduction is better and thus interfaces (edges) are sharper. In total they introduce a trade off between damping Gibbs phenomena and sharpness of the interfaces.

4.2 Periodised Box Splines

Finite elements with one quadrature point are directly included in this framework. They are obtained by choosing a suitable Box spline $f=B_{\mathbf{M},\Xi}$ as ansatz function for the space of translates $V_{\mathbf{M}}^{f}$. In the bottom row of Figure 5 the logarithmic error between the numerical solution and the analytical solution is depicted, using the matrices $\mathbf{M}_{1}$, $\mathbf{M}_{2}$, $\mathbf{M}_{3}$ from above. For all computations the Box spline $B_{\mathbf{M},\Xi}$ with

[TABLE]

is used. This ansatz function corresponds to a finite element of third order with reduced integration. The bracket sums in (14) are precomputed with $33^{2}$ terms each.

A comparison of the relative error in the effective matrix $e_{\mathrm{eff}}$ between the Box splines (bottom row) and the Dirichlet kernel (middle row) shows that for matrices $\mathbf{M}_{1}$ and $\mathbf{M}_{2}$ the error in the effective matrix can be reduced by about $25\%$, whereas the error is approximately tripled for $\mathbf{M}_{3}$. The $\ell^{2}$-error $e_{\ell^{2}}$ is increased in all three cases. The pattern matrix $\mathbf{M}_{3}$ was optimized for the Dirichlet kernel with respect to the $\ell^{2}$-norm. For Box splines this pattern leads to an aliasing effect reducing the quality of the solution.

The choice of the Box spline was not optimized and the quality of the effective matrix might be further increased by a different choice of $\Xi$. This corresponds to a fine tuning with respect to dominant directions in the pattern unit cell \mathbf{M}^{-1}\bigl{[}-\frac{1}{2},\frac{1}{2}\bigr{)}^{d}.

5 Conclusion

The introduced framework unifies and analyses for the first time the truncated Fourier series approach and finite element ansatz functions. In this framework the periodised Green operator possesses the same properties as the Green operator of the Galerkin method of Vondřejc et.al [27]. The projection operator $\Gamma^{0}\mathcal{C}^{0}$ emerges as a special case of Dirichlet kernel translates. The periodised Green operator $\Gamma^{\mathrm{p}}\mathcal{C}^{0}$ can further be characterised to be a projection if and only if the space $V_{\mathbf{M}}^{f}$ is the one derived from a Dirichlet kernel. The finite element method emerges for certain Box splines and thus the constant finite elements of Brisard and Dormieux [8] are included. For these Box splines the infinite Bracket sums have to be precomputed up to a given precision, but the algorithm has the same complexity as the Fourier framework.

Finite elements with more sophisticated quadrature rules can also be viewed in terms of this framework. However, their performance is an open question for future work. How to choose a certain Box spline especially with respect to anisotropies present in the data is a point for further analysis. Convergence of the numerical discretisation towards the continuous case and the convergence of the algorithm will be dealt with in a paper on our road map. Both seem to be strongly suggested by the numerical results. Finally the periodic multiresolution analysis can be used to extend this framework in order to exploit sparsity properties of given data. This is also a point for future work.

Acknowledgement

The authors would like to thank Bernd Simeon and Gabriele Steidl for their valuable comments on preliminary versions of this manuscript and fruitful discussions.