###### Abstract

The computational homogenization of hyperelastic solids in the geometrically
nonlinear context has yet to be treated with sufficient efficiency in order to
allow for real-world applications in true multiscale settings. This problem is
addressed by a problem specific surrogate model founded on a reduced basis
approximation of the deformation gradient on the microscale. The setup phase
is based upon a snapshot POD on deformation gradient fluctuations, in contrast
to the widespread displacement-based approach. To this end, the space of
relevant macroscopic stretch tensors is sampled in an efficient manner in
order to lower the computational offline costs. Numerical results show speed-
up factors in the order of 5–100 and significantly improved robustness while
retaining good accuracy. An open-sourced demonstrator tool with 50 lines of
code emphasizes the simplicity and efficiency of the method.

###### keywords:

computational homogenization; large strain; finite deformation; geometric
nonlinearity; reduced basis; reduced order model; sampling

24 1 39 Received: 31 March 2019 Finite strain homogenization using a reduced
basis and efficient sampling Oliver Kunc 1,† and Felix Fritzen 2,†,‡ Oliver
Kunc and Felix Fritzen These authors contributed equally to this work.
Corresponding author. 74Q05, 74B20, 74S30

##  1 Introduction

###  1.1 Purpose

The description of solid mechanics under finite strains is of particular
interest in both academia and industry. It allows for accurate descriptions of
rotations and stretches under mild assumptions. Thus, many geometric effects
can be captured. For instance, alignments and rearrangements of the respective
structures may trigger pronounced stiffening or softening effects.

In such cases, in which rotations and deformations are not suitable for
linearization, dissipative effects also play a notable role for many
materials. Regardless of the kind of dissipation involved in a certain
process, hyperelasticity usually persists to a certain extend. Therefore, it
is worthwile investigating this comparatively simple case at first, before
introducing history dependence into the description. Prominent examples of
materials that require a hyperelastic description at finite strains include
carbon black-filled rubber (Rendek and Lion, 2010) and amorphous glassy
polymers (Nguyen et al., 2016), to name just two.

The main purpose of this work is the _computationally efficient quasi-static
homogenization of hyperelastic solids with full account for geometric
nonlinearities_. The employed methodology is twofold: first, a Reduced Basis
(RB) model for the microscopic problem is established. Once set up, it enables
more efficient evaluations of the homogenized material response as compared to
the Finite Element Method (FEM). Second, an efficient strategy for sampling of
the space of macroscopic kinematic states is proposed. This renders the setup
phase of the RB model more rational.

###  1.2 State of the art

Efficiently determining the overall solid-mechanical properties of
microstructures has been investigated for decades, and a large body of
literature is available. Comprehensive review articles, such as Geers and
Yvonnet (2016) and Saeb et al. (2016), summarize the progress. Here, attention
is restrained to few methods most similar or relevant to the present work.

The FE2 method Feyel (1999) is theoretically capable of performing realistic
two-scale simulations with arbitrary accuracy. Therefore, it serves as a
reference method in the context of first-order homogenization based on the
assumption of separated length scales. In The FE2, the evaluations of the
unknown macroscopic constitutive law are approximated by microscopic FE
simulations. However, this comes along with computational costs that quickly
exceed the capabilities of common workstations, both at present and in the
foreseeable future. Roughly speaking, the computational effort required on the
microscale multiplies with that of the macroscale, hence the method’s name. It
is thus worthwhile to develop order reduction methods for the microscopic
problem.

A common approach within the field of computational homogenization (and well
beyond) is to extract essential information from provided in silico data. To
this end, schemes based on the Proper Orthogonal Decomposition (POD) compute
correlations within snapshot data. Such methods include the R3M (Yvonnet and
He, 2007) and can be further enhanced by the use of, e.g., the EIM, as in
Radermacher and Reese (2016). Numerical comparisons of various schemes were
conducted in Radermacher and Reese (2013) and Soldner et al. (2017). To the
best of the authors’ knowledge, all published POD-based methods addressing the
finite strain hyperelastic problem choose to reduce the number of degrees of
freedom (DOF) of the displacement field. This results in sometimes significant
speed-ups. Another important feature is that they allow for reconstruction of
the microscopic displacement fields.

Still, the solution of the reduced equations remains a complex task. It
requires evaluations of material laws and numerical integration over the
microstructure. Promising progress is made in the field of efficient
integration schemes, see for instance An et al. (2009) and Hernández et al.
(2017). A main reason for the speed-up of these methods is the reduced number
of function evaluations.

The highest speed-ups are achievable if the computational effort of the
determination of effective microstructural responses can be fully decoupled
from underlying microstructural discretizations. Such homogenization methods
directly approximate the effective material law by means of a dedicated
numerical scheme. Technically, this can be seen as the direct surrogation of
unknown functions, e.g. of the effective free energy or stress. For instance,
the Material Map Temizer and Zohdi (2007) interpolates the coefficients of an
assumed macroscopic material model. Another example is the NEXP method Yvonnet
et al. (2013), the effective stored energy density is approximated using a
tensor product of one-dimensional splines. The authors treated the case of
small strains by introducing the RNEXP method Fritzen and Kunc (2018), where
the effective stored energy is interpolated by a dedicated kernel scheme.

However, interpolatory and regressional methods suffer the inherent drawback
of not providing any explicit information on the microscale. For instance,
microscopic displacement or stress fields cannot be reconstructed from the
solutions of macroscopic interpolation. Another important open question is how
to provide the supporting data points for the interpolation in an efficient
manner. The data at these points is usually provided by the solution of a
full-order model (FOM) and come along with the corresponding numerical costs.
Hence the positions of the data points in the parameter space should be chosen
carefully, as unnecessary or redundant solutions of the FOM should be avoided.
On the other hand, too sparsely seeded points might not capture the
homogenized properties of the microstructure appropriately.

###  1.3 Main contributions and outline

The present work generalizes parts of the previous work Fritzen and Kunc
(2018) to the finite strain regime. It aims at reducing the computational
complexity of the problem associated with the microstructure by means of a
_Reduced Basis approximation of the microscopic deformation gradient_. The
basis is obtained with the aid of a _POD of snapshots of fluctuation fields of
the deformation gradient_. Thus, the application of the RB model _does not
necessitate the computation of gradients of displacement fields_ , and even
does not require the displacements to be available at all. In other words,
microscopic displacement fields are completely avoided. However, they _can be
reconstructed_ from the RB approximation of the deformation gradient, uniquely
up to rigid body motion.

Another key advantage is the _sleek implementation_ of the method. A
_demonstrator_ containing a minimum working example of the RB model with 50
lines of MATLAB/Octave code, (Kunc, 2019).

As for the setup phase, the snapshot data is created by means of an _efficient
sampling procedure_. To this end, a _linear parameter space_ is identified,
within which the sampling sites are placed based upon _physical
interpretation_. This allows for _controllability_ of the resolution of
certain key characteristics of the effective material response while keeping
the total number of samples within bounds.

The Reduced Basis method is presented in Section 2. The basis identification
is based on the sampling strategy developed in Section 3. Numerical examples
are presented in Section 4. Both the numerical and the theoretical findings
are summarized and discussed in Section 5.

###  1.4 Notation

The set of real numbers and the subset of positive numbers greater than zero
are denoted by ${\mathbb{R}}$ and ${\mathbb{R}}_{+}$, respectively. Matrices
are marked by two underlines and vectors by one underline, e.g.
$\underline{\underline{A}}$, $\underline{a}$. Vectors are assumed to be
columns and the dot product of two vectors of the same size is understood as
the Euclidean scalar product,
$\underline{x}\cdot\underline{y}=\underline{x}^{\sf T}\underline{y}$. First
order and second order tensors in coordinate free description are denoted by
bold letters, e.g. $A$, $a$. No conclusion can be drawn on the order of a
tensor based on its capitalization. Here, the underlying space is always the
Euclidean space ${\mathbb{R}}^{3}$ with its standard basis. First order and
second order tensors can be represented as vectors and matrices, e.g.
$\mbox{\boldmath$A$}\widehat{=}\underline{A}\in{\mathbb{R}}^{3}$ and
$\mbox{\boldmath$B$}\widehat{=}\underline{\underline{B}}\in{\mathbb{R}}^{3\times
3}$ respectively. Norms of vectors and matrices respectively denote the
Euclidean and the Frobenius norm. The norm of a tensor of second order equals
the norm of its matrix representation for the chosen basis. Fourth order
tensors are denoted by blackboard bold symbols other than ${\mathbb{R}}$, e.g.
${\mathbb{C}}$ and ${\mathbb{I}}$. Components of tensors of order $M$ are with
respect to the Euclidean tensorial basis
$\mbox{\boldmath$e$}^{(1)}\otimes\mbox{\boldmath$e$}^{(2)}\otimes\dots\otimes\mbox{\boldmath$e$}^{(M)}$,
e.g. $A_{ij}$, $B_{ij}$ for second order tensors $A$, $B$ and
$C_{ijkl}$,$C^{\prime}_{ijkl}$ for ${\mathbb{C}}$, ${\mathbb{C}}^{\prime}$.
The following contractions are defined:

$\displaystyle\mbox{\boldmath$A$}\cdot\mbox{\boldmath$B$}$
$\displaystyle=\sum_{i,j=1}^{3}A_{ij}B_{ij}\,,$
$\displaystyle\mbox{\boldmath$A$}\cdot{\mathbb{C}}$
$\displaystyle=\sum_{i,j,k,l=1}^{3}A_{ij}C_{ijkl}\,\mbox{\boldmath$e$}^{(k)}\otimes\mbox{\boldmath$e$}^{(l)}\,,$
$\displaystyle{\mathbb{C}}\cdot\mbox{\boldmath$B$}$
$\displaystyle=\sum_{i,j,k,l=1}^{3}C_{ijkl}B_{kl}\,\mbox{\boldmath$e$}^{(i)}\otimes\mbox{\boldmath$e$}^{(j)}\,,$
$\displaystyle{\mathbb{C}}\cdot{\mathbb{C}}^{\prime}$
$\displaystyle=\sum_{i,j,k,l,m,n=1}^{3}C_{ijmn}C_{mnkl}^{\prime}\,\mbox{\boldmath$e$}^{(i)}\otimes\mbox{\boldmath$e$}^{(j)}\otimes\mbox{\boldmath$e$}^{(k)}\otimes\mbox{\boldmath$e$}^{(l)}\,.$

Let $\varOmega\subset{\mathbb{R}}^{3}$ be the domain occupied by a physical
body undergoing elastic deformations, and let $\varOmega_{0}$ be its initial
configuration. Then $x$ and $X$ describe the coordinates of material points
within the current configuration $\varOmega$ and within the reference state
$\varOmega_{0}$ respectively. Their difference is the displacement
$\mbox{\boldmath$u$}=\mbox{\boldmath$x$}-\mbox{\boldmath$X$}$, see Figure 1.
The gradient of a vector field
$\mbox{\boldmath$v$}=\mbox{\boldmath$v$}(\mbox{\boldmath$X$})$ is defined as a
_right gradient_ and denoted by
$\displaystyle\frac{\displaystyle\partial\mbox{\boldmath$v$}}{\displaystyle\partial\mbox{\boldmath$X$}}=\mbox{\boldmath$v$}\otimes\nabla_{\rm
X}$. The divergence of a second order tensor field is the vector field
resulting from _row-wise_ divergence. The boundaries of the respective
configurations are denoted by $\partial\varOmega$ and $\partial\varOmega_{0}$.
The set of square-integrable Lebesgue functions on the reference domain is
tagged $L^{2}(\varOmega_{0})$.

Figure 1: Initial ($\varOmega_{0}$) and current ($\varOmega$) configuration,
with elementary kinematic quantities.

The displacement gradient
$\mbox{\boldmath$H$}=\mbox{\boldmath$u$}\otimes\nabla_{\rm X}$ and the
deformation gradient
$\mbox{\boldmath$F$}=\mbox{\boldmath$x$}\otimes\nabla_{\rm X}$ are related
through $\mbox{\boldmath$F$}=\mbox{\boldmath$H$}+\mbox{\boldmath$I$}$, where
$I$ is the second order identity tensor in three dimensions. The determinant
$J=\det{\mbox{\boldmath$F$}}$ measures the relative volumetric change due to
the present deformation.

Unimodular quantities, i.e. second order tensors with determinant one, may be
emphasized by a hat, e.g.
$\widehat{\mbox{\boldmath$F$}}=J^{-1/3}\mbox{\boldmath$F$}$. This
multiplicative decomposition is sometimes attributed to Flory Flory (1961) and
also goes by the name Dilatational-Deviatoric Multiplicative Split (DDMS).

In the two-scale context, overlined symbols represent quantities on the
macroscopic scale, e.g. $\overline{\mbox{\boldmath$A$}}$,
$\overline{\mbox{\boldmath$a$}}$, while symbols without overline correspond to
their microscopic counterpart, e.g. $A$, $a$. Equivalently, macroscopic
quantities are called _global_ and microscopic ones are called _local_. The
volume average of a general local field $\varphi$

$\displaystyle\left<\varphi\right>$
$\displaystyle=\left<\varphi(\bullet)\right>=\frac{1}{|\varOmega_{0}|}\int_{\varOmega_{0}}\varphi(\bullet)\,{\rm
d}V$ (1)

is essential to the theory. The dependence of a microscopic quantity $A$ on
both the microscopic coordinates $X$ and a macroscopic quantity
$\overline{\mbox{\boldmath$B$}}$ is denoted by
$\mbox{\boldmath$A$}=\mbox{\boldmath$A$}(\mbox{\boldmath$X$};\overline{\mbox{\boldmath$B$}})$.
In such a case, the components of the macroscopic quantity
$\overline{\mbox{\boldmath$B$}}$ are called _parameters_ of the microscopic
function $\mbox{\boldmath$A$}(\bullet;\overline{\mbox{\boldmath$B$}})$. The
application of the volume averaging operator is abbreviated by
$\left<\mbox{\boldmath$A$}\right>=\left<\mbox{\boldmath$A$}(\bullet;\overline{\mbox{\boldmath$B$}})\right>$.
The case of a concatenated function
$f(\mbox{\boldmath$A$})=f(\mbox{\boldmath$A$}(\mbox{\boldmath$X$};\overline{\mbox{\boldmath$B$}}))$
is analogous, i.e.
$\left<f\right>=\left<f(\mbox{\boldmath$A$})\right>=\left<f(\mbox{\boldmath$A$}(\bullet;\overline{\mbox{\boldmath$B$}}))\right>$,
regardless of the tensorial order of the image of the function $f$.

###  1.5 Material models

In this work hyperelastic materials are investigated. They are characterized
by _stored energy density functions_ $W=W(\mbox{\boldmath$F$})$. The first
Piola-Kirchhoff stress

$\displaystyle\mbox{\boldmath$P$}(\mbox{\boldmath$F$})$
$\displaystyle=\displaystyle\frac{\displaystyle\partial
W}{\displaystyle\partial\mbox{\boldmath$F$}}(\mbox{\boldmath$F$})$ (2) and the
corresponding fourth-order stiffness tensor
$\displaystyle{\mathbb{C}}(\mbox{\boldmath$F$})$
$\displaystyle=\displaystyle\frac{\displaystyle\partial^{2}W}{\displaystyle\partial{\mbox{\boldmath$F$}}^{2}}(\mbox{\boldmath$F$})$
(3)

characterize the material response.

Henceforth, for reasons of readability, the stored energy density function $W$
will be spoken of as an energy and the terms _stored_ and _density_ will not
always be mentioned explicitly. In the infinitesimal strain framework,
hyperelastic energies have been formulated that model deformation plasticity
(e.g., Fritzen and Kunc, 2018; Yvonnet et al., 2013; Bilger et al., 2005).
Although these models are only valid for purely proportional loading
conditions, they provide means to simulate highly nonlinear material behavior
in certain scenarios comparably easy within the context of hyperelasticity.
Note that genuine dissipative processes require additional state describing
variables with corresponding evolution laws.

The proposed method is _suitable for any type of hyperelastic constitutive
law_. As the modeling of complex material behavior is not the main focus of
this study, the Neo-Hookean law

$\displaystyle W(\mbox{\boldmath$F$})=W_{\rm
DDMS}(J,\widehat{\mbox{\boldmath$F$}})$
$\displaystyle=\frac{K}{4}\left[(J-1)^{2}+({\rm
ln}\,J)^{2}\right]+\frac{G}{2}\left({\rm
tr}(\,\widehat{\mbox{\boldmath$F$}}^{\sf
T}\widehat{\mbox{\boldmath$F$}})-3\right)$ (4)

is used, with $K$ the bulk modulus and $G$ the shear modulus. The volumetric
part of the energy is taken from Doll and Schweizerhof (1999). Using the DDMS,
a decoupled dependence on the volumetric and isochoric part of the deformation
is assumed, which is a common way to model the distinct material behavior with
respect to these two contributions, see e.g. (Simo, 1988).

###  1.6 Problem setting of first order homogenization

Assuming stationarity and separability of scales, the following coupled,
deformation-driven problems can be derived by means of asymptotic expansion of
the displacement $u$ and subsequent first order approximation. This procedure
is carried out in Pruchnicki (1998) with much detail. Here, the technical
process is omitted and only the resulting equations are stated.

####  1.6.1 Macroscopic problem

Balance of linear momentum

$\displaystyle{\rm
Div}_{\rm\overline{X}}(\overline{\mbox{\boldmath$P$}})+\overline{\mbox{\boldmath$b$}}$
$\displaystyle=\mbox{\boldmath$0$},$ (5)

where $\overline{\mbox{\boldmath$b$}}$ denote bulk forces, and balance of
angular momentum

$\displaystyle\overline{\mbox{\boldmath$F$}}^{-1}\overline{\mbox{\boldmath$P$}}$
$\displaystyle=\overline{\mbox{\boldmath$P$}}^{\sf
T}\overline{\mbox{\boldmath$F$}}^{\sf-T}\text{,}$ (6)

along with well-posed boundary conditions constitute the _macroscopic boundary
value problem_. This system of equations is closed by means of the
hyperelasticity law, cf. (2),

$\displaystyle\overline{\mbox{\boldmath$P$}}(\overline{\mbox{\boldmath$F$}})$
$\displaystyle=\displaystyle\frac{\displaystyle\partial\overline{W}}{\displaystyle\partial\overline{\mbox{\boldmath$F$}}}.$
(7)

Note that, in general, $\overline{W}$ is a priori not available and
$\eqref{eq:intro:macro:P}$ is thus a purely formal relation. For reasons of
readability, the dependence of any quantity on the macroscopic material
coordinate $\overline{\mbox{\boldmath$X$}}$ is usually spared, e.g.
$\overline{\mbox{\boldmath$F$}}=\overline{\mbox{\boldmath$F$}}(\overline{\mbox{\boldmath$X$}})$,
$\overline{\mbox{\boldmath$H$}}=\overline{\mbox{\boldmath$H$}}(\overline{\mbox{\boldmath$X$}})$,
or
$\overline{\mbox{\boldmath$u$}}=\overline{\mbox{\boldmath$u$}}(\overline{\mbox{\boldmath$X$}})$.

####  1.6.2 Microscopic problem

The _microscopic boundary value problem_ is given by the balance equations

$\displaystyle{\rm Div}_{\rm X}(\mbox{\boldmath$P$})$
$\displaystyle=\mbox{\boldmath$0$}$ (8)
$\displaystyle\mbox{\boldmath$F$}^{-1}\mbox{\boldmath$P$}$
$\displaystyle=\mbox{\boldmath$P$}^{\sf T}\mbox{\boldmath$F$}^{\sf-T}$ (9)

and suitable boundary conditions, e.g. as discussed in Miehe (2003). In this
work, _periodic displacement fluctuation boundary conditions_ are employed.
The microscopic displacements take the form

$\displaystyle\mbox{\boldmath$u$}(\mbox{\boldmath$X$};\overline{\mbox{\boldmath$F$}})$
$\displaystyle=\overline{\mbox{\boldmath$u$}}+\overline{\mbox{\boldmath$H$}}\mbox{\boldmath$X$}+\mbox{\boldmath$w$}(\mbox{\boldmath$X$};\overline{\mbox{\boldmath$F$}}).$
(10)

Therein, the macroscopic displacement is independent of microscopic
quantities. The second term,
$\overline{\mbox{\boldmath$H$}}\mbox{\boldmath$X$}$, corresponds to a
homogeneous deformation of the microstructure. The third term,
$\mbox{\boldmath$w$}(\mbox{\boldmath$X$})$, is a _displacement fluctuation_
with the _zero mean property_
$\langle\mbox{\boldmath$w$}\rangle=\mbox{\boldmath$0$}$. Hence, the
deformation gradient reads

$\displaystyle\mbox{\boldmath$F$}(\mbox{\boldmath$X$};\overline{\mbox{\boldmath$F$}})$
$\displaystyle=\overline{\mbox{\boldmath$F$}}+\widetilde{\mbox{\boldmath$H$}}(\mbox{\boldmath$X$};\overline{\mbox{\boldmath$F$}})=\overline{\mbox{\boldmath$F$}}+\widetilde{\mbox{\boldmath$F$}}(\mbox{\boldmath$X$};\overline{\mbox{\boldmath$F$}}),$
(11)

where the fluctuational part
$\widetilde{\mbox{\boldmath$H$}}=\mbox{\boldmath$w$}\otimes\nabla_{\rm
X}=\widetilde{\mbox{\boldmath$F$}}$, has the zero mean property

$\displaystyle\left<\widetilde{\mbox{\boldmath$F$}}\right>$
$\displaystyle=\mbox{\boldmath$0$}.$ (12)

Thus, the volume average of the local deformation gradient equals its
macroscopic counterpart,

$\displaystyle\overline{\mbox{\boldmath$F$}}$
$\displaystyle=\left<\mbox{\boldmath$F$}\right>.$ (13)

This motivates the identification of $\overline{\mbox{\boldmath$F$}}$ as _the
boundary condition_ to the microscopic problem (8). As for the material
response, the following relationships can be deduced:

$\displaystyle\overline{W}$ $\displaystyle=\left<W\right>$ (14)
$\displaystyle\overline{\mbox{\boldmath$P$}}$
$\displaystyle=\left<\mbox{\boldmath$P$}\right>$ (15)
$\displaystyle\overline{{\mathbb{C}}}$
$\displaystyle\neq\left<{\mathbb{C}}\right>$ (16)

Equations (13) and (15) are called kinematic and static _coupling relations_ ,
respectively. The inequality (16) generally holds for heterogeneous
microstructures, even in the most simple case of infinitesimal strains and
linear elasticity. More precisely, the volume average overestimates the
effective stiffness in the spectral sense,

$\displaystyle\mbox{\boldmath$x$}\cdot\left<{\mathbb{C}}\right>\cdot\mbox{\boldmath$x$}$
$\displaystyle\geq\mbox{\boldmath$x$}\cdot\overline{{\mathbb{C}}}\cdot\mbox{\boldmath$x$}\qquad\forall\text{
second order tensors }\mbox{\boldmath$x$}.$ (17)

##  2 Reduced Basis homogenization for hyperelasticity

###  2.1 Formulation

The _Reduced Basis (RB)_ scheme is based on a direct approximation of the
microscopic deformation gradient $F$ from Equation (11) without the need to
explicitly have the corresponding displacement at hand. The initial
approximation is given by

$\displaystyle\mbox{\boldmath$F$}_{\xi}(\mbox{\boldmath$X$};\overline{\mbox{\boldmath$F$}},\underline{\xi})$
$\displaystyle=\overline{\mbox{\boldmath$F$}}+\sum_{i=1}^{N}\mbox{\boldmath$B$}^{(i)}(\mbox{\boldmath$X$})\xi_{i}.$
(18)

The parameters $\overline{\mbox{\boldmath$F$}}$ and
$\underline{\xi}\in{\mathbb{R}}^{N}$ are the boundary condition to (8) and the
reduced coefficient vector, respectively. Note that the macroscopic
coordinates $\overline{\mbox{\boldmath$X$}}$ is not assumed to be an explicit
part of the parameter set, i.e. the same approximation is made throughout the
macrostructure. The set $\\{\mbox{\boldmath$B$}^{(i)}\\}_{i=1}^{N}$ is
linearly independent within the space $L^{2}(\varOmega_{0})$ and is called _RB
of $F$_. In a later context, it will also be referred to as the set of _ansatz
functions_. In order to enforce the relationship

$\displaystyle\left<\mbox{\boldmath$F$}_{\xi}\right>$
$\displaystyle=\overline{\mbox{\boldmath$F$}}$ (19)

regardless of the reduced coefficients $\underline{\xi}$, the basis functions
are asserted to have the _fluctuation_ property, i.e. for $i=1,\dots,N$

$\displaystyle\left<\mbox{\boldmath$B$}^{(i)}\right>$
$\displaystyle=\mbox{\boldmath$0$}.$ (20)

For now, the RB is assumed to be given.

The ansatz (18) allows for evaluation of the energy at the macroscale as a
function of the macroscopic kinematic variable
$\overline{\mbox{\boldmath$F$}}$ and of the reduced coefficients
$\underline{\xi}$,

$\displaystyle\overline{W}_{\xi}(\overline{\mbox{\boldmath$F$}},\underline{\xi})$
$\displaystyle=\left<W(\mbox{\boldmath$F$}_{\xi})\right>.$ (21)

By the principle of minimum energy, the optimal coefficients

$\displaystyle\underline{\xi}^{*}(\overline{\mbox{\boldmath$F$}})$
$\displaystyle=\arg\min_{\underline{\xi}\in{\mathbb{R}}^{N}}\overline{W}_{\xi}(\overline{\mbox{\boldmath$F$}},\underline{\xi})$
(22)

are sought-after. The unconstrained and unique solvability of this task is
assumed for the moment and will be addressed in Section 2.4. The solution of
(22) defines the _RB approximation_ of the deformation gradient

$\displaystyle\mbox{\boldmath$F$}^{\rm
RB}(\mbox{\boldmath$X$};\overline{\mbox{\boldmath$F$}})$
$\displaystyle=\overline{\mbox{\boldmath$F$}}+\sum_{i=1}^{N}\mbox{\boldmath$B$}^{(i)}(\mbox{\boldmath$X$})\xi_{i}^{*}(\overline{\mbox{\boldmath$F$}}).$
(23)

The microscopic energy, stress and stiffness within the microstructure are
then approximated by

$\displaystyle W^{\rm RB}(\mbox{\boldmath$X$};\overline{\mbox{\boldmath$F$}})$
$\displaystyle=W(\mbox{\boldmath$F$}^{\rm
RB}(\mbox{\boldmath$X$};\overline{\mbox{\boldmath$F$}})),$ (24)
$\displaystyle\mbox{\boldmath$P$}^{\rm
RB}(\mbox{\boldmath$X$};\overline{\mbox{\boldmath$F$}})$
$\displaystyle=\displaystyle\frac{\displaystyle\partial
W}{\displaystyle\partial\mbox{\boldmath$F$}}(\mbox{\boldmath$F$}^{\rm
RB}(\mbox{\boldmath$X$};\overline{\mbox{\boldmath$F$}})),$ (25)
$\displaystyle{\mathbb{C}}^{\rm
RB}(\mbox{\boldmath$X$};\overline{\mbox{\boldmath$F$}})$
$\displaystyle=\displaystyle\frac{\displaystyle\partial^{2}W}{\displaystyle\partial{\mbox{\boldmath$F$}}^{2}}(\mbox{\boldmath$F$}^{\rm
RB}(\mbox{\boldmath$X$};\overline{\mbox{\boldmath$F$}})),$ (26)

respectively. The resulting effective energy is readily given by

$\displaystyle\overline{W}^{\rm RB}(\overline{\mbox{\boldmath$F$}})$
$\displaystyle=\left<W^{\rm RB}\right>.$ (27)

Also, the effective responses $\overline{\mbox{\boldmath$P$}}^{\rm
RB}(\overline{\mbox{\boldmath$F$}})$ and $\overline{{\mathbb{C}}}^{\rm
RB}(\overline{\mbox{\boldmath$F$}})$ may now be calculated. However, before
going into detail on that, it is advantageous to first elaborate on the
solution of the minimization problem (22). This short survey will reveal
essential properties of some occurring quantities which are important for the
determination of the effective material response.

The necessary, first order optimality conditions to (22) define the components
of the residual vector $\underline{r}\in{\mathbb{R}}^{N}$,

$\displaystyle r_{i}(\overline{\mbox{\boldmath$F$}},\underline{\xi})$
$\displaystyle=\displaystyle\frac{\displaystyle\partial\overline{W}_{\xi}}{\displaystyle\partial\xi_{i}}(\overline{\mbox{\boldmath$F$}},\underline{\xi})=\left<\displaystyle\frac{\displaystyle\partial
W}{\displaystyle\partial\mbox{\boldmath$F$}}(\mbox{\boldmath$F$}_{\xi})\cdot\displaystyle\frac{\displaystyle\partial\mbox{\boldmath$F$}_{\xi}}{\displaystyle\partial\xi_{i}}\right>=\left<\mbox{\boldmath$P$}(\mbox{\boldmath$F$}_{\xi})\cdot\mbox{\boldmath$B$}^{(i)}\right>\overset{!}{=}0\qquad(i=1,\dots,N).$
(28)

######  Note 1.

The solution stress field $\mbox{\boldmath$P$}^{\rm RB}$ is
$L^{2}(\varOmega_{0})$-orthogonal to the RB ansatz functions
$\\{\mbox{\boldmath$B$}^{(i)}\\}_{i=1}^{N}$.

A viable choice for the solution of the minimization problem (22) is the
Newton-Raphson scheme, which necessitates the Jacobian
$\underline{\underline{D}}\in{\rm Sym}({\mathbb{R}}^{N\times N})$ with the
components

$\displaystyle D_{ij}$
$\displaystyle=\left<\mbox{\boldmath$B$}^{(i)}\cdot\displaystyle\frac{\displaystyle\partial\mbox{\boldmath$P$}}{\displaystyle\partial\xi_{j}}(\mbox{\boldmath$F$}_{\xi})\right>=\left<\mbox{\boldmath$B$}^{(i)}\cdot{\mathbb{C}}(\mbox{\boldmath$F$}_{\xi})\cdot\mbox{\boldmath$B$}^{(j)}\right>=D_{ji}\qquad(i,j=1,\dots,N)\,.$
(29)

Then, the $k$th iteration to the solution
$\underline{\xi}^{*}(\overline{\mbox{\boldmath$F$}})$ reads

$\displaystyle\underline{\xi}^{(k)}$
$\displaystyle=\underline{\xi}^{(k-1)}-\left(\underline{\underline{D}}^{(k-1)}\right)^{-1}\underline{r}^{(k-1)}\qquad(k\geq
1)\,.$ (30)

The initial guess $\underline{\xi}^{(0)}$ can be zero or a more sophisticated
alternative.

The deduction of the effective material response by means of Note 1 and the
Jacobian $\underline{\underline{D}}$ is given in Appendix B. The following
list summarizes the homogenized quantities arising from the $F$-RB:

$\displaystyle\overline{\mbox{\boldmath$F$}}$
$\displaystyle=\left<\mbox{\boldmath$F$}^{\rm RB}\right>$
$\displaystyle\text{see }\eqref{eq:RB:upscaling}$
$\displaystyle\overline{W}^{\rm RB}$ $\displaystyle=\left<W^{\rm RB}\right>$
$\displaystyle\text{see }\eqref{eq:RB:energy}$
$\displaystyle\overline{\mbox{\boldmath$P$}}^{\rm RB}$
$\displaystyle=\left<\mbox{\boldmath$P$}^{\rm RB}\right>$ see Appendix B.1
$\displaystyle\overline{{\mathbb{C}}}^{\rm RB}$
$\displaystyle=\left<{\mathbb{C}}^{\rm
RB}\right>-\sum_{i,j=1}^{N}D_{ij}^{-1}\left<{\mathbb{C}}^{\rm
RB}\cdot\mbox{\boldmath$B$}^{(i)}\right>\otimes\left<\mbox{\boldmath$B$}^{(j)}\cdot{\mathbb{C}}^{\rm
RB}\right>$ see Appendix B.2

In total, the problem of determining both the local field $F$ and the
homogenized material properties depends only on $N$ degrees of freedom, namely
on the $N$ coefficients $\xi_{i}$. Usually, the number of DOF $N$ is in the
range of 10 to 150, which compares to the full order model’s complexity that
can easily reach more than 105 DOF.

Despite this impressive reduction of the number of DOF, the computational
costs associated with the assembly of the residual $\underline{r}$ and of the
Jacobian $\underline{\underline{D}}$ still relate to the number of quadrature
points of the microstructural discretization.

This method extends corresponding methods for the geometrically linear case,
where the infinitesimal strain tensor $\mbox{\boldmath$\varepsilon$}={\rm
sym}(\mbox{\boldmath$H$})$ is considered. For more information on these
topics, the reader is referred to the authors’ previous work Fritzen and Kunc
(2018) or standard literature, such as (Castañeda and Suquet, 1998, part
II.C).

###  2.2 Identification of the Reduced Basis

The basis $\\{\mbox{\boldmath$B$}^{(i)}\\}_{i=1}^{N}$ is computed by means of
a classical snapshot POD, see (Sirovich, 1987). In contrast to many other POD
based reduction methods, it is important to point out that here, the primal
variable is _not_ taken to be the displacement field. Instead, the POD is
performed on fluctuations $\widetilde{\mbox{\boldmath$F$}}$ of the deformation
gradient.

During the snapshot POD, first, _snapshots_ are created by means of high-
fidelity solutions to the nonlinear microscopic problem (8) with different
_snapshot boundary conditions_ $\overline{\mbox{\boldmath$F$}}^{(j)}$,
$j=1,\dots,N_{\rm s}$, which are also called _training boundary conditions_.
Each of these boundary conditions leads to a solution field
$\mbox{\boldmath$F$}^{(j)}(\mbox{\boldmath$X$};\overline{\mbox{\boldmath$F$}})$.
Typically, the Finite Element Method (FEM) or solvers making use of the Fast
Fourier Transform (e.g. Kabel et al., 2014) are used to this end. The RB
method presented here is independent of the discretization method utilized to
obtain full field solutions. It is merely necessary to know the quadrature
weights and the related discrete values of the solutions
$\mbox{\boldmath$F$}^{(j)}(\mbox{\boldmath$X$};\overline{\mbox{\boldmath$F$}}^{(j)})$.
For better readability, the continuous notation is maintained for the moment.
The corresponding fluctuation fields are computed by means of local
subtraction of the macroscopic deformation gradient,

$\displaystyle\widetilde{\mbox{\boldmath$F$}}^{(j)}(\mbox{\boldmath$X$};\overline{\mbox{\boldmath$F$}}^{(j)})$
$\displaystyle=\mbox{\boldmath$F$}^{(j)}(\mbox{\boldmath$X$};\overline{\mbox{\boldmath$F$}}^{(j)})-\overline{\mbox{\boldmath$F$}}^{(j)}\qquad(j=1,\dots,N_{\rm
s}).$ (31)

Each of these $N_{\rm s}$ fluctuation fields
$\widetilde{\mbox{\boldmath$F$}}^{(j)}$ represents one snapshot.

Second, the most dominant features of the snapshots are extracted. This is
done by means of the eigendecomposition of the correlation matrix. It consists
of the mutual $L^{2}(\varOmega_{0})$ scalar products of the snapshots,
$\left<\widetilde{\mbox{\boldmath$F$}}^{(i)}\cdot\widetilde{\mbox{\boldmath$F$}}^{(j)}\right>$
($i,j=1,\dots,N_{\rm s}$). The remaining procedure is standard, see for
instance (Sirovich, 1987) or (Quarteroni et al., 2016): the $N_{\rm s}$
eigenvalues $\lambda_{j}$, corresponding to the eigenvectors
$\underline{E}^{(j)}\in{\mathbb{R}}^{N_{\rm s}}$, are sorted in descending
order and truncated by only considering the first $N$ values,
$\lambda_{1}\geq\ldots\geq\lambda_{N}$. The decision on a particular threshold
index $N$ is based on consideration of the Schmidt-Eckhard-Young-Mirsky
theorem. Finally, the RB is constructed as

$\displaystyle\mbox{\boldmath$B$}^{(i)}(\mbox{\boldmath$X$})$
$\displaystyle=\sum_{j=1}^{N_{\rm
s}}\frac{1}{\sqrt{\lambda_{i}}}\left(\underline{E}^{(i)}\right)_{j}\widetilde{\mbox{\boldmath$F$}}^{(j)}(\mbox{\boldmath$X$})\quad(i=1,\ldots,N)$
(32)

where the factor $1/\sqrt{\lambda_{i}}$ accounts for $L^{2}(\varOmega_{0})$
normalization of the basis elements. We conclude that the RB is a collection
of $L^{2}(\varOmega_{0})$ orthonormal $\widetilde{\mbox{\boldmath$H$}}$-like
quantities.

###  2.3 Mathematical motivation of the Reduced Basis model

Next, the obtained deformation gradient field $\mbox{\boldmath$F$}^{\rm RB}$
and the associated stress field $\mbox{\boldmath$P$}^{\rm RB}$ are shown to
weakly solve the original problem (8), ${\rm Div}_{\rm
X}(\mbox{\boldmath$P$})=\mbox{\boldmath$0$}$.

Let $\delta\mbox{\boldmath$w$}\in H^{1}_{0}(\varOmega_{0})$ be an admissible
test function, i.e. a once weakly differentiable periodic displacement
fluctuation field, and let
$\delta\mbox{\boldmath$H$}=\delta\mbox{\boldmath$w$}\otimes\nabla_{\rm X}$
denote its gradient. Then the well-known weak form of (8) is equivalent to the
principle of virtual work,

$\displaystyle\delta\overline{W}=$
$\displaystyle\int_{\varOmega_{0}}\mbox{\boldmath$P$}\cdot\delta\mbox{\boldmath$H$}\,{\rm
d}V=0.$ (33)

The residual $\underline{r}$ from (28) coincides with the left-hand side of
the weak form, if the test function $\delta\mbox{\boldmath$w$}$ is chosen
suitably. In fact, for each basis element $\mbox{\boldmath$B$}^{(i)}$ there is
a corresponding displacement fluctuation field $\delta\mbox{\boldmath$w$}$
that is defined uniquely up to rigid body motion, hence the equivalence of
(28) and (33).

It follows that the Reduced Basis scheme is a _Galerkin_ procedure, taking the
displacement fields corresponding to the RB elements
$\mbox{\boldmath$B$}^{(i)}$ as both ansatz and test functions. Hence, the _RB
model is equivalent to the FEM_ , but with basis functions that are globally
supported in $\varOmega_{0}\backslash\partial\varOmega_{0}$. In other words,
the basis functions of the RB method span subspace of dimension $N$ within the
high-dimensional space of FE basis functions. Although this property is shared
with RB schemes based on POD of displacement snapshots, a notable difference
is that this novel approach directly operates on fields entering the
constitutive equations.

###  2.4 Details on the coefficient optimization

The coefficient optimization task (22) leads to a weak solution of the
microscopic boundary value problem, as was just shown. Hence, the well-
established theories on which the FEM is built, e.g. the calculus of
variations, are applicable to the presented method just as much. This implies
that the well-known issues with suitable convexity conditions and with
existence and uniqueness of minimizers apply to the RB method, too. We focus
on ad hoc numerical treatments of these issues. For more details on the
theoretical part, the reader is referred to standard literature, such as Ball
(1976), and recent developments in this matter, e.g. Schneider (2017).

A constraint to the optimization problem is the physical condition

$\displaystyle J$ $\displaystyle={\rm
det}(\mbox{\boldmath$F$}_{\xi}(\mbox{\boldmath$X$}))>0\quad\forall\,\mbox{\boldmath$X$}\in\varOmega_{0},$
(34)

which means that no singular ($J=0$) or self-penetrating ($J<0$) deformations
shall occur. This reduces the set of admissible coefficients $\underline{\xi}$
to a subset of ${\mathbb{R}}^{N}$ that is non-trivial to access. The
positiveness of the microscopic determinant of the deformation gradient
introduces a constraint to the, thus far, unconstrained minimum problem (22)
representing the weak form of (8) in the RB setting.

In case of a violation of the inequality (34), the implementation of the RB
method is prone to failure as soon as the constitutive law (4) is evaluated.
This occurs only in the context of volume averaging, i.e. for the assembly of
the residual, the Jacobian, or the effective energy, stress, or stiffness. The
numerical quadrature approximating the volume averaging operation is

$\displaystyle\left<\bullet\right>$
$\displaystyle=\frac{1}{|\varOmega_{0}|}\int_{\varOmega_{0}}\bullet(\mbox{\boldmath$X$})\,{\rm
d}V\approx\frac{1}{|\varOmega_{0}|}\sum_{p=1}^{N_{\rm
qp}}\bullet(\mbox{\boldmath$X$}_{p})\,w_{p}.$ (35)

Here, $N_{\rm qp}$, $\mbox{\boldmath$X$}_{p}$, and $w_{p}$ respectively denote
the number of quadrature points, their positions and the corresponding
positive weights. Moreover, even if the inequality (34) is satisfied
everywhere, the local field $\mbox{\boldmath$F$}_{\xi}$ might possibly have
some positive but overly small values of the determinant, $0<{\rm
det}(\mbox{\boldmath$F$}_{\xi}(\mbox{\boldmath$X$}))\ll 1$, that are
unphysical. In such a case, the energy optimization, cf. (22), would be
dominated by these nearly singular points. Instead of allowing the
optimization to focus almost exclusively on these exceptional points, we
interpret unphysically small values of the determinant as limitations to the
reliability of the RB method. On the other hand, large values ${\rm
det}(\mbox{\boldmath$F$}_{\xi}(\mbox{\boldmath$X$}))\gg 1$ are not too
detrimental to the functionality of the scheme, although such values are just
as questionable.

Thus, the following _weighted numerical quadrature rule_ is introduced:

$\displaystyle\left<\bullet\right>$
$\displaystyle\approx\left(\sum_{p=1}^{N_{\rm
qp}}\bullet(\mbox{\boldmath$X$}_{p})\,\phi(J_{p})\,w_{p}\right)/\left(\sum_{p=1}^{N_{\rm
qp}}\phi(J_{p})w_{p}\right)$ (36)

Therein, the almost smooth _cutoff function
$\phi:{\mathbb{R}}\to{\mathbb{R}}_{\geq 0}$_ is empirically defined by

$\displaystyle\phi(J)$ $\displaystyle=\begin{cases}1&\text{if $J>0.6$}\\\
0.5\,{\rm erf}(30\,J-15)+0.5&\text{if $0.4<J\leq 0.6$}\\\ 0&\text{if $J\leq
0.4$}\end{cases}\,\text{,}$ (37)

which makes use of the well-known error function. The cutoff function is
depicted in Figure 2. This reliability indicator could in principle be
modified, e.g. the steepness, the smoothness, and its center could be
adjusted. Thus, it should be regarded as an example only.

Figure 2: Cutoff function $\phi$. Its value is used as a reliability factor in
the numerical quadrature.

This weighted numerical quadrature rule is used henceforward for all
_numerical_ volume averaging operations. Its application will not be noted
explicitly. However, the _theoretical_ derivation of the RB method as
described in Section 2.1 is not affected by this, i.e. volume averaging
operations remain exact as far as the theory is concerned. The two most
important consequences of this numerical tweak are:

  * •

The RB method is robust with respect to outlier values of the determinant. The
modified quadrature rule extends the set of coefficient vectors
$\underline{\xi}$ for which effective quantities can be computed, albeit
approximately, to the whole space ${\mathbb{R}}^{N}$.

  * •

The significance of local fields varies with the value of the cutoff function.
When $\phi$ attains values less than one, information is considered
accordingly less reliable. In this sense, microscopic information is filtered
based on a trust region for $J$ defined by $\phi$ can be seen as a reliability
indicator.

##  3 Sampling

###  3.1 General considerations

The proposed sampling strategy builds on the previous work Fritzen and Kunc
(2018), in which the authors proposed an analogous sampling procedure for the
small strain setting. However, substantial modifications are required in order
to account for the finite rotational part, $\overline{\mbox{\boldmath$R$}}$,
of the macroscopic deformation gradient, $\overline{\mbox{\boldmath$F$}}$, and
the nonlinearity of the volumetric part of the deformation, $J$, with respect
to the local displacements, $u$. For the setup of the Reduced Basis model as
described in Section 2.2, the space of macroscopic deformation gradients

$\displaystyle\overline{{\cal F}}$ $\displaystyle=\\{\text{second order
tensors
}\overline{\mbox{\boldmath$F$}}\,|\,\det{\overline{\mbox{\boldmath$F$}}}>0\\}$
(38)

needs to be sampled, i.e. the discrete sampling set $\overline{{\cal F}}_{\rm
s}=\\{\overline{\mbox{\boldmath$F$}}^{(m)}\\}_{m=1}^{N_{\rm
s}}\subset\overline{{\cal F}}$ is to be defined. Two contradictory
requirements need to be satisfied when constructing $\overline{{\cal F}}_{\rm
s}$:

  1. 1.

The samples should be _densely and homogeneously distributed_ within the space
of all admissible macroscopic kinematic configurations. This is owing to the
desire that the POD may extract correlation information from a wholistic and
unbiased set. In other words, the samples should be as _uniformly random_ as
possible _within the anticipated query domain of the surrogate_.

  2. 2.

The sample number $N_{\rm s}$ should not exceed a certain limit. Only with
this property may the RB be identified within the bounds of available
computational resources (e.g. memory and CPU time).

###  3.2 Large strain sampling strategy

The set of admissible macroscopic deformation gradients $\overline{{\cal F}}$
is a subset of a nine-dimensional space
($\overline{\underline{\underline{F}}}\in{\mathbb{R}}^{3\times
3}\sim{\mathbb{R}}^{9}$) restricted by the inequality

$\displaystyle\det{\overline{\mbox{\boldmath$F$}}}>0.$ (39)

Therefore, a regular grid in the components of
$\overline{\mbox{\boldmath$F$}}$ might lead to a prohibitively large amount of
samples, and even to a violation of (39). For instance, such a grid with a
rather moderate resolution of just 10 samples of each component would require
1 billion solutiotns of the full order model. Also, the subsequent POD would
involve a snapshot matrix and/or correlation matrix of accordingly huge
dimensionality.

In order to decrease the dimension of the sampling space, recall the polar
decomposition of the deformation gradient,
$\overline{\mbox{\boldmath$F$}}=\overline{\mbox{\boldmath$R$}}\,\overline{\mbox{\boldmath$U$}}$,
where $\overline{\mbox{\boldmath$R$}}$ is a rotation and
$\overline{\mbox{\boldmath$U$}}$ is the symmetric positive definite (s.p.d.)
stretch tensor. Material objectivity implies the energy function to be
independent of the frame of reference,

$\displaystyle\overline{W}(\overline{\mbox{\boldmath$R$}}\,\overline{\mbox{\boldmath$U$}})$
$\displaystyle=\overline{W}(\overline{\mbox{\boldmath$U$}}),$ (40) and the
transformation behavior
$\displaystyle\overline{\mbox{\boldmath$P$}}(\overline{\mbox{\boldmath$R$}}\,\overline{\mbox{\boldmath$U$}})$
$\displaystyle=\overline{\mbox{\boldmath$R$}}\,\overline{\mbox{\boldmath$P$}}(\overline{\mbox{\boldmath$U$}})$
(41) for the first Piola-Kirchhoff stress and
$\displaystyle\overline{C}_{ijkl}(\overline{\mbox{\boldmath$R$}}\,\overline{\mbox{\boldmath$U$}})$
$\displaystyle=\sum_{m,n=1}^{3}\delta_{ik}\overline{U}^{-1}_{lm}\overline{P}_{mj}(\overline{\mbox{\boldmath$U$}})+\overline{R}_{im}\overline{C}_{mjnl}(\overline{\mbox{\boldmath$U$}})\overline{R}_{kn}\qquad(i,j,k,l=1,2,3)$
(42)

for the components of the corresponding stiffness tensor ${\mathbb{C}}$, see
Appendix A. These known facts lead to the following

######  Note 2.

In order to collect representative samples of the hyperelastic response
functions $\overline{W}$, $\overline{\mbox{\boldmath$P$}}$, and
$\overline{{\mathbb{C}}}$, it suffices to evaluate them on samples of the
stretch tensor $\overline{\mbox{\boldmath$U$}}\sim{\mathbb{R}}^{6}$ instead of
evaluating them on samples of the deformation gradient
$\overline{\mbox{\boldmath$F$}}\sim{\mathbb{R}}^{9}$.

This effectively reduces the number of dimensions of $\overline{{\cal F}}$
from nine to six. The same dimensionality is attained when considering the
response functions with respect to a symmetric measure of strain, e.g. as is
done in Miehe (1996) where the tangent stiffness is efficiently computed using
a perturbation technique. However, such measures are unsuitable for reduction
by means of the proposed RB method.

The remaining six-dimensional space of s.p.d. tensors is not linear but a
convex cone (which does not include the zero element). Moreover, linearly
combining elements within this space quickly leads to values of
$\overline{J}=\det{\overline{\mbox{\boldmath$U$}}}$ describing unphysically
large changes in volume. For instance,
$\overline{\mbox{\boldmath$U$}}=1.3\cdot\mbox{\boldmath$I$}$, with the second-
order identity tensor $I$, equates to more than 100% volumetric extension,
which is well beyond the regime of usual hyperelastic materials that are often
nearly incompressible. On the other hand, 100% deviatoric strain is within the
range of many standard materials, such as rubber. Hence, _in order to describe
the space of practically relevant stretch tensors_ , we propose to apply the
DDMS to the macroscopic stretch tensor,

$\displaystyle\overline{\mbox{\boldmath$U$}}$
$\displaystyle=\overline{J}^{1/3}\,\widehat{\overline{\mbox{\boldmath$U$}}}.$
(43)

Let $\widehat{\overline{{\cal U}}}$ denote the manifold of unimodular
macroscopic stretch tensors
$\widehat{\overline{\mbox{\boldmath$U$}}}=(\overline{J})^{-1/3}\overline{\mbox{\boldmath$U$}}$.
The multiplicative split (43) is the basis for The set of practically relevant
macroscopic stretch tensors $\overline{{\cal U}}$ can be sampled via sampling
of both the macroscopic determinants,

$\displaystyle\left\\{\overline{J}^{(m)}\right\\}_{m=1}^{N_{\rm dil}}$
$\displaystyle\subset{\mathbb{R}}_{+}\text{ ,}$ and the macroscopic unimodular
stretch tensors,
$\displaystyle\left\\{\widehat{\overline{\mbox{\boldmath$U$}}}^{(j)}\right\\}_{j=1}^{N_{\rm
dev}}$ $\displaystyle\subset\widehat{\overline{{\cal U}}}.$

The sampling set is determined by the product set

$\displaystyle\left\\{\left(\overline{J}^{(m)}\right)^{1/3}\widehat{\overline{\mbox{\boldmath$U$}}}^{(j)}\right\\}_{m,j=1}^{m=N_{\rm
dil},\,j=N_{\rm dev}}\subset\overline{{\cal U}}.$ (44)

The choice of the dilatational samples is relatively straight-forward. For
many common materials, the expected range of $\overline{J}$ is rather narrow,
so a few equisized or adaptive sub-intervals around $\overline{J}=1$ deliver
sufficient resolution.

For the space of unimodular s.p.d. matrices representing
$\widehat{\overline{\mbox{\boldmath$U$}}}\in\widehat{\overline{{\cal U}}}$,
basic results of Lie group theory can be exploited. We restrict to stating
well-known facts that are necessary at this point. For more details, the
interested reader is referred to standard text books, such as Faraut (2008).

######  Corollary .

Let

$\displaystyle\rm{symsl}$
$\displaystyle=\left\\{\underline{\underline{X}}\in{\mathbb{R}}^{3\times
3}\,\big{|}\,\underline{\underline{X}}=\underline{\underline{X}}^{\sf T},{\rm
tr}(\underline{\underline{X}})=0\right\\}$ be the _tangent space_ and
$\displaystyle\rm{SymSL}_{+}$
$\displaystyle=\Big{\\{}\underline{\underline{U}}\in{\mathbb{R}}^{3\times
3}\,\big{|}\,\underline{\underline{U}}=\underline{\underline{U}}^{\sf
T},\break\det{\underline{\underline{U}}}=1,\underline{x}^{\sf
T}\,\underline{\underline{U}}\,\underline{x}>0\,\forall\underline{x}\in{\mathbb{R}}^{3}\Big{\\}}$

be the _manifold_ of unimodular s.p.d. matrices. The matrix exponential maps
the tangent space bijectively onto the manifold,

$\displaystyle\rm{exp}:\quad\rm{symsl}\to\rm{SymSL}_{+}\,.$

The proof of this statement is standard, e.g. by means of the eigenvalue
decomposition, and does not necessitate the reference to the abstract setting
of Lie groups. In fact, all of the following arguments could be given without
the notion of tangent spaces and manifolds. However, this notion is a
fundamental concept in nonlinear mechanics. For a very descriptive and
comprehensive work on this topic, the reader is refered to Neff et al. (2016).
We choose to build upon this concept, as it comes along with vivid
interpretations of the function spaces $\overline{{\cal U}}$ and
$\widehat{\overline{U}}$.

The set $\rm SymSL_{+}$ is the set of matrix representations of
$\widehat{\overline{\mbox{\boldmath$U$}}}\in\widehat{\overline{{\cal U}}}$.
Notably, the tangent space $\rm{symsl}$ is linear. Hence, by virtue of the
matrix exponential, the sampling of the nonlinear manifild
$\widehat{\overline{{\cal U}}}$ can be reduced to the sampling of a linear
space. Moreover, the restrictions of symmetry and zero trace render the
tangent space five-dimensional. This property is by definition shared with the
manifold $\rm SymSL_{+}$.

The two-dimensional case is now addressed for the sake of visualization. In
this setting, the nonlinearity of the manifold and the structure of the
sampling set $\overline{{\cal U}}$ can be illustrated figuratively. With the
subscript (2) denoting two-dimensional quantities, a basis of the tangent
space is given by

$\displaystyle\underline{\underline{X}}^{(1)}_{(2)}$
$\displaystyle=\sqrt{\frac{1}{2}}\begin{bmatrix}1&0\\\
0&-1\end{bmatrix}\;\text{ and }$
$\displaystyle\underline{\underline{X}}^{(2)}_{(2)}$
$\displaystyle=\sqrt{\frac{1}{2}}\begin{bmatrix}0&1\\\ 1&0\end{bmatrix}.$

The stretch tensors are obtained through

$\displaystyle\underline{\underline{\overline{U}}}_{(2)}$
$\displaystyle=\overline{J}^{\frac{1}{2}}{\rm
exp}\left(t\left(\alpha\underline{\underline{X}}^{(1)}_{(2)}+\beta\underline{\underline{X}}^{(2)}_{(2)}\right)\right)=\begin{bmatrix}a&b\\\
b&d\end{bmatrix}.$

A visualization of such samples is depicted in Figure 3. There, for the sake
of visual clarity, the determinant $\overline{J}$ is sampled by four
equidistant (and rather unrealistic) values between 0.1 and 4. The value
$t\in[-2,2]$ is called _deviatoric amplitude_. A densely uniform sampling
$\varphi\in[0,2\pi)$ yields the coefficients $\alpha=\cos{\varphi}$ and
$\beta=\sin{\varphi}$.

Figure 3: Visualizations of the family of
$\overline{\mbox{\boldmath$U$}}$-manifolds with constant determinant
$\overline{J}\in\\{0.1,\ldots,4.0\\}$ for the two-dimensional case from two
different perspectives. The green line represents the set
$\\{\lambda\mbox{\boldmath$I$}\,|\,\lambda>0\\}$.

This emphasizes the important role of the DDMS in the context of sampling, as
utilized in (44): it allows for the identification of a physically meaningful
sampling domain that is much smaller than the surrounding space of all
admissible stretch tensors. On a side note, the isodet surfaces are
perpendicular to the line representing the dilatational stretch tensors.

The proposed sampling procedure for $\overline{{\cal U}}$ in three dimensions
is given in Algorithm 1. For this purpose, an orthonormal basis
$\underline{\underline{X}}^{(1)},\dots,\underline{\underline{X}}^{(5)}$ of the
tangent space $\rm symsl$ is fixed, cf. Appendix C. The numbers of different
kind of samples $N_{\rm det}$, $N_{\rm dir}$, and $N_{\rm amp}$ relate to the
quantities $N_{\rm dil}$ and $N_{\rm dev}$ introduced in (44) by $N_{\rm
det}=N_{\rm dil}$ and $N_{\rm dir}N_{\rm amp}=N_{\rm dev}$.

1 $\overline{J}_{\rm min},\overline{J}_{\rm max}$ minimum and maximum
determinant with $\overline{J}_{\rm min}<1<\overline{J}_{\rm max}$

2$t_{\rm max}>0$ maximum deviatoric amplitude

3$N_{\rm det}$ number of macroscopic determinants

4$N_{\rm dir}$ number of deviatoric directions

5$N_{\rm amp}$ number of deviatoric amplitudes $N_{\rm det}N_{\rm dir}N_{\rm
amp}$ samples of $\overline{\mbox{\boldmath$U$}}$ Place $N_{\rm det}$
determinants regularly between the extremal values,
$\displaystyle\overline{J}_{\rm
min}=\overline{J}^{(1)}<\ldots<1<\ldots<\overline{J}^{\left(N_{\rm
det}\right)}=\overline{J}_{\rm max}\,.$

6 Generate any approximately uniform distribution of $N_{\rm dir}$ directions
in ${\mathbb{R}}^{5}$, e.g. as in Fritzen and Kunc (2018),
$\displaystyle\left\\{\underline{N}^{(n)}\right\\}_{n=1}^{N_{\rm
dir}}\subset\left\\{\underline{N}\in{\mathbb{R}}^{5}\,:\,\lVert\underline{N}\rVert=1\right\\}.$

7 Place $N_{\rm amp}$ deviatoric amplitudes regularly between 0 and the
expected maximum value,  $\displaystyle 0\leq t^{(1)}<\ldots<t^{\left(N_{\rm
amp}\right)}=t_{\rm max}\,.$

Return the set of samples of $\overline{\mbox{\boldmath$U$}}$:
$\displaystyle\left\\{\left(\overline{J}^{(m)}\right)^{1/3}{\rm
exp}\left(t^{(p)}\left[\sum_{k=1}^{5}\left(\underline{N}^{(n)}\right)_{k}\underline{\underline{X}}^{(k)}\right]\right)\right\\}_{m,n,p=1}^{m=N_{\rm
det},\,n=N_{\rm dir},\,p=N_{\rm amp}}\subset\overline{{\cal U}}$ (45)

Algorithm 1 Sampling of the macroscopic stretch tensor.

The order of Steps 1 to 3 is interchangeable. Details on these parts are now
presented:

  1. Step 1.

Uniform seeding of the determinants is actually not required, but any pattern
implying the sampling determinants $\\{\overline{J}^{(k)}\\}_{k=1}^{N_{\rm
det}}$ to be dense in $[\overline{J}_{\rm min},\overline{J}_{\rm max}]$ as
$N_{\rm det}\to\infty$ works without loss of generality. In this way, the
dilatational response may be resolved adaptively.

  2. Step 2.

The generation of uniform point distributions on spheres is a research topic
on its own, see (Brauchart and Grabner, 2015) for an overview. The method
described in Fritzen and Kunc (2018) is based on energy minimization, which is
also used in the present work. Some point sets of various sizes are included
in the example program (Kunc, 2019). More detailed investigations on this
topic and an open-source code of a point generation program are part of
another work, (Kunc and Fritzen, 2018). Alternatively, Equal Area Points
(Leopardi, 2006) may be used as a rough but quickly computable approximation
of such point sets.

  3. Step 3.

As in Step 1, the uniform placement of the deviatoric amplitudes, $t^{(p)}$,
may be substituted by adaptive alternatives. In (Fritzen and Kunc, 2018), we
have suggested to use an exponential distance function.

The result of Steps 2 and 3, i.e. the sampling of the tangent space, is
exemplary visualized in Figure 4 for the two-dimensional case
($\mbox{\boldmath$u$}\in{\mathbb{R}}^{2}$) and for $N_{\rm dir}=5$ and $N_{\rm
amp}=3$. There, an adaptive spacing of the deviatoric amplitudes is applied.
This might be beneficial for capturing strongly changing material behavior
near the relaxed state.

Figure 4: Example of a discretization of the two-dimensional tangent space.
The samples are placed along the equidistant (in higher dimensions: uniformly
distributed) directions with non-uniformly increasing amplitude.

In general, the vector $\underline{N}\in{\mathbb{R}}^{5}$,
$\lVert\underline{N}\rVert=1$, corresponding to a macroscopic stretch tensor
characterizes the _direction of the applied stretch
$\overline{\mbox{\boldmath$U$}}$_, which we also coin the _load case_.

###  3.3 Application of the stretch tensor trained Reduced Basis model

Since the RB model is trained on only the rotationally invariant part of
$\overline{\mbox{\boldmath$F$}}$ but should be applied to general deformation
gradients, the transformation rules (41) and (42) are incorporated during the
evaluation of the surrogate model. Details on the procedure are given in
Algorithm 2.

1 $\overline{\mbox{\boldmath$F$}}$ macroscopic deformation gradient
$\overline{\mbox{\boldmath$P$}}^{\rm RB}(\overline{\mbox{\boldmath$F$}})$,
$\overline{{\mathbb{C}}}^{\rm RB}(\overline{\mbox{\boldmath$F$}})$ effective
material response Compute polar decomposition
$\overline{\mbox{\boldmath$F$}}=\overline{\mbox{\boldmath$R$}}\,\overline{\mbox{\boldmath$U$}}$.

2 Compute approximations of effective stress
$\overline{\mbox{\boldmath$P$}}^{\rm RB}(\overline{\mbox{\boldmath$U$}})$ and
effective stiffness $\overline{{\mathbb{C}}}_{F}^{\rm
RB}(\overline{\mbox{\boldmath$U$}})$.

Transform effective responses to correct frame
$\overline{\mbox{\boldmath$P$}}^{\rm RB}(\overline{\mbox{\boldmath$F$}})$,
$\overline{{\mathbb{C}}}^{\rm RB}(\overline{\mbox{\boldmath$F$}})$, using
$\overline{\mbox{\boldmath$R$}}$, cf. (41), (42).

Algorithm 2 Online phase of the stretch tensor trained Reduced Basis method

##  4 Numerical examples

###  4.1 Reduced Basis for a fibrous microstructure

The applicability of the proposed RB method in combination with the sampling
scheme is now numerically studied for a fibrous microstructure roughly
resembling polymers with woven reinforcements. The goal is to prove the
efficiency of the $F$-RB scheme in principle and under “worst-case”
conditions, the latter meaning that the phase contrast is chosen rather large.
At this stage, yet, it is explicitly not aspired to provide fully realistic
examples. These would require investigations on the proper size of the
microstructure and should employ dissipative material laws, both of which are
not within the scope of this work.

A cubical microstructure with two fibrous inclusions is considered, see Figure
5  _(a)_ and cf. (Kim and Swan, 2003) for a related example. The inclusions
are periodic and occupy approximately 0.7% of the volume. The mesh contains
35516 nodes in 25633 quadratic TET10 elements (C3D10).

*(a)[.25][c]*(b)[.25][c] *(c)[.25][c]*(d)[.25][c]

Figure 5: _(a)_ fibrous microstructure. _(b)-(d):_ random $F$-RB elements

For the matrix, the material constants are chosen to be ${K_{\rm m}=400}$ MPa
and ${G_{\rm m}=0.4}$ MPa, resembling rubber-like material properties. For the
fibres, the values are ${K_{\rm f}=800}$ MPa and ${G_{\rm f}=240}$ MPa. The
latter parameters approximate the behavior of stiffer polymers, such as
polyethylene. The phase contrast is 600 with respect to the shear moduli and
the Poisson ratios are ${\nu_{\rm m}=0.4995}$ and ${\nu_{\rm f}=0.3636}$.

The training boundary conditions are defined with $N_{\rm dir}=128$ deviatoric
directions, $\underline{N}^{(n)}$, and $N_{\rm amp}=10$ deviatoric amplitudes,
$t^{(p)}$, which are regularly distributed in the interval $[0.05,0.5]$. In
order to also consider response to volumetric extension in the training data,
an additional set of 10 boundary conditions of the form
$\left(\overline{J}^{(m)}\right)^{1/3}\\!\mbox{\boldmath$I$}$ is included in
the training set, with the determinant $\overline{J}^{(m)}$ being linearly
increased from 1 to 1.02.

The validation conditions are chosen along $N_{\rm dir}=64$ new deviatoric
directions at the same deviatoric amplitudes as the training data. In order to
render these validation conditions more different from the training set, 0.5%
compression is proportionally applied.

The results for various values $N$ of the RB-size are compared with the
results of FE simulations with the same boundary conditions. To this end, the
error measures

$\displaystyle{\rm err}_{\rm W}$ $\displaystyle=\frac{\lVert\overline{W}^{\rm
RB}-\overline{W}^{\rm FEM}\rVert}{\lVert\overline{W}^{\rm
FEM}\rVert}\text{\hbox to0.0pt{$\quad\qquad$and\hss}}$ $\displaystyle{\rm
err}_{\rm P}$ $\displaystyle=\frac{\lVert\overline{\mbox{\boldmath$P$}}^{\rm
RB}-\overline{\mbox{\boldmath$P$}}^{\rm
FEM}\rVert}{\lVert\overline{\mbox{\boldmath$P$}}^{\rm FEM}\rVert}$ (46)

are employed. Figures 6 and 7 visualize the results.

Figure 6: Cumulative energy error distribution per direction for the _RB of
the fibrous microstructure_ under validation boundary conditions. Figure 7:
Cumulative stress error distribution per direction for the _RB of the fibrous
microstructure_ under validation boundary conditions.

The distribution of the energy error, ${\rm err}_{\rm W}$, improves
monotonically as the RB is enriched from $N=8$ to $N=128$ elements. This
enrichment corresponds to the inclusion of additional subtrahends in the
computation of $\overline{{\mathbb{C}}}^{\rm RB}$, improving the spectral
over-estimation by the volume average of the stiffness, cf. (17). It is also
noteworthy that the error tends to be higher for larger magnitudes of the
applied kinematic boundary condition, although that is not always the case.

In contrast to this, the stress error ${\rm err}_{\rm P}$ distribution
monotonically improves only up to a certain threshold value of the number of
basis elements, which is $N=64$ in this example. For the greater bases with
$N=96$ and $N=128$, the quality of the results deteriorates as far as the
stress error is concerned. This is most likely due to an excessively
oscillatory nature of the higher order modes: at some critical level $1\ll
i=N_{\rm crit}<N$, the POD constructs eigenvectors $\underline{E}^{(i)}$ with
the $L^{2}(\varOmega_{0})$-norm $\sqrt{\lambda_{i}}\ll 1$. Therefore, the POD
would construct basis vectors out of _numerical fluctuations_ , which would be
unphysical. While the enrichment of the optimization space with unphysical
information cannot increase the minimum energy error $\rm err_{\rm W}$, it
might lead to fluctuations in the displacement field that have significant
impact on the overall stress response. This is especially the case for
numerical fluctuations within the stiff inclusion phase where low overall
displacement errors still could lead to notable impact on the effective
stress.

Nonetheless, all observed errors are less than 20% and stay below 3% for the
optimal sampled size $N=64$. For half the basis size, $N=32$, the errors max
at approximately 5%, which is still acceptable considering the uncertainties
involved in realistic two-scale simulations. Note that the maximum errors
strongly depend on the maximum load amplitude, which is chosen very large in
this example (50% deviatoric strain and 0.5% compression).

The runtimes of the RB model for different sizes $N$ are depicted in Figure 8.
A nearly linear growth of the runtime with respect to the basis size can be
asserted. It is noteworthy that the online time of the RB method is strongly
dominated by the assembly of the Jacobian $\underline{\underline{D}}$.
Therefore, a _Quasi-Newton implementation_ was chosen, resulting in only two
assemblies per load increment.

Figure 8: Laptop computer runtimes of the RB model for the fibrous
microstructure for various sizes $N$ of the basis. Each data point represents
the time needed for all 10 load increments. The spread of the individual times
of the 64 validation cases around these average values is negligible.

Speed-ups become impressive when very large load increments are considered. In
all examples observed thus far, _the RB converges to the final load amplitude
in a single increment, requiring 7-13 Quasi-Newton iterations with only 2-4
assemblies of the Jacobi matrix $\underline{\underline{D}}$ and a runtime of
10-50 seconds._ This is in strong contrast to the FEM that is very sensitive
to large load increments as they come along with a high probability of a
violation of the condition ${\rm det}(\mbox{\boldmath$F$})>0$. By means of
standard implementation, such occurrences are usually treated with cutbacks of
the load increment, being detrimental to the runtime of the FEM.

No rigorous speed-up analysis is intended at this point. Both the codes of the
FEM and of the RB method are fairly efficient in-house C/C++ developments and
perform exact line searches. While the FEM has not yet been equipped with a
Quasi-Newton procedure, the linear solver makes use of parallelization. This
is contrary to the current implementation of the RB method. Depending on the
microstructure (especially the geometry, material nonlinearities, and phase
contrast), the loading conditions, and the size $N$ of the RB, _speed-up
factors are in the order of 5-100_.

###  4.2 Reduced Basis for a stiffening microstructure

The second example takes the “worst-case” approach further by considering a
non-cubical microstructure with even higher phase contrast and significant
topological nonlinearity. To this end, a cuboid microstructure occupying the
domain $[-0.5,0.5]\times[-0.3,0.3]\times[-0.05,0.05]\subset{\mathbb{R}}^{3}$
and containing a hash-like inclusion is investigated, see Figure 9(a). The
mesh is periodic and contains 33923 nodes in 21726 quadratic TET10 elements
(C3D10). The reinforcement makes up approximately 13.3% of the volume. Due to
this large volume fraction, a pronounced geometry-induced nonlinearity of the
effective response is expected under uniaxial loading conditions along the
$x$-axis. As it is elongated, the hash-like part is straightened and thus
increasingly aligned with the external loading, see Figure 9(b). Such effects
might be desirable when designing microstructures for functional materials.

*(a)[.33][c]*(b)[.33][c] *(c)[.33][c]

Figure 9: _(a)_ Cuboidal microstructure with hash-like inclusion phase. _(b)_
Deformed state under uniaxial tension loading. Only inclusion is shown,
coloring indicates $\overline{\mbox{\boldmath$P$}}_{xx}$. _(c)_ Straight
inclusion substitute microstructure, leading to a comparable effective stress.

The material parameters are ${K_{\rm m}=19.867}$ MPa, ${G_{\rm m}=0.4}$ MPa,
${K_{\rm f}=19867}$ MPa, and ${G_{\rm f}=240}$ MPa, implying a Poisson ratio
of $0.49$ in both materials and a phase contrast of $1000$.

The training boundary conditions are the deviatoric ones of the set considered
in Section 4.1, i.e. $N_{\rm dir}=128$ deviatoric directions and $N_{\rm
amp}=10$ regularly spaced deviatoric amplitudes from the intervall
$[0.05,0.5]$. No dilatational training cases are considered, i.e. only points
from a five-dimensional submanifold of the space $\overline{{\cal U}}$ are
sampled.

Uniaxial tension boundary conditions are applied for the validation. More
precisely, the stretch component $\overline{\mbox{\boldmath$U$}}_{\rm xx}$ is
increased from $1.0$ to $1.5$ in $10$ increments of equal size. The other
components are chosen such that all but the $\rm xx$-component of the
effective stress $\overline{\mbox{\boldmath$P$}}$ vanish.

Figure 10 depicts the results for different sizes $N$ of the RB. The influence
of the stiffening effect on the stress curve is emphasized by the black dashed
line corresponding to a similar microstructure with a straight, cuboid
inclusion that leads to the same final stress value under these boundary
conditions, see Figure 9(c).

Figure 10: Stress curves for a microstructure with geometric stiffening (cf.
Figure 9 (a)), comparing the FEM result to the RB for various number of basis
elements $N$. A similar microstructure without geometric stiffening but with
the same final stress value (cf. Figure 9 (c)) leads to the black dashed
curve.

In this example, the geometric stiffening effect is captured by the RB with
high accuracy, with as few as $N=24$ basis elements. For moderate stretches,
even an RB size of $N=16$ suffices. These results are somewhat more impressive
when noticing that the applied boundary condition contains more than 1.2%
volumetric compression, i.e. the validation loading actually lies outside the
submanifold covered during training.

In order to examine the action of the cutoff function $\phi$, the following
two indicators are introduced:

$\displaystyle c_{\rm qp}$ $\displaystyle=\text{\\# quadrature points with
}\left(\phi(J)<1\right)$ (47) $\displaystyle V_{\rm excl}$
$\displaystyle=\left(\lvert\varOmega_{0}\rvert-\sum_{p=1}^{N_{\rm
qp}}\phi(J_{p})w_{p}\right)/\lvert\varOmega_{0}\rvert$ (48)

The first quantity, $c_{\rm qp}$, counts the number of quadrature points at
which the cutoff function has an influence. The second one, $V_{\rm excl}$,
measures the relative excluded volume, interpreting the value of $\phi$ as a
scaling of the corresponding quadrature weight. The values of these indicators
are depicted in Figure 11.

Figure 11: _(Left)_ Number of quadrature points at which the cutoff function
$\phi$ attains a value less than one. _(Right)_ Relative excluded volume.

Most notably, the cutoff function does not have any effect before the eighth
load increment in this example. Only for large load amplitudes does this
numerical stability tweak become necessary. Even then, the number of points at
which it has an influence is small, considering the total number of quadrature
points, $N_{\rm qp}=86904$. This example is representative for all conducted
numerical experiments.

##  5 Discussion

###  5.1 Discussion of the Reduced Basis method

####  5.1.1 Relation of the RB homogenization to analytical estimates

Zero coefficients, $\underline{\xi}=\underline{0}$, correspond to the Taylor
homogenization, i.e. to the nonlinear counterpart of the Voigt estimate
(Voigt, 1966), which provides an upper bound for the material response, cf.
(17). Starting with the initial guess $\underline{\xi}^{(0)}=\underline{0}$,
the evolution of the coefficients corresponds to a relaxation of this overly
stiff response into a more natural state. In view of improved computational
efficiency, a non-zero initial guess $\underline{\xi}^{(0)}$ combined with an
exact line search has proven reasonable and easy to realize. For instance,
such a guess might stem from previous load steps or an
interpolation/extrapolation of available coefficient data.

####  5.1.2 Reconstruction of displacement fields

Given an RB approximation of the deformation gradient,
$\mbox{\boldmath$F$}^{\rm RB}$, one can reconstruct the corresponding
displacement field uniquely up to rigid body motion. This is possible due to
the linear dependence of the deformation gradient fluctuations on the
displacement fluctuations. Recall the definition of the RB in (32),

$\displaystyle\mbox{\boldmath$B$}^{(i)}$ $\displaystyle=\sum_{j=1}^{N_{\rm
s}}\frac{1}{\sqrt{\lambda_{i}}}\left(\underline{E}^{(i)}\right)_{j}\widetilde{\mbox{\boldmath$F$}}^{(j)}\quad(i=1,\ldots,N).$
The corresponding displacement fluctuations are
$\displaystyle\widetilde{\mbox{\boldmath$u$}}_{\rm B}^{(i)}$
$\displaystyle=\sum_{j=1}^{N_{s}}\frac{1}{\sqrt{\lambda_{i}}}\left(\underline{E}^{(i)}\right)_{j}\widetilde{\mbox{\boldmath$u$}}^{(j)}\quad(i=1,\ldots,N).$
(49) The displacement fluctuation fields
$\widetilde{\mbox{\boldmath$u$}}^{(j)}$ are defined by
$\widetilde{\mbox{\boldmath$u$}}^{(j)}(\mbox{\boldmath$X$})=\mbox{\boldmath$u$}^{(j)}(\mbox{\boldmath$X$})-{\overline{\mbox{\boldmath$H$}}}^{(j)}\mbox{\boldmath$X$}$,
where the displacement fields $\mbox{\boldmath$u$}^{(j)}$ are the solutions
computed during training, and
${\overline{\mbox{\boldmath$H$}}}^{(j)}={\overline{\mbox{\boldmath$U$}}}^{(j)}-\mbox{\boldmath$I$}$.
Thus, a displacement field compatible to the RB result
$\mbox{\boldmath$F$}^{\rm
RB}(\mbox{\boldmath$X$};\overline{\mbox{\boldmath$F$}})$ is given by
$\displaystyle\mbox{\boldmath$u$}^{\rm
RB}(\mbox{\boldmath$X$};\overline{\mbox{\boldmath$F$}})$
$\displaystyle=\overline{\mbox{\boldmath$H$}}\mbox{\boldmath$X$}+\sum_{i=1}^{N}\xi^{*}_{i}(\overline{\mbox{\boldmath$F$}})\widetilde{\mbox{\boldmath$u$}}_{\rm
B}^{(i)}(\mbox{\boldmath$X$}).$ (50)

The missing term
$\overline{\mbox{\boldmath$u$}}(\overline{\mbox{\boldmath$X$}})$, cf. (10),
cannot be reconstructed.

####  5.1.3 Relation to classical displacement-based POD methods

In a certain sense, the entries of the correlation matrix used in the offline
phase, cf. Section 2.2, are _“weighted” scalar products_ of the displacement
fluctuation fields $\widetilde{\mbox{\boldmath$w$}}^{(i)}$ within the Sobolev
space $H^{1}_{0}(\varOmega_{0})$. “Weighted” is to be understood in that the
zeroth order derivative is multiplied by zero. Classical displacement-based
POD methods compute correlations of the fluctuations
$\widetilde{\mbox{\boldmath$w$}}^{(i)}$ within the Lebesque space
$L^{2}(\varOmega_{0})$. The change to $H^{1}_{0}(\varOmega_{0})$-like scalar
products is physically motivated by the fact that the local energy
$W=W(\mbox{\boldmath$F$})$ explicitly depends on the gradient of the
displacement, $\mbox{\boldmath$F$}=\mbox{\boldmath$u$}\otimes\nabla_{\rm
X}+\mbox{\boldmath$I$}$, but _not_ on the displacement, $u$.

####  5.1.4 Advantages compared to general displacement-based schemes

The proposed method is advantageous compared to both displacement-based POD
methods and the classical FEM for the following reasons:

  * •

_No gradients need to be computed from displacement fields_ , which
displacement-based schemes always require prior to the evaluation of the
material law.

  * •

The residual $\underline{r}$ and the Jacobian $\underline{\underline{D}}$ are
_algorithmically sleek and trivial to implement_.

  * •

The _absence of element formulations_ in the assembly of the residual
$\underline{r}$ and of the Jacobi matrix $\underline{\underline{D}}$
contributes not only to the simplicity of the method but is also in favor of
massively parallel implementation. This implementation is outstanding for the
current method, but has been conducted for related problems in the small
strain setting in (Fritzen and Hodapp, 2016).

Although the storage of the basis $\\{\mbox{\boldmath$B$}^{(i)}\\}_{i=1}^{N}$
requires $9N_{\rm qp}N$ double precision values, the basis is compact enough
to be completely loaded into random access memory. For instance, the bases
considered in Section 4 occupy only ~200 Mb of memory for $N=32$. Moreover,
the data is aligned efficiently allowing for _linear memory access_.

We now address the _algorithmic complexity_ associated with the proposed
$F$-RB method and with the $u$-RB method that was employed in previous works,
such as (Radermacher and Reese, 2016) and (Yvonnet and He, 2007). To this end,
the fully discretized versions of the residual $\underline{r}$ and of the
Jacobian $\underline{\underline{D}}$ as well as discrete quantities associated
with the $u$-RB method are introduced in the following listing.

  * •

$\underline{P}(\mbox{\boldmath$X$}_{p})\in{\mathbb{R}}^{9}$: nine values of
the stress
$\mbox{\boldmath$P$}(\mbox{\boldmath$F$}_{\xi}(\mbox{\boldmath$X$}_{p}))$ at
the quadrature point $\mbox{\boldmath$X$}_{p}\in\varOmega_{0}$

  * •

$\underline{\underline{C}}(\mbox{\boldmath$X$}_{p})\in{\mathbb{R}}^{9\times
9}$: symmetric stiffness tensor

  * •

$\underline{\underline{B}}(\mbox{\boldmath$X$}_{p})\in{\mathbb{R}}^{9\times
N}$: $F$-RB matrix containing the nine values of each basis elements
$\mbox{\boldmath$B$}^{(i)}$ as columns

  * •

$w_{p}$: the quadrature weight at $\mbox{\boldmath$X$}_{p}$

  * •

$N_{\rm DOF}$: three times the number of nodes, $\sim N_{\rm qp}$

  * •

$\underline{\underline{B}}^{\rm u}\in{\mathbb{R}}^{N_{\rm DOF}\times N}$:
$u$-RB matrix of which the columns contain the nodal displacement values

  * •

$\underline{r}^{\,\rm FE}\in{\mathbb{R}}^{N_{\rm DOF}}$: FE global residual
vector

  * •

$\underline{\underline{K}}^{\,\rm FE}\in{\mathbb{R}}^{N_{\rm DOF}\times N_{\rm
DOF}}$: global FE stiffness matrix

RB method | quantity | complexity
---|---|---
$F$-based |  $\displaystyle\underline{r}=\sum_{p=1}^{N_{\rm qp}}\underline{\underline{B}}(\mbox{\boldmath$X$}_{p})^{\sf T}\underline{P}(\mbox{\boldmath$X$}_{p})w_{p}$ | ${\cal O}(9NN_{\rm qp})$
$\displaystyle\underline{\underline{D}}=\sum_{p=1}^{N_{\rm qp}}\underline{\underline{B}}(\mbox{\boldmath$X$}_{p})^{\sf T}\underline{\underline{C}}(\mbox{\boldmath$X$}_{p})\underline{\underline{B}}(\mbox{\boldmath$X$}_{p})w_{p}$ | ${\cal O}([9^{2}N+9N^{2}]N_{\rm qp})$
$u$-based | $\displaystyle\underline{r}=\left(\underline{\underline{B}}^{\rm u}\right)^{\sf T}\underline{r}^{\rm FE}$ |  ${\cal O}(NN_{\rm DOF})$ \+ assembly of $\underline{r}^{\rm FE}$
$\displaystyle\underline{\underline{D}}=\left(\underline{\underline{B}}^{\rm u}\right)^{\sf T}\underline{\underline{K}}^{\rm FE}\underline{\underline{B}}^{\rm u}$ |  ${\cal O}(N^{2}N_{\rm DOF})$ \+ assembly of $\underline{\underline{K}}^{\rm FE}$
Table 1: Algorithmic complexities of the well-established $u$-based RB method
and the novel $F$-based RB method. The assembly of the FE residual and of the
FE stiffness depend on $N_{\rm qp}$.

Table 1 compares the _algorithmic complexity_ of the presented $F$-RB method
with that of standard $u$-RB schemes. First of all, both methods share a
_quadratic_ dependence of their Jacobi matrices on the number of modes, $N$.
Therefore, the assembly of the Jacobian is usually the most costly operation.
Secondly, both approaches’ complexities suffer a linear dependency on the
number of quadrature points, $N_{\rm qp}$. In the displacement-based approach,
this is included in the assembly of the residual and of the stiffness, which
relate to the factor 9 and 92, respectively. Thirdly, the novel $F$-RB scheme
spares the computational overhead associated with FE formulations
$\underline{r}^{\rm FE}$ and $\underline{\underline{K}}^{\rm FE}$.

####  5.1.5 Outlook

Future research should aim at an application of the introduced Reduced Basis
method within realistic two-scale simulations, in analog to (Fritzen and Kunc,
2018; Fritzen and Hodapp, 2016; Rambausek et al., 2019; Kochmann et al.,
2017). Further, modifications of the cutoff function, $\phi$, should be
investigated. A function with compact support might be more appropriate and
possibly even beneficial to the accuracy of the method. The long-term
perspective is to extend this framework to the context of dissipative
materials.

###  5.2 Discussion of the sampling strategy

The proposed sampling strategy is meant to serve as a template. As exemplified
in Section 4.1, the samplings can be modified and still lead to a coverage of
the set of macroscopic boundary conditions that is sufficient for the problem
at hand. The example of Section 4.2 took this idea further and showed that it
might not even be necessary to sample the macroscopic determinant. Hence, the
sampling can sometimes be reduced to the five-dimensional subspace of
isochoric macroscopic stretch tensors.

In any case, the exact choice of both the inputs to Algorithm 1 and the
distributions of the deviatoric amplitudes and the macroscopic determinants
remains based on knowledge and sophisticated guesses, at least at the current
state of the art. Further research on this matter might possibly lead to a
refined alternative to Algorithm 1, possibly involving the evaluation of error
estimators.

Conceptualization, Oliver Kunc and Felix Fritzen; Funding acquisition, Felix
Fritzen; Methodology, Oliver Kunc and Felix Fritzen; Project administration,
Felix Fritzen; Software, Oliver Kunc; Supervision, Felix Fritzen; Writing –
original draft, Oliver Kunc and Felix Fritzen.

This research was funded by Deutsche Forschungsgemeinschaft (DFG) within the
Emmy-Noether programm under grant DFG-FR2702/6.

The authors declare no conflict of interest.

The following abbreviations are used in this manuscript:

POD | Proper Orthogonal Decomposition
---|---
DOF | degree(s) of freedom
RB | Reduced Basis
FE(M) | Finite Element (Method)
s.p.d. | symmetric positive definite
DDMS | Dilatational-Deviatoric Multiplicative Split
FOM | full-order model

no

##  Appendix A Material objectivity

The components of the stiffness tensor ${\mathbb{C}}$ show the following
transformation behavior.

$\displaystyle
C_{ijkl}(\mbox{\boldmath$F$})=\displaystyle\frac{\displaystyle\partial
P_{ij}(\mbox{\boldmath$F$})}{\displaystyle\partial F_{kl}}$
$\displaystyle\overset{\eqref{eq:sampling:objectivity:P}}{=}\sum_{m=1}^{3}\displaystyle\frac{\displaystyle\partial
R_{im}P_{mj}(\mbox{\boldmath$U$})}{\displaystyle\partial
F_{kl}}=\sum_{m=1}^{3}\displaystyle\frac{\displaystyle\partial
R_{im}}{\displaystyle\partial
F_{kl}}P_{mj}(\mbox{\boldmath$U$})+R_{im}\displaystyle\frac{\displaystyle\partial
P_{mj}(\mbox{\boldmath$U$})}{\displaystyle\partial F_{kl}}$ (51)
$\displaystyle\overset{\phantom{\eqref{eq:sampling:objectivity:P}}}{=}\sum_{m,n,o=1}^{3}\delta_{ik}U^{-1}_{lm}P_{mj}(\mbox{\boldmath$U$})+R_{im}\displaystyle\frac{\displaystyle\partial
P_{mj}(\mbox{\boldmath$U$})}{\displaystyle\partial
U_{no}}\displaystyle\frac{\displaystyle\partial U_{no}}{\displaystyle\partial
F_{kl}}$ (52)
$\displaystyle\overset{\phantom{\eqref{eq:sampling:objectivity:P}}}{=}\sum_{m,n,o=1}^{3}\delta_{ik}U^{-1}_{lm}P_{mj}(\mbox{\boldmath$U$})+R_{im}C_{mjno}(\mbox{\boldmath$U$})R_{kn}\delta_{ol}$
(53)
$\displaystyle\overset{\phantom{\eqref{eq:sampling:objectivity:P}}}{=}\sum_{m,n=1}^{3}\delta_{ik}U^{-1}_{lm}P_{mj}(\mbox{\boldmath$U$})+R_{im}C_{mjnl}(\mbox{\boldmath$U$})R_{kn}$
(54)

Here, $\delta_{ij}$ is the Kronecker symbol, and $i,j,k,l=1,2,3$.

##  Appendix B Effective material responses of the RB

Let ${\mathbb{I}}$ denote the fourth order identity tensor and let the
arguments of the $F$-RB approximation (23) be omitted, i.e. here
$\mbox{\boldmath$F$}^{\rm RB}=\mbox{\boldmath$F$}^{\rm
RB}(\mbox{\boldmath$X$};\overline{\mbox{\boldmath$F$}})$. Its derivative with
respect to the boundary condition $\overline{\mbox{\boldmath$F$}}$ is

$\displaystyle\displaystyle\frac{\displaystyle\partial\mbox{\boldmath$F$}^{\rm
RB}}{\displaystyle\partial\overline{\mbox{\boldmath$F$}}}$
$\displaystyle={\mathbb{I}}+\sum_{i=1}^{N}\mbox{\boldmath$B$}^{(i)}\otimes\displaystyle\frac{\displaystyle\partial\xi_{i}^{*}}{\displaystyle\partial\overline{\mbox{\boldmath$F$}}}(\overline{\mbox{\boldmath$F$}}).$
(55)

###  B.1 Effective stress

$\displaystyle\begin{split}\overline{\mbox{\boldmath$P$}}^{\rm
RB}&=\displaystyle\frac{\displaystyle\partial\overline{W}^{\rm
RB}}{\displaystyle\partial\overline{\mbox{\boldmath$F$}}}=\displaystyle\frac{\displaystyle\partial}{\displaystyle\partial\overline{\mbox{\boldmath$F$}}}\left<W^{\rm
RB}\right>=\left<\displaystyle\frac{\displaystyle\partial W^{\rm
RB}}{\displaystyle\partial\mbox{\boldmath$F$}}\cdot\displaystyle\frac{\displaystyle\partial\mbox{\boldmath$F$}^{\rm
RB}}{\displaystyle\partial\overline{\mbox{\boldmath$F$}}}\right>\\\
&\overset{\eqref{eq:app:effective:aux}}{=}\left<\mbox{\boldmath$P$}^{\rm
RB}\right>+\sum_{i=1}^{N}\left<\mbox{\boldmath$P$}^{\rm
RB}\cdot\left(\mbox{\boldmath$B$}^{(i)}\otimes\displaystyle\frac{\displaystyle\partial\xi_{i}^{*}}{\displaystyle\partial\overline{\mbox{\boldmath$F$}}}(\overline{\mbox{\boldmath$F$}})\right)\right>\\\
&\overset{\text{Note
\eqref{note:orthogonality}}}{=}\left<\mbox{\boldmath$P$}^{\rm
RB}\right>\end{split}$ (56)

###  B.2 Effective stiffness

$\displaystyle\begin{split}\overline{{\mathbb{C}}}^{\rm
RB}&=\displaystyle\frac{\displaystyle\partial\overline{\mbox{\boldmath$P$}}^{\rm
RB}}{\displaystyle\partial\overline{\mbox{\boldmath$F$}}}\overset{\eqref{eq:app:effective:stress}}{=}\displaystyle\frac{\displaystyle\partial}{\displaystyle\partial\overline{\mbox{\boldmath$F$}}}\left<\mbox{\boldmath$P$}^{\rm
RB}\right>=\left<\displaystyle\frac{\displaystyle\partial{}^{2}W^{\rm
RB}}{\displaystyle\partial\overline{\mbox{\boldmath$F$}}\partial\mbox{\boldmath$F$}}\right>=\left<\displaystyle\frac{\displaystyle\partial^{2}W^{\rm
RB}}{\displaystyle\partial{\mbox{\boldmath$F$}}^{2}}\cdot\displaystyle\frac{\displaystyle\partial\mbox{\boldmath$F$}^{\rm
RB}}{\displaystyle\partial\overline{\mbox{\boldmath$F$}}}\right>\\\
&\overset{\eqref{eq:app:effective:aux}}{=}\left<{\mathbb{C}}^{\rm
RB}\right>+\sum_{i=1}^{N}\left<{\mathbb{C}}^{\rm
RB}\cdot\left(\mbox{\boldmath$B$}^{(i)}\otimes\displaystyle\frac{\displaystyle\partial\xi_{i}^{*}}{\displaystyle\partial\overline{\mbox{\boldmath$F$}}}(\overline{\mbox{\boldmath$F$}})\right)\right>\\\
&=\left<{\mathbb{C}}^{\rm RB}\right>+\sum_{i=1}^{N}\left<{\mathbb{C}}^{\rm
RB}\cdot\mbox{\boldmath$B$}^{(i)}\right>\otimes\displaystyle\frac{\displaystyle\partial\xi_{i}^{*}}{\displaystyle\partial\overline{\mbox{\boldmath$F$}}}(\overline{\mbox{\boldmath$F$}})\end{split}$
(57)

For
$\displaystyle\frac{\displaystyle\partial\xi_{i}^{*}}{\displaystyle\partial\overline{\mbox{\boldmath$F$}}}(\overline{\mbox{\boldmath$F$}})$,
we demand that the residual
$r_{i}(\overline{\mbox{\boldmath$F$}},\underline{\xi})$ from (28) is stable
with respect to the boundary condition $\overline{\mbox{\boldmath$F$}}$ when
converged to
$r_{i}(\overline{\mbox{\boldmath$F$}},\underline{\xi}^{*}(\overline{\mbox{\boldmath$F$}}))=0$,

$\displaystyle\begin{split}\displaystyle\frac{\displaystyle\partial
r_{i}}{\displaystyle\partial\overline{\mbox{\boldmath$F$}}}(\overline{\mbox{\boldmath$F$}},\underline{\xi}^{*}(\overline{\mbox{\boldmath$F$}}))&=\left<\mbox{\boldmath$B$}^{(i)}\cdot\displaystyle\frac{\displaystyle\partial\mbox{\boldmath$P$}^{\rm
RB}}{\displaystyle\partial\overline{\mbox{\boldmath$F$}}}\right>\\\
&=\left<\mbox{\boldmath$B$}^{(i)}\cdot\left(\displaystyle\frac{\displaystyle\partial\mbox{\boldmath$P$}^{\rm
RB}}{\displaystyle\partial\mbox{\boldmath$F$}}\displaystyle\frac{\displaystyle\partial\mbox{\boldmath$F$}^{\rm
RB}}{\displaystyle\partial\overline{\mbox{\boldmath$F$}}}\right)\right>\\\
&=\left<\mbox{\boldmath$B$}^{(i)}\cdot{\mathbb{C}}^{\rm
RB}\right>+\sum_{j=1}^{N}\left<\mbox{\boldmath$B$}^{(i)}\cdot{\mathbb{C}}^{\rm
RB}\cdot\displaystyle\frac{\displaystyle\partial\mbox{\boldmath$F$}^{\rm
RB}}{\displaystyle\partial\xi_{j}^{*}}\displaystyle\frac{\displaystyle\partial\xi_{j}^{*}}{\displaystyle\partial\overline{\mbox{\boldmath$F$}}}\right>\\\
&=\left<\mbox{\boldmath$B$}^{(i)}\cdot{\mathbb{C}}^{\rm
RB}\right>+\sum_{j=1}^{N}\underbrace{\left<\mbox{\boldmath$B$}^{(i)}\cdot{\mathbb{C}}^{\rm
RB}\cdot\mbox{\boldmath$B$}^{(j)}\right>}_{D_{ij}}\displaystyle\frac{\displaystyle\partial\xi_{j}^{*}}{\displaystyle\partial\overline{\mbox{\boldmath$F$}}}\quad\overset{!}{=}\mbox{\boldmath$0$}.\end{split}$
(58)

It follows that

$\displaystyle\begin{split}\displaystyle\frac{\displaystyle\partial\xi_{j}^{*}}{\displaystyle\partial\overline{\mbox{\boldmath$F$}}}(\overline{\mbox{\boldmath$F$}})&=-\sum_{i=1}^{N}\left(\underline{\underline{D}}^{-1}\right)_{ij}\left<\mbox{\boldmath$B$}^{(i)}\cdot{\mathbb{C}}^{\rm
RB}\right>.\end{split}$ (59)

##  Appendix C Basis for symmetric traceless second order tensors

$\displaystyle\begin{split}\underline{\underline{X}}^{(1)}&=\sqrt{\frac{1}{6}}\begin{bmatrix}2&\phantom{-}0&\phantom{-}0\\\
0&-1&\phantom{-}0\\\
0&\phantom{-}0&-1\end{bmatrix}\qquad\underline{\underline{X}}^{(2)}=\sqrt{\frac{1}{2}}\begin{bmatrix}0&\phantom{-}0&\phantom{-}0\\\
0&\phantom{-}1&\phantom{-}0\\\ 0&\phantom{-}0&-1\end{bmatrix}\\\
\underline{\underline{X}}^{(3)}&=\sqrt{\frac{1}{2}}\begin{bmatrix}0&\phantom{-}1&\phantom{-}0\\\
1&\phantom{-}0&\phantom{-}0\\\
0&\phantom{-}0&\phantom{-}0\end{bmatrix}\qquad\underline{\underline{X}}^{(4)}=\sqrt{\frac{1}{2}}\begin{bmatrix}0&\phantom{-}0&\phantom{-}1\\\
0&\phantom{-}0&\phantom{-}0\\\
1&\phantom{-}0&\phantom{-}0\end{bmatrix}\qquad\underline{\underline{X}}^{(5)}=\sqrt{\frac{1}{2}}\begin{bmatrix}0&\phantom{-}0&\phantom{-}0\\\
0&\phantom{-}0&\phantom{-}1\\\
0&\phantom{-}1&\phantom{-}0\end{bmatrix}\end{split}$ (60)

References

## References

  * Rendek and Lion (2010) Rendek, M.; Lion, A.  Amplitude dependence of filler-reinforced rubber: Experiments, constitutive modelling and FEM – Implementation.  International Journal of Solids and Structures 2010, 47, 2918 – 2936.  doi:black10.1016/j.ijsolstr.2010.06.021.
  * Nguyen et al. (2016) Nguyen, V.D.; Lani, F.; Pardoen, T.; Morelle, X.; Noels, L.  A large strain hyperelastic viscoelastic-viscoplastic-damage constitutive model based on a multi-mechanism non-local damage continuum for amorphous glassy polymers.  International Journal of Solids and Structures 2016, 96, 192 – 216.  doi:black10.1016/j.ijsolstr.2016.06.008.
  * Geers and Yvonnet (2016) Geers, M.; Yvonnet, J.  Multiscale modeling of microstructure-property relations.  MRS Bulletin 2016, 41, 610–616.  doi:black10.1557/mrs.2016.165.
  * Saeb et al. (2016) Saeb, S.; Steinmann, P.; Javili, A.  Aspects of Computational Homogenization at Finite Deformations: A Unifying Review From Reuss’ to Voigt’s Bound.  Applied Mechanics Reviews 2016, 68, 050801.  doi:black10.1115/1.4034024.
  * Feyel (1999) Feyel, F.  Multiscale FE2 elastoviscoplastic analysis of composite structures.  Computational Materials Science 1999, 16, 344 – 354\.  doi:black10.1016/S0927-0256(99)00077-4.
  * Yvonnet and He (2007) Yvonnet, J.; He, Q.C.  The reduced model multiscale method (R3M) for the non-linear homogenization of hyperelastic media at finite strains.  Journal of Computational Physics 2007, 223, 341 – 368\.  doi:black10.1016/j.jcp.2006.09.019.
  * Radermacher and Reese (2016) Radermacher, A.; Reese, S.  POD-based model reduction with empirical interpolation applied to nonlinear elasticity.  International Journal for Numerical Methods in Engineering 2016, 107, 477–495.  doi:black10.1002/nme.5177.
  * Radermacher and Reese (2013) Radermacher, A.; Reese, S.  Proper orthogonal decomposition-based model reduction for non-linear biomechanical analysis.  International Journal of Materials Engineering Innovation 2013, 4, 149–165.  doi:black10.1504/IJMATEI.2013.054393.
  * Soldner et al. (2017) Soldner, D.; Brands, B.; Zabihyan, R.; Steinmann, P.; Mergheim, J.  A numerical study of different projection-based model reduction techniques applied to computational homogenisation.  Computational Mechanics 2017, 60, 613–625.  doi:black10.1007/s00466-017-1428-x.
  * An et al. (2009) An, S.; Kim, T.; James, D.L.  Optimizing cubature for efficient integration of subspace deformations.  ACM Trans Graph 2009, 27, 165:1–165:10.  doi:black10.1145/1409060.1409118.
  * Hernández et al. (2017) Hernández, J.; Caicedo, M.; Ferrer, A.  Dimensional hyper-reduction of nonlinear finite element models via empirical cubature.  Computer Methods in Applied Mechanics and Engineering 2017, 313, 687 – 722.  doi:black10.1016/j.cma.2016.10.022.
  * Temizer and Zohdi (2007) Temizer, I.; Zohdi, T.  A numerical method for homogenization in non-linear elasticity.  Computational Mechanics 2007, 40, 281–298.  doi:black10.1007/s00466-006-0097-y.
  * Yvonnet et al. (2013) Yvonnet, J.; Monteiro, E.; He, Q.C.  Computational homogenization method and reduced database model for hyperelastic heterogeneous structures.  Journal for Multiscale Computational Engineering 2013, 11, 201–225.  doi:black10.1615/IntJMultCompEng.2013005374.
  * Fritzen and Kunc (2018) Fritzen, F.; Kunc, O.  Two-stage data-driven homogenization for nonlinear solids using a reduced order model.  European Journal of Mechanics - A/Solids 2018, 69, 201–220.  doi:black10.1016/j.euromechsol.2017.11.007.
  * Kunc (2019) Kunc, O.  GitHub repository ReducedBasisDemonstrator.  https://github.com/EMMA-Group/ReducedBasisDemonstrator.
  * Flory (1961) Flory, P.  Thermodynamic relations for high elastic materials.  Transactions of the Faraday Society 1961, 57, 829–838.  doi:black10.1039/tf9615700829.
  * Bilger et al. (2005) Bilger, N.; Auslender, F.; Bornert, M.; Michel, J.C.; Moulinec, H.; Suquet, P.; Zaoui, A.  Effect of a nonuniform distribution of voids on the plastic response of voided materials: a computational and statistical analysis.  International Journal of Solids and Structures 2005, 42, 517–538.  doi:black10.1016/j.ijsolstr.2004.06.048.
  * Doll and Schweizerhof (1999) Doll, S.; Schweizerhof, K.  On the Development of Volumetric Strain Energy Functions.  Journal of Applied Mechanics 1999, 67, 17–21.  doi:black10.1115/1.321146.
  * Simo (1988) Simo, J.  A framework for finite strain elastoplasticity based on maximum plastic dissipation and the multiplicative decomposition: Part I. Continuum formulation.  Computer Methods in Applied Mechanics and Engineering 1988, 66, 199 – 219.  doi:black10.1016/0045-7825(88)90076-X.
  * Pruchnicki (1998) Pruchnicki, E.  Hyperelastic homogenized law for reinforced elastomer at finite strain with edge effects.  Acta Mechanica 1998, 129, 139–162.  doi:black10.1007/BF01176742.
  * Miehe (2003) Miehe, C.  Computational micro-to-macro transitions for discretized micro-structures of heterogeneous materials at finite strains based on the minimization of averaged incremental energy.  Computer Methods in Applied Mechanics and Engineering 2003, 192, 559 – 591.  doi:black10.1016/S0045-7825(02)00564-9.
  * Castañeda and Suquet (1998) Castañeda, P.P.; Suquet, P.  Nonlinear Composites.  Advances in Applied Mechanics 1998, 34, 172 – 302.  doi:black10.1016/S0065-2156(08)70321-1.
  * Sirovich (1987) Sirovich, L.  Turbulence and the Dynamics of Coherent Structures. Part 1: Coherent structures.  Quarterly of Applied Mathematics 1987, 45, 561 – 571\.  doi:black10.1090/qam/910462.
  * Kabel et al. (2014) Kabel, M.; Böhlke, T.; Schneider, M.  Efficient fixed point and Newton–Krylov solvers for FFT-based homogenization of elasticity at large deformations.  Computational Mechanics 2014, 54, 1497–1514.  doi:black10.1007/s00466-014-1071-8.
  * Quarteroni et al. (2016) Quarteroni, A.; Manzoni, A.; Negri, F.  Reduced Basis Methods for Partial Differential Equations: An Introduction; La matematica per il 3 + 2, vol. 92, Springer, 2016.  doi:black10.1007/978-3-319-15431-2.
  * Ball (1976) Ball, J.M.  Convexity conditions and existence theorems in nonlinear elasticity.  Archive for Rational Mechanics and Analysis 1976, 63, 337–403.  doi:black10.1007/BF00279992.
  * Schneider (2017) Schneider, M.  Beyond polyconvexity: an existence result for a class of quasiconvex hyperelastic materials.  Mathematical Methods in the Applied Sciences 2017, 40, 2084–2089.  doi:black10.1002/mma.4123.
  * Miehe (1996) Miehe, C.  Numerical computation of algorithmic (consistent) tangent moduli in large-strain computational inelasticity.  Computer Methods in Applied Mechanics and Engineering 1996, 134, 223 – 240.  doi:black10.1016/0045-7825(96)01019-5.
  * Faraut (2008) Faraut, J.  Analysis on Lie groups: an introduction; Cambridge University Press, 2008.  doi:black10.1017/cbo9780511755170.
  * Neff et al. (2016) Neff, P.; Eidel, B.; Martin, R.J.  Geometry of Logarithmic Strain Measures in Solid Mechanics.  Archive for Rational Mechanics and Analysis 2016, 222, 507–572.  doi:black10.1007/s00205-016-1007-x.
  * Brauchart and Grabner (2015) Brauchart, J.S.; Grabner, P.J.  Distributing many points on spheres: Minimal energy and designs.  Journal of Complexity 2015, 31, 293 – 326.  doi:black10.1016/j.jco.2015.02.003.
  * Kunc and Fritzen (2018) Kunc, O.; Fritzen, F.  Generation of energy-minimizing point sets on spheres and their application in mesh-free interpolation and differentiation.  Advances in Computational Mathematics 2018.  Under review.
  * Leopardi (2006) Leopardi, P.  A partition of the unit sphere into regions of equal area and small diameter.  Electronic Transactions on Numerical Analysis 2006, 25, 309–327.
  * Kim and Swan (2003) Kim, H.J.; Swan, C.C.  Algorithms for automated meshing and unit cell analysis of periodic composites with hierarchical tri-quadratic tetrahedral elements.  International Journal for Numerical Methods in Engineering 2003, 58, 1683–1711.  doi:black10.1002/nme.828.
  * Voigt (1966) Voigt, W., Lehrbuch der Kristallphysik; Vieweg+Teubner Verlag, 1966.  doi:black10.1007/978-3-663-15884-4.
  * Fritzen and Hodapp (2016) Fritzen, F.; Hodapp, M.  The Finite Element Square Reduced (FE2R) method with GPU acceleration: towards three-dimensional two-scale simulations.  International Journal for Numerical Methods in Engineering 2016, 107, 853–881.  doi:black10.1002/nme.5188.
  * Rambausek et al. (2019) Rambausek, M.; Göküzüm, F.S.; Nguyen, L.T.K.; Keip, M.A.  A two-scale FE-FFT approach to nonlinear magneto-elasticity.  International Journal for Numerical Methods in Engineering 2019, 117, 1117–1142.  doi:black10.1002/nme.5993.
  * Kochmann et al. (2017) Kochmann, J.; Wulfinghoff, S.; Ehle, L.; Mayer, J.; Svendsen, B.; Reese, S.  Efficient and accurate two-scale FE-FFT-based prediction of the effective material behavior of elasto-viscoplastic polycrystals.  Computational Mechanics 2017.  doi:black10.1007/s00466-017-1476-2.