# Comptes Rendus
## Mécanique

Andrea Bonito, Diane Guignard and Angelique Morvant

**Numerical approximations of thin structure deformations**

Volume 351, Special Issue S1 (2023), p. 181-217

Online since: 21 August 2023

**Part of Special Issue:** The scientific legacy of Roland Glowinski

**Guest editors:** Grégoire Allaire (CMAP, Ecole Polytechnique, Institut Polytechnique de Paris, Palaiseau, France),

Jean-Michel Coron (Laboratoire Jacques-Louis Lions, Sorbonne Université)

and Vivette Girault (Laboratoire Jacques-Louis Lions, Sorbonne Université)

https://doi.org/10.5802/crmeca.201

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

This article is licensed under the

CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.

http://creativecommons.org/licenses/by/4.0/
# The scientific legacy of Roland Glowinski / L'héritage scientifique de Roland Glowinski

## Numerical approximations of thin structure deformations

**Andrea Bonito<sup>a</sup>, Diane Guignard<sup>*,b</sup> and Angelique Morvant<sup>a</sup>**

*E-mails*: bonito@math.tamu.edu (A. Bonito), dguignar@uottawa.ca (D. Guignard), mae4102@tamu.edu (A. Morvant)

*In memory of Roland Glowinski*

**Abstract.** We review different models for thin structures using bending as principal mechanism to undergo large deformations. Each model consists of minimizing a fourth order energy, potentially subject to a nonconvex constraint. Equilibrium deformations are approximated using local discontinuous Galerkin finite elements. The discrete energies relies on a discrete Hessian operator defined on discontinuous functions with better approximation properties than the piecewise Hessian. Discrete gradient flows are used to drive the minimization process. They are chosen for their robustness and ability to preserve the nonconvex constraint. Several numerical experiments are presented to showcase the variety of shapes achievable with these models.

**Keywords.** Nonlinear elasticity, plate deformation, folding, prestrain metric, discontinuous Galerkin, reconstructed Hessian, numerical simulations.

**Funding.** Andrea Bonito was partially supported by the NSF Grant DMS-2110811. Diane Guignard was partially supported by NSERC Grant RGPIN-2021-04311. Angelique Morvant was partially supported by the NSF Grants DMS-1817691 and DMS-2110811.

*Manuscript received 24 May 2023, accepted 25 May 2023.*

## 1. Introduction

Deformations of thin materials are widely observed in nature, from the snapping of the venus flytrap [1, 2] to the natural growth of soft tissues in leaves and flowers [3–6]. They also appear in a variety of man-made applications [7, 8]. Some thin structures, called *bilayers*, are made of two thin layers of different materials that react differently to external stimuli (changes in temperature or humidity, electrical or chemical stimuli, etc.). Examples include the bimetal strips in thermostats, microactuators [9, 10], and plywood panels in climate-responsive architectures [11] inspired from the morphology of conifer cones. Bending of thin sheets can also occur when folding is present; examples of this include origami and flexible structures [12, 13]. We refer to the review paper [6] for additional examples.

All the phenomena listed above can be modelled as thin limits of hyperelastic materials. A thin structure of thickness *s* endowed with a hyperelastic energy *E*<sub>s</sub> is approximated by its two

<footnote><sup>a</sup> Department of Mathematics, Texas A&M University, College Station, TX 77845, USA.</footnote>

<footnote><sup>b</sup> Department of Mathematics and Statistics, University of Ottawa, Ottawa, ON K1N 6N5, Canada.</footnote>

<footnote>* Corresponding author.</footnote>
dimensional midplane endowed with a limiting energy as the thickness vanishes. A hierarchy of reduced models for elastic deformations is described in [14], in which the type of the model depends on the scaling of the elastic energy $E_s$ with the thickness $s$ of the plate. The scaling $E_s \sim s$ corresponds to stretching, the scaling $E_s \sim s^3$ corresponds to bending, and the scaling $E_s \sim s^5$ corresponds to the von Kármán theory of plate bending. Each of these reduced energies is the $\Gamma$-limit of $E_s$. We refer to [15] for the membrane theory, to [16–18] for the bending theory, and to [14, 19] for the von Kármán theory. In particular, cluster points of sequences of minimizers of $\{E_s\}_{s>0}$ are minimizers of the $\Gamma$-limit.

In this work, we are mainly interested in the bending regime, namely when the hyperelastic three dimensional energy

$$E_s(\mathbf{u}) = \int_{\Omega \times (-\frac{s}{2}, \frac{s}{2})} W \left( \nabla \mathbf{u}(\mathbf{x}) G^{-\frac{1}{2}}(\mathbf{x}) \right) \mathrm{d}\mathbf{x} \quad (1)$$

scales like the cube of the thickness of the plate $\Omega \times (-\frac{s}{2}, \frac{s}{2})$. Here, $\Omega$ is the midplane of the thin structure, $\mathbf{u}: \Omega \times (-\frac{s}{2}, \frac{s}{2}) \to \mathbb{R}^3$ is a deformation of the plate, $G$ is a prescribed Riemannian metric (deformations satisfying the target metric $G$ are stress-free), and $W$ is some energy density function satisfying [17]

* $W \in C^0(\mathbb{R}^{3 \times 3})$, $W \in C^2$ in a neighborhood of $SO(3)$;

* $W$ is frame indifferent: $W(F) = W(RF)$ for all $F \in \mathbb{R}^{3 \times 3}$ and all $R \in SO(3)$;

* $W(F) \ge \mathrm{Cdist}^2(F, SO(3))$ for all $F \in \mathbb{R}^{3 \times 3}$ and $W(F) = 0$ if $F \in SO(3)$,

where $SO(3)$ denotes the special orthogonal group of rotations in $\mathbb{R}^3$. In this work, we restrict our considerations to the St. Venant–Kirchhoff stored energy, see (8).

The prestrain metric $G$ characterizes the material. From the third property of the energy density, we see that when $G$ is the identity matrix, the minimum of the hyperelastic energy (1) (without boundary conditions or external forces) is zero and is achieved by rigid motions. That is, at equilibrium, the plate is stress free and flat. When the Riemannian curvature tensor of $G$ is not identically zero, then the minimum of the hyperelastic energy is strictly positive [20]. In this case, there is no stress-free configuration. The plate is then said to be non-Euclidean or pretrained and generally exhibits more complex equilibrium shapes.

In 1850, Kirchhoff [21] obtained a reduced energy that can be formally obtained assuming that the material deformation reads

$$\mathbf{u}(x_1, x_2, x_3) = \mathbf{y}(x_1, x_2) + x_3 \mathbf{v}(x_1, x_2), \quad (2)$$

where $\mathbf{v}$ is the normal to the deformed midplane $\mathbf{y}(\Omega)$. This assumption, used by Love [22], is usually referred to as (nonlinear) Kirchhoff–Love ansatz in the literature, although it is not due to Kirchhoff [17]. When the thickness $s$ of the plate is small, it is convenient to derive reduced models for the deformation $\mathbf{y}$ of the midplane $\Omega$. In [23], see also [5], the Kirchhoff–Love assumption (2) is made for pretrained plates. The resulting model consists of the sum of a stretching and a bending component multiplied by different powers of the plate thickness $s$. We derive a similar energy referred to as the preasymptotic energy to express that it is obtained upon assuming that $s$ is small but not vanishing. To do this, we follow [20] and assume that

$$G(x_1, x_2, x_3) = \begin{pmatrix} g(x_1, x_2) & \mathbf{0} \\ \mathbf{0} & 1 \end{pmatrix}, \quad (3)$$

i.e. $G$ is uniform within the thickness and no stretching occurs in the direction orthogonal to the plane. The function $g: \Omega \to \mathbb{R}^{2 \times 2}$ is assumed to be symmetric and uniformly positive definite.

The vanishing thickness limit of pretrained materials has been studied in several works. When $G$ is the identity matrix, as mentioned above, a reduced energy was obtained formally by Kirchhoff [21]. Much later, an ansatz-free derivation was obtained via $\Gamma$-convergence in the
seminal work [17]. The limiting energy involves the second fundamental form of the midplane deformation and is finite only for isometric immersions.

For the derivation in the general case, namely when $G$ is not the identity matrix, we refer to [5, 20, 24–26]. The first $\Gamma$-convergence results for prestrained plates were obtained in [20] for target metrics as in (3), and later extended to more general metrics in [24]. In both cases, the metric is assumed to be independent of the thickness and uniform throughout the thickness. For the particular metric given in (3), the limiting energy involves the second fundamental form of the midplane deformation and is finite only for isometric immersion of the two dimensional target metric $g$. In [27] a formal derivation based on a modified Kirchhoff–Love assumption (2) of the reduced model is provided.

Bilayer materials are composed of two thin layers of materials with different properties. External (thermal, electrical or chemical) stimuli correspond to a prestrain metric

$$G = (I_3 \pm \zeta s N)^T (I_3 \pm s \zeta N), \quad (4)$$

where the product $s\zeta$ describes the lattice mismatch between the layers, and $N$ encodes the inhomogeneity and anisotropy of the bilayer. When actuated, the material deforms to relax its internal stress and can typically undergo large deformations even with relatively small stimuli. The reduced model for bilayer plates was derived via $\Gamma$-convergence in [28] and formally explained in [29] for a metric of the form (4). The limiting energy penalizes deviations of the second fundamental form of the midplane from a spontaneous curvature $Z$ which depends on the mismatch between the two layers (i.e. on $\zeta$ and $N$). In this model, as in the others when $G = I_3$, the energy is finite only when the deformations are isometries.

The ability of the plate to fold along creases was incorporated in [30]. It considers hyperelastic materials with $G = I_3$ assumed to be weakened in a neighborhood of a curve $\Sigma \subset \Omega$ modelling a crease. Assuming an asymptotic behavior of the defect width and strength (see (28)), the limiting model reduces to the original model without the crease, except that jumps on the deformations gradients are allowed across the crease without suffering any energetic penalty. Note that in principle, any plate – with an isometry constraint or prestrain, a single layer or a bilayer – can have folding. However, only the case of a single layer plate with isometry has been derived rigorously so far.

## Motivations and Novelties

One of the goals of this paper is to review different models for the deformation of plates, including single layer and bilayer plates, plates with an isometry constraint as well as general prestrain metric constraints, and plates with folding along some curves inside the domain. Although emphasis is made on the pure bending of plates, the preasymptotic regime, in which both the bending and the stretching/shearing of the plate are considered, is also analyzed. To give an intuitive derivation, the reduced two dimensional energies are obtained formally using a (modified) Kirchhoff–Love assumption for the deformation of the three dimensional plate. However, as mentioned above, ansatz-free limiting energies have been obtained via $\Gamma$-convergence for most cases.

A second goal is to collect in one place the numerical methods proposed recently by the authors and collaborators for the approximation of near minimizers of the discrete energies. These methods are based on a *local discontinuous Galerkin* (LDG) approach for the space discretization and a gradient flow for the minimization of the discrete energy. This methodology has already been successfully applied to a wide range of problems, see for instance [27, 30, 31]. Moreover, the $\Gamma$-convergence as the mesh size goes to zero of the discrete energies obtained using LDG is analyzed in [31–33] and [30] for the preasymptotic, prestrain, bilayer, and folding case,
respectively. Note that other discretizations are also possible. We refer to [34–36] for approaches using the symmetric interior penalty discontinuous Galerkin method, and to [29, 37–39] for a method based on Kirchhoff elements.

This work also offers several novel contributions. We present a recently developed algorithm for the preasymptotic model, extend the existing algorithms to time-dependent input data, present a uniform treatment of folding effects, and introduce an accelerated algorithm that drastically improves the performance of the unconstrained discrete gradient flows. We perform new and challenging simulations highlighting the various types of deformations that can be achieved with the models, as well as the versatility of the code.

## Outline

The rest of this paper is organized as follows. In Section 2, we formally derive the various two dimensional models for the deformation of thin structures and discuss how to incorporate forcing and boundary conditions into these models. In Section 3, we present the spatial discretization of the two dimensional energies based on the LDG method, including the introduction of a discrete Hessian operator – the key component of the method. The discrete gradient flows put forward to minimize the discrete energies are described in Section 4, which also includes strategies to produce good initial deformations and to potentially improve the convergence rate of the iterative processes. Section 5 is devoted to numerical experiments demonstrating the variety of deformations that can be achieved.

## Notation

In what follows, we will denote by $I_n$ the identity matrix in $\mathbb{R}^{n \times n}$. Moreover, for a multivariate function, we write $\partial_i$ the partial derivative with respect to the $i^{th}$ variable. Uppercase letters and bold lowercase letters will be used for matrix-valued and vector-valued functions, respectively, and subindices will denote their components. For instance, the $(i, j)$ component of $M: \mathbb{R}^m \to \mathbb{R}^{n_1 \times n_2}$ will be denoted $M_{ij}: \mathbb{R}^m \to \mathbb{R}$, $1 \le i \le n_1$ and $1 \le j \le n_2$, while the $j^{th}$ component of $v: \mathbb{R}^m \to \mathbb{R}^n$ will be denoted $v_j: \mathbb{R}^m \to \mathbb{R}$, $1 \le j \le n$. Furthermore, for $\Phi: \mathbb{R}^m \to \mathbb{R}^{n \times m}$, $\Phi_j: \mathbb{R}^m \to \mathbb{R}^n$ will denote the $j^{th}$ row of $\Phi$, $1 \le j \le n$. For $\Phi: \mathbb{R}^m \to \mathbb{R}^{m \times m}$, the divergence operator is applied row-wise, namely $[\text{div}(\Phi)]_i = \sum_{j=1}^m \partial_j(\Phi_{ij})$ for $i=1, 2, ..., m$. Finally, we use the notation $A: B$ to denote the Euclidean scalar product between two tensors $A, B$ and $|\cdot|$ for the corresponding Frobenius norm. For the particular case of vectors, the Euclidean scalar product between $v, w \in \mathbb{R}^n$ is instead denoted by $v \cdot w$.

## 2. Mathematical Models

In this section, we provide a formal but intuitive derivation of two dimensional models for three dimensional thin hyperelastic structures. References to rigorous derivations of the bending model with isometry and prestrain constraint and of the bilayer model via $\Gamma$-convergence can be found in Sections 2.2, 2.3, and 2.4, respectively. The thin structure is denoted $\Omega_s := \Omega \times (-\frac{s}{2}, \frac{s}{2}) \subset \mathbb{R}^3$, where $s > 0$ stands for the thickness of the structure and $\Omega \subset \mathbb{R}^2$ is an open, bounded domain with Lipschitz boundary which represents the midplane (see the diagram on the left of Figure 1). A deformation of the plate is denoted by $\mathbf{u}: \Omega_s \subset \mathbb{R}^3 \to \mathbb{R}^3$ while its restriction to the midplane is denoted $\mathbf{y}: \Omega \subset \mathbb{R}^2 \to \mathbb{R}^3$ (see the diagram on the right of Figure 1). For later use, we also introduce the unit normal vector to the surface $\mathbf{y}(\Omega)$ at the point $\mathbf{y}(\mathbf{x}')$

$$ \mathbf{v}(\mathbf{x}') := \frac{\partial_1 \mathbf{y}(\mathbf{x}') \times \partial_2 \mathbf{y}(\mathbf{x}')}{|\partial_1 \mathbf{y}(\mathbf{x}') \times \partial_2 \mathbf{y}(\mathbf{x}')|}, \quad \mathbf{x}' \in \Omega, $$
and the first and second fundamental forms of the surface $\mathbf{y}(\Omega)$

$$I(\mathbf{y}) := (\nabla' \mathbf{y})^T \nabla' \mathbf{y} \quad \text{and} \quad \Pi(\mathbf{y}) := -(\nabla' \mathbf{y})^T \nabla' \mathbf{v} = (\partial_{ij} \mathbf{y} \cdot \mathbf{v})^2_{i,j=1}. \qquad (5)$$

Here $\nabla'$ denotes the gradient with respect to the variable $\mathbf{x}' \in \Omega \subset \mathbb{R}^2$.

<figure><img src="image_2.png" /><figcaption>**Figure 1.** The undeformed reference plate (left) and the deformed midplane (right).</figcaption></figure>

In all the models presented below, the structure is endowed with an energy favoring deformation with Cauchy-Green strain tensor $(\nabla \mathbf{u})^T \nabla \mathbf{u}$ matching a given target metric $G: \Omega_s \to \mathbb{R}^{3 \times 3}$ (symmetric and positive definite), i.e. for which the strain tensor

$$\epsilon(\nabla \mathbf{u}) := \frac{(\nabla \mathbf{u})^T \nabla \mathbf{u} - G}{2} \qquad (6)$$

is as close to zero as possible. Note that $\epsilon(\nabla \mathbf{u}) \equiv 0$ signifies that the first fundamental form of the deformed surface $\mathbf{u}(\Omega_s)$ is $G$ and that $\mathbf{u}$ is an isometric immersion of $G$. Whether a given $n$ dimensional Riemannian manifold has an isometric immersion into $\mathbb{R}^N$ is a long standing problem in differential geometry, see for instance [40, 41]. By the Nash-Kuiper embedding theorem [42–44], there always exists an isometric immersion $\mathbf{u}: \Omega_s \to \mathbb{R}^3$ of $G$ in $W^{1,\infty}(\Omega_s)$. However, the deformation $\mathbf{u}$ cannot be orientation preserving unless the Riemannian curvature tensor of $G$ is identically zero, see [6] and the references therein.

## 2.1. Bending of Isotropic Thin Structures

We start by considering thin structures endowed with the isotropic hyperelastic energy

$$E_s(\mathbf{u}) := \int_{\Omega_s} W(\nabla \mathbf{u}(\mathbf{x})) \mathrm{d}\mathbf{x}, \qquad (7)$$

where $W: \mathbb{R}^{3 \times 3} \to \mathbb{R}$ is the St. Venant-Kirchhoff stored energy density function defined by

$$W(F) := \mu |\epsilon(F)|^2 + \frac{\lambda}{2} \mathrm{tr}(\epsilon(F))^2, \quad \epsilon(F) := \frac{F^T F - I_3}{2}. \qquad (8)$$

Here $F := \nabla \mathbf{u}$ is the deformation gradient assumed to satisfy the orientation condition $\det(\nabla \mathbf{u}) > 0$, $\epsilon(F)$ is the Green-Lagrange strain tensor corresponding to $G = I_3$ in (6), and $\lambda$ and $\mu$ are the Lamé constants. As we shall see below, deformations for which $E_s(\mathbf{u})$ scales like $s$ correspond to a stretching of the midplane (membrane theory), while bending occurs when $E_s(\mathbf{u})$ scales like $s^3$ (bending theory).

Note that in absence of other constraints such as boundary conditions or external forcing, energy (7) is minimized on $SO(3)$, i.e. when the deformations are isometries:

$$(\nabla \mathbf{u})^T \nabla \mathbf{u} = I_3 \quad \text{in } \Omega_s.$$
## 2.2. Limiting Model and Modified Kirchhoff Assumption

Two dimensional models for hyperelastic thin structures are derived in [17]. It is shown that when bending is the chief mechanism of deformations, the energy $E_s$ must scale like $s^3$ and

$$ \Gamma - \lim_{s \to 0^+} s^{-3} E_s = E $$

where for $y: \Omega \to \mathbb{R}^3$

$$ E(y) = \begin{cases} \frac{1}{24} \int_{\Omega} \left( 2\mu |\Pi(y(x'))|^2 + \frac{2\mu\lambda}{2\mu + \lambda} \mathrm{tr}(\Pi(y(x'))^2) \right) \mathrm{d}x', & \text{when } \mathrm{I}(y) = I_2, \\ +\infty, & \text{otherwise.} \end{cases} \quad (9) $$

Part of the $\Gamma$-convergence argument provided in [17] is based on the construction of a recovery sequence given for $x' \in \Omega$ and $x_3 \in (-\frac{s}{2}, \frac{s}{2})$ by

$$ \mathbf{u}(\mathbf{x}', x_3) = \mathbf{y}(\mathbf{x}') + \left( x_3 + \underbrace{\frac{\lambda}{2\mu + \lambda} \mathrm{tr}(\Pi(\mathbf{y}(\mathbf{x}')))}_{=: \beta(\mathbf{x}')} \frac{1}{2} x_3^2 \right) \mathbf{v}(\mathbf{x}') \quad (10) $$

for a midplane deformation $y: \Omega \to \mathbb{R}^3$ with finite energy $E(y) < \infty$. This indicates that the two dimensional model (9) requires that the fibers orthogonal to the midplane are not stretched homogeneously on $\Omega$. This differentiates the theory in [17] from the standard Cosserat (or Kirchhoff-Love) theory based on the ansatz

$$ \mathbf{u}(\mathbf{x}', x_3) = \mathbf{y}(\mathbf{x}') + x_3 \mathbf{v}(\mathbf{x}') $$

and leading to the same two dimensional energy improperly scaled with coefficient $\lambda$ instead of $\frac{2\mu\lambda}{2\mu+\lambda}$ in front of the trace term [16, 17].

### 2.2.1. Preasymptotic

We follow [27] and provide a formal derivation of the preasymptotic energy for thin structures with non-vanishing but small thickness under assumption (10) on the deformations. We are interested in deformations where bending is the chief mechanism and thus also fix $0 < s \ll 1$ and assume that

$$ E_s(\mathbf{u}) \le \Lambda s^3 \quad (11) $$

for a constant $\Lambda$ independent of $s$. We recall that $\nabla'$ stands for the gradient with respect to $x' \in \Omega \subset \mathbb{R}^2$ so that

$$ \nabla \mathbf{u} = \left[ \nabla' \mathbf{y} + x_3 \nabla' \mathbf{v} + \frac{1}{2} x_3^2 \nabla' (\beta \mathbf{v}), \mathbf{v} + x_3 \beta \mathbf{v} \right]. $$

Because $\partial_i y \cdot v = 0$ for $i = 1, 2$ and $v^T v = 1$, we have

$$ \mathbf{v}^T \nabla' \mathbf{y} = \mathbf{0}, \quad (\nabla' \mathbf{y})^T \mathbf{v} = \mathbf{0}, \quad \mathbf{v}^T \nabla' \mathbf{v} = \mathbf{0}, \quad (\nabla' \mathbf{v})^T \mathbf{v} = \mathbf{0}, $$

and thus

$$ \begin{aligned} (\nabla \mathbf{u})^T \nabla \mathbf{u} = & \begin{bmatrix} (\nabla' \mathbf{y})^T \nabla' \mathbf{y} & \mathbf{0} \\ \mathbf{0} & 1 \end{bmatrix} + x_3 \begin{bmatrix} (\nabla' \mathbf{y})^T \nabla' \mathbf{v} + (\nabla' \mathbf{v})^T \nabla' \mathbf{y} & \mathbf{0} \\ \mathbf{0} & 2\beta \end{bmatrix} \\ & + x_3^2 \begin{bmatrix} (\nabla' \mathbf{v})^T \nabla' \mathbf{v} + \frac{1}{2} (\nabla' \mathbf{y})^T \nabla' (\beta \mathbf{v}) + \frac{1}{2} (\nabla' (\beta \mathbf{v}))^T \nabla' \mathbf{y} & \frac{1}{2} (\nabla' (\beta \mathbf{v}))^T \mathbf{v} \\ \frac{1}{2} \mathbf{v}^T \nabla' (\beta \mathbf{v}) & \beta^2 \end{bmatrix} \\ & + \mathcal{O}(x_3^3). \end{aligned} $$
We rewrite the above expression using the first and second fundamental forms (5) of the de-
formed midplane $\mathbf{y}(\Omega)$. In addition, we take advantage of the symmetry of the second funda-
mental form to obtain

$$
(\nabla \mathbf{u})^T \nabla \mathbf{u} - I_3 = A_1 + 2x_3 A_2 + x_3^2 A_3 + \mathcal{O}(x_3^3), \quad (12)
$$

where

$$
A_1 := \begin{bmatrix} \mathbf{I}(\mathbf{y}) - I_2 & \mathbf{0} \\ \mathbf{0} & 0 \end{bmatrix}, \quad A_2 := \begin{bmatrix} -\Pi(\mathbf{y}) & \mathbf{0} \\ \mathbf{0} & \beta \end{bmatrix},
$$

$$
A_3 := \begin{bmatrix} (\nabla' \mathbf{v})^T \nabla' \mathbf{v} + \frac{1}{2} (\nabla' \mathbf{y})^T \nabla' (\beta \mathbf{v}) + \frac{1}{2} (\nabla' (\beta \mathbf{v}))^T \nabla' \mathbf{y} & \frac{1}{2} (\nabla' (\beta \mathbf{v}))^T \mathbf{v} \\ \frac{1}{2} \mathbf{v}^T \nabla' (\beta \mathbf{v}) & \beta^2 \end{bmatrix}.
$$

Using this relation in the energy density (8) with $F = \nabla \mathbf{u}$, we obtain

$$
W(\nabla \mathbf{u}) = \frac{\mu}{4} (|A_1|^2 + 4x_3 A_1 : A_2 + 4x_3^2 |A_2|^2 + 2x_3^2 A_1 : A_3 + \mathcal{O}(x_3^3)) \\
\qquad + \frac{\lambda}{8} (\mathrm{tr}(A_1)^2 + 4x_3 \mathrm{tr}(A_1) \mathrm{tr}(A_2) + 4x_3^2 \mathrm{tr}(A_2)^2 + 2x_3^2 \mathrm{tr}(A_1) \mathrm{tr}(A_3) + \mathcal{O}(x_3^3)).
$$

Note that the terms with odd power of $x_3$ vanish when integrating $x_3$ on $(-\frac{s}{2}, \frac{s}{2})$ to derive an expression of the energy (7). Thus, we get

$$
E_s(\mathbf{u}) = \frac{\mu}{4} \int_{\Omega} \left( s|A_1|^2 + \frac{s^3}{3}|A_2|^2 + \frac{s^3}{6} A_1 : A_3 \right) d\mathbf{x}' \\
\qquad + \frac{\lambda}{8} \int_{\Omega} \left( s \operatorname{tr}(A_1)^2 + \frac{s^3}{3} \operatorname{tr}(A_2)^2 + \frac{s^3}{6} \operatorname{tr}(A_1) \operatorname{tr}(A_3) \right) d\mathbf{x}' + \mathcal{O}(s^4).
$$

The terms with $A_1 : A_3$ and $\operatorname{tr}(A_1)\operatorname{tr}(A_3)$ are of higher order. Indeed, thanks to the scaling assumption (11), we have

$$
\int_{\Omega} A_1 : A_3 dx' \leq \left( \int_{\Omega} |A_1|^2 dx' \right)^{\frac{1}{2}} \left( \int_{\Omega} |A_3|^2 dx' \right)^{\frac{1}{2}} \leq \left( \int_{\Omega} |A_3|^2 dx' \right)^{\frac{1}{2}} \left( \frac{4}{\mu} \Lambda \right)^{\frac{1}{2}} s,
$$

where $(\int_{\Omega} |A_3|^2 dx')^{\frac{1}{2}}$ is independent of $s$. A similar argument holds for the trace terms. As a consequence, in view of the definitions of $A_1$ and $A_2$ in (12), the expression (5) for the first and second fundamental forms, and the definition of $\beta$ in (10), we arrive at the final expression for the three dimensional energy per unit volume $s^{-1} E_s(\mathbf{u})$:

$$
s^{-1} E_s(\mathbf{u}) = \tilde{E}^S(\mathbf{y}) + s^2 \tilde{E}^B(\mathbf{y}) + \mathcal{O}(s^3), \quad (13)
$$

where

$$
\tilde{E}^S(\mathbf{y}) := \frac{1}{8} \int_{\Omega} \left( 2\mu |\mathbf{I}(\mathbf{y}) - I_2|^2 + \lambda \operatorname{tr} (\mathbf{I}(\mathbf{y}) - I_2)^2 \right) d\mathbf{x}'
$$

$$
\tilde{E}^B(\mathbf{y}) := \frac{1}{24} \int_{\Omega} \left( 2\mu |\Pi(\mathbf{y})|^2 + \frac{2\mu\lambda}{2\mu + \lambda} \operatorname{tr}(\Pi(\mathbf{y}))^2 \right) d\mathbf{x}'
$$

denote the stretching and bending energies, respectively.

2.2.2. Limiting Bending Model

We are now interested in taking the limit when $s \to 0^+$. We recall that we are assuming (11), which, as we shall see, implies that the limiting plate deformation $\mathbf{y}$ cannot stretch nor shear the midplane $\Omega$ but can bend it to reduce its energy.

The starting point to derive the limiting energy as the thickness $s$ vanishes is the preasymptotic expression (13) for the energy. The bending regime condition (11) implies that $\tilde{E}^S(\mathbf{y}) = 0$ and

$$
E(\mathbf{y}) := \lim_{s \to 0^+} s^{-3} E_s(\mathbf{u}) = \frac{1}{24} \int_{\Omega} \left( 2\mu |\Pi(\mathbf{y})|^2 + \frac{2\mu\lambda}{2\mu + \lambda} \operatorname{tr}(\Pi(\mathbf{y}))^2 \right) d\mathbf{x}'. \quad (14)
$$
Note that when $\tilde{E}^S(\mathbf{y}) = 0$,

$$I(\mathbf{y}) = (\nabla' \mathbf{y})^T \nabla' \mathbf{y} = I_2 \quad \text{a.e. in } \Omega,$$

i.e. is an isometry and in particular $\mathbf{y} \in [W^{1,\infty}(\Omega)]^3$. Furthermore, for isometries the following relations for the second fundamental form hold (see for instance [37]):

$$|\Pi(\mathbf{y})| = |\nabla' \nabla' \mathbf{y}| = \mathrm{tr}(\Pi(\mathbf{y})). \qquad (15)$$

This allows to rewrite the limiting energy as

$$E(\mathbf{y}) = \frac{\alpha}{2} \int_{\Omega} |\nabla' \nabla' \mathbf{y}|^2 d\mathbf{x}', \qquad (16)$$

when $\mathbf{y} \in [H^2(\Omega) \cap W^{1,\infty}(\Omega)]^3$ is an isometry, i.e. $I(\mathbf{y}) = I_2$, and where $\alpha := \frac{\mu(\mu+\lambda)}{3(2\mu+\lambda)}$.

## 2.3. Single Layer with a General Metric Constraint

Without additional constraints such as boundary conditions, the minimum of energy (7) with density (8) is achieved by the identity deformation $\mathbf{u}(\mathbf{x}) = \mathbf{x}$, which is an isometry. In particular, the reference flat configuration is stress-free. We now consider pretrained materials characterized by the presence of internal stresses in the flat configuration. Examples are nematic glasses [45, 46], natural growth of soft tissues [3, 4], and manufactured polymer gels [7, 47, 48].

We follow [20, 24] (see also [27]) and modify the energy (7) to read

$$E_s^{\text{pre}}(\mathbf{u}) := \int_{\Omega_s} W \left( \nabla \mathbf{u}(\mathbf{x}) G(\mathbf{x})^{-\frac{1}{2}} \right) d\mathbf{x}, \qquad (17)$$

where $G: \Omega_s \to \mathbb{R}^{3 \times 3}$ is a given symmetric positive definite target metric and $G^{-\frac{1}{2}}$ denotes the positive definite symmetric square root of the inverse $G^{-1}$ of $G$. Note that the previous case is recovered when $G = I_3$. Furthermore, when $G$ is a metric immersion, i.e. there exists a deformation $\mathbf{u}: \Omega_s \to \mathbb{R}^3$ such that

$$(\nabla \mathbf{u})^T \nabla \mathbf{u} = G \quad \text{a.e. in } \Omega_s, \qquad (18)$$

we have

$$(\nabla \mathbf{u}(\mathbf{x}) G(\mathbf{x})^{-\frac{1}{2}})^T \nabla \mathbf{u}(\mathbf{x}) G(\mathbf{x})^{-\frac{1}{2}} = I_3,$$

and thus $\mathbf{u}$ is a minimizer of $E_s^{\text{pre}}$ with $E_s^{\text{pre}}(\mathbf{u}) = 0$.

Following [20, 23], we further assume that $G$ has the form

$$G(\mathbf{x}', x_3) = \begin{pmatrix} g(\mathbf{x}') & \mathbf{0} \\ \mathbf{0} & 1 \end{pmatrix}, \quad \mathbf{x}' \in \Omega, x_3 \in (-s/2, s/2), \qquad (19)$$

where $g: \Omega \to \mathbb{R}^{2 \times 2}$ is symmetric and uniformly positive definite. In other words, the metric $G$ is independent of the variable $x_3$, is uniform throughout the thickness, and no stretching is allowed in the vertical direction.

The $\Gamma$-convergence of this model towards the two dimensional energy given below is obtained in [20, 24]. Proceeding as in Section 2.2.1, an intuitive derivation of the two dimensional $\Gamma$-limit can be obtained, using again assumption (10). Indeed, when $0 < s \ll 1$, neglecting the $\mathcal{O}(s^3)$ terms the three dimensional prestrain energy, (17) satisfies

$$s^{-1} E_s^{\text{pre}}(\mathbf{u}) = E^S(\mathbf{y}) + s^2 E^B(\mathbf{y}) := E_s^{\text{pre}}(\mathbf{y}), \qquad (20)$$

where

$$E^S(\mathbf{y}) := \frac{1}{8} \int_{\Omega} \left( 2\mu \left| g^{-\frac{1}{2}}(\mathbf{I}(\mathbf{y}) - \mathbf{g}) g^{-\frac{1}{2}} \right|^2 + \lambda \mathrm{tr} \left( g^{-\frac{1}{2}}(\mathbf{I}(\mathbf{y}) - \mathbf{g}) g^{-\frac{1}{2}} \right)^2 \right) d\mathbf{x}' \qquad (21)$$
and

$$E^B(\mathbf{y}) := \frac{1}{24} \int_{\Omega} \left( 2\mu \left| g^{-\frac{1}{2}} \Pi(\mathbf{y}) g^{-\frac{1}{2}} \right|^2 + \frac{2\mu\lambda}{2\mu + \lambda} \operatorname{tr} \left( g^{-\frac{1}{2}} \Pi(\mathbf{y}) g^{-\frac{1}{2}} \right)^2 \right) \mathrm{d}\mathbf{x}'. \quad (22)$$

Note that we slightly abuse the notation by using $E_s^{\text{pre}}$ to refer to both the three dimensional and two dimensional prestrain energy. Which energy $E_s^{\text{pre}}$ is referring to will be clear from the context.

For the limit $s^{-3}E_s(\mathbf{u})$ to be finite as the thickness $s$ goes to zero, the two dimensional deformation $\mathbf{y}: \Omega \to \mathbb{R}^3$ must satisfy the metric constraint

$$(\nabla \mathbf{y})^T \nabla \mathbf{y} = g \quad (23)$$

and in that case

$$\lim_{s \to 0^+} s^{-3} E_s(\mathbf{u}) = \frac{1}{24} \int_{\Omega} \left( 2\mu \left| g^{-\frac{1}{2}} \Pi(\mathbf{y}) g^{-\frac{1}{2}} \right|^2 + \frac{2\mu\lambda}{2\mu + \lambda} \operatorname{tr} \left( g^{-\frac{1}{2}} \Pi(\mathbf{y}) g^{-\frac{1}{2}} \right)^2 \right) \mathrm{d}\mathbf{x}'. \quad (24)$$

**Remark 1.** Nash's theorem guarantees that there exists an isometric immersion of $g$ into $\mathbb{R}^{10}$ [49], but whether or not there exists an isometric immersion into $\mathbb{R}^3$ depends on $g$. We refer to [40] for a discussion of positive and negative results for metrics with specific properties. Note that the existence of $\mathbf{y} \in H^2(\Omega)$ satisfying (23) is equivalent to the boundedness condition (11) on $E_s$. We refer to [20] for the case where $G$ has the form (19) and to [50] for a general Riemannian metric.

Comparing the limiting energy (24) with (14), we find that the isometry case is recovered for $g = I_2$. However, in the case of isometries, the energy is further reduced to (16) thanks to (15). When $g \neq I_2$, the second fundamental form of $\mathbf{y}$ in (24) cannot be substituted by the Hessian of $\mathbf{y}$. Nevertheless, [27, Proposition 1] and [33, Proposition A.1] guarantee that for $\mathbf{y} \in [H^2(\Omega)]^3$ and $g \in [H^1(\Omega) \cap L^\infty(\Omega)]^{2 \times 2}$, up to an additive term only depending on $g$, the limiting energy (24) is equal to

$$\frac{1}{24} \sum_{m=1}^{3} \int_{\Omega} \left( 2\mu \left| g^{-\frac{1}{2}} (\nabla' \nabla' y_m) g^{-\frac{1}{2}} \right|^2 + \frac{2\mu\lambda}{2\mu + \lambda} \operatorname{tr} \left( g^{-\frac{1}{2}} (\nabla' \nabla' y_m) g^{-\frac{1}{2}} \right)^2 \right) \mathrm{d}\mathbf{x}',$$

where $y_m$ is the $m$th component of $\mathbf{y}$. In addition, in view of the metric constraint (23) we have $\mathbf{y} \in [W^{1,\infty}(\Omega)]^3$. As a consequence, the equilibrium deformations of the thin limit of pretrained materials are deformations $\mathbf{y} \in [H^2(\Omega) \cap W^{1,\infty}(\Omega)]^3$ satisfying $I(\mathbf{y}) = g$ and minimizing

$$E^{\text{pre}}(\mathbf{y}) := \frac{1}{24} \sum_{m=1}^{3} \int_{\Omega} \left( 2\mu \left| g^{-\frac{1}{2}} (\nabla' \nabla' y_m) g^{-\frac{1}{2}} \right|^2 + \frac{2\mu\lambda}{2\mu + \lambda} \operatorname{tr} \left( g^{-\frac{1}{2}} (\nabla' \nabla' y_m) g^{-\frac{1}{2}} \right)^2 \right) \mathrm{d}\mathbf{x}'. \quad (25)$$

## 2.4. Bilayer Plates

We now discuss an anisotropic model where the three dimensional hyperelastic structure is a compound of two thin layers with different mechanical properties. Each material exhibits a different response under external stimuli, generating a material mismatch compensated by the bending of the structure. Temperature, humidity, pH, and electric current are typical external stimuli. We refer for instance to [9, 51–53] for devices actuated by electric current and to [11, 54] for humidity controlled materials.

From a mathematical point of view, this corresponds to a prestrain tensor of the form

$$G(\mathbf{x}', x_3) = (I_3 \pm \zeta s N)^T (I_3 \pm \zeta s N) = I_3 \pm 2\zeta s N + \zeta^2 s^2 N^2, \quad \pm x_3 \ge 0, \quad (26)$$

where the inhomogeneity and anisotropy of the bilayer material is encoded in the tensor

$$N := \begin{pmatrix} N_{\mathbf{x}'\mathbf{x}} & \mathbf{m} \\ \mathbf{m}^T & n \end{pmatrix}$$

with $N_{\mathbf{x}'\mathbf{x}}: \Omega \to \mathbb{R}^{2 \times 2}$ uniformly symmetric, $n: \Omega \to \mathbb{R}$, and $\mathbf{m}: \Omega \to \mathbb{R}^2$. Note that the expression (26) does not reduce to the previous setting with target metric (18).
The limiting model is derived in [28] via Γ-convergence when $N = I_3$ and can again be formally recovered using a modified Kirchhoff-Love assumption, namely

$$ \mathbf{u}(\mathbf{x}', x_3) = \mathbf{y}(\mathbf{x}') + x_3(1 \pm c\zeta s)\mathbf{v}(\mathbf{x}') + \frac{1}{2}x_3^2 d(\mathbf{x}')\mathbf{v}(\mathbf{x}') $$

with

$$ c := \frac{2\mu + 3\lambda}{2\mu + \lambda} \quad \text{and} \quad d(\mathbf{x}') := -\frac{\lambda}{2\mu + \lambda} \operatorname{tr}(\Pi(\mathbf{x}')). $$

Proceeding as in the previous section, the limiting energy as the thickness $s$ goes to zero is

$$ \frac{1}{24} \int_{\Omega} \left( 2\mu |\Pi(\mathbf{y}) - 3\zeta N_{\mathbf{x}'\mathbf{x}'}|^2 + \frac{2\mu\lambda}{2\mu + \lambda} \operatorname{tr}(\Pi(\mathbf{y}) - 3\zeta N_{\mathbf{x}'\mathbf{x}'}')^2 \right) d\mathbf{x}' + C $$

provided that $\mathbf{y}$ is an isometry, where the constant $C$ depends on $\mu$, $\lambda$, $\zeta$, and $N$, but not on the deformation $\mathbf{y}$. We refer to [29] for more details in the case $\mu = 6$ and $\lambda = 0$.

In what follows, we consider an equivalent limiting two dimensional energy

$$ \frac{\alpha}{2} \int_{\Omega} |\Pi(\mathbf{y}) - Z|^2 d\mathbf{x}' $$

for $\alpha > 0$ and with $Z(\mathbf{x}') := 3\zeta N_{\mathbf{x}'\mathbf{x}'}$. Note that $Z$ acts as an intrinsic spontaneous curvature tensor favoring deformations with principal curvatures matching the eigenvalues of $Z$. Now recall that the entries of the second fundamental form $\Pi(\mathbf{y})$ are given by

$$ [\Pi(\mathbf{y})]_{ij} = \partial_{ij}\mathbf{y} \cdot \mathbf{v} = \partial_{ij}\mathbf{y} \cdot \frac{\partial_1 \mathbf{y} \times \partial_2 \mathbf{y}}{|\partial_1 \mathbf{y} \times \partial_2 \mathbf{y}|}, \quad i, j = 1, 2, $$

and is thus nonlinear in $\mathbf{y}$. Taking advantage of relation (15) valid for isometries, we have

$$ |\Pi(\mathbf{y})| = |\nabla'\nabla'\mathbf{y}| \quad \text{and} \quad |\partial_1 \mathbf{y} \times \partial_2 \mathbf{y}| = 1, $$

which allows us to write

$$ \frac{\alpha}{2} \int_{\Omega} |\Pi(\mathbf{y}) - Z|^2 d\mathbf{x}' = \frac{\alpha}{2} \int_{\Omega} |\nabla'\nabla'\mathbf{y}|^2 d\mathbf{x}' - \alpha \sum_{i,j=1}^{2} \int_{\Omega} \partial_{ij}\mathbf{y} \cdot (\partial_1\mathbf{y} \times \partial_2\mathbf{y}) Z_{ij} d\mathbf{x}' + \frac{\alpha}{2} \int_{\Omega} |Z|^2 d\mathbf{x}'. $$

Deformations of bilayer materials are thus deformations $\mathbf{y} \in [H^2(\Omega) \cap W^{1,\infty}(\Omega)]^3$ satisfying $\Pi(\mathbf{y}) = I_2$ and minimizing

$$ E^{\text{bil}}(\mathbf{y}) := \frac{\alpha}{2} \int_{\Omega} |\nabla'\nabla'\mathbf{y}|^2 d\mathbf{x}' - \alpha \sum_{i,j=1}^{2} \int_{\Omega} \partial_{ij}\mathbf{y} \cdot (\partial_1\mathbf{y} \times \partial_2\mathbf{y}) Z_{ij} d\mathbf{x}'. \quad (27) $$

## 2.5. Folding

We incorporate the capability for the plates to fold along a given crease. In particular, we make structural assumptions leading to a limiting model for which folding does not require any energy. We restrict our description to single layer plates with isometry constraint (i.e. $g = I_2$). The derivation of the folding model described here was originally proposed in [30]. Extensions to other models, although feasible, have not been treated in the existing literature. However, numerical simulation are available in [31] for the bending of a bilayer plate with folding.

Folding models are obtained by assuming that the thin structure has a material defect in a cylindrical neighborhood of a $C^2$ curve $\Sigma \subset \Omega$ (the crease) splitting $\Omega$ in two parts and intersecting the boundary transversely. Let $\Omega_s = \Omega \times (-s/2, s/2)$ be as above and let $\Sigma_{s,r} := B(\Sigma, r) \times (-s/2, s/2)$ be the location of the defect. Here, for $r > 0$ to be determined, $B(\Sigma, r) := \cup_{\mathbf{x}' \in \Sigma} \{\mathbf{z}' \in \Omega : |\mathbf{z}' - \mathbf{x}'| < r\}$ (see Figure 2).

The hyperelastic energy (7) in presence of a defect reads

$$ E_s^{\text{fold}}(\mathbf{u}) := \int_{\Omega_s} f_{\varepsilon,r}(\mathbf{x}') W(\nabla \mathbf{u}(\mathbf{x})) d\mathbf{x}, $$
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Material with defect. (Left) Midplane domain $\Omega$ with a crease $\Sigma$ and its tubular neighborhood $B(\Sigma, r)$. (Right) Thin domain $\Omega_s$ with inflated crease $\Sigma_{s,r}$.</figcaption></figure>

where $f_{\varepsilon,r} : \Omega \to (0,1]$ defined by

$$f_{\varepsilon,r}(\mathbf{x}') := \varepsilon\chi_{B(\Sigma,r)}(\mathbf{x}') + 1 - \chi_{B(\Sigma,r)}(\mathbf{x}')$$

is 1 away from $B(\Sigma, r)$ and $\varepsilon$ in $B(\Sigma, r)$ for $0 < \varepsilon \ll 1$. Here $\chi_{B(\Sigma,r)}$ is the characteristic function of $B(\Sigma, r)$.

We assume the following asymptotic behaviors of $r$ and $\varepsilon$:

$$\limsup_{s \to 0^+} \frac{s^2}{\varepsilon} < \infty, \quad \limsup_{s \to 0^+} \frac{s}{r} < \infty, \quad \text{and} \quad \limsup_{s \to 0^+} \frac{\varepsilon r}{s^2} = 0. \qquad (28)$$

They guarantee that the limiting deformations are globally Lipschitz and that folding is performed at no energy cost. The $\Gamma$-limit of $s^{-3}E_s^{\text{fold}}$ under assumptions (28) and in the bending regime $s^{-3}E_s^{\text{fold}} \le \Lambda < \infty$, is

$$E^{\text{fold}}(\mathbf{y}) := \frac{1}{24} \int_{\Omega \setminus \Sigma} \left( 2\mu |\Pi(\mathbf{y})|^2 + \frac{2\mu\lambda}{2\mu + \lambda} \text{tr}(\Pi(\mathbf{y}))^2 \right) d\mathbf{x}' = \frac{\alpha}{2} \int_{\Omega \setminus \Sigma} |\nabla' \nabla' \mathbf{y}|^2 d\mathbf{x}',$$

provided $\mathbf{y} \in [H^2(\Omega \setminus \Sigma) \cap W^{1,\infty}(\Omega)]^3$ is an isometry and where $\alpha = \alpha(\mu, \lambda)$ is defined as in (16); see [30]. In particular, $E^{\text{fold}}$ reduces to $E$ in (16) when no creases are present.

Although not rigorously justified, in the numerical section below we provide experiments with multiple folding arcs not necessarily intersecting the boundary. Also, we assume that the prestrain and bilayer energies can be modified similarly to account for possible folding along a crease $\Sigma$. Without modifying the notation but with the understanding that $\Sigma = \emptyset$ in the original cases (25) and (27), we redefine for $\mathbf{y} \in [H^2(\Omega \setminus \Sigma) \cap W^{1,\infty}(\Omega)]^3$

$$E^{\text{pre}}(\mathbf{y}) := \frac{1}{24} \sum_{m=1}^{3} \int_{\Omega \setminus \Sigma} \left( 2\mu \left|g^{-\frac{1}{2}}(\nabla' \nabla' y_m)g^{-\frac{1}{2}}\right|^2 + \frac{2\mu\lambda}{2\mu + \lambda} \text{tr}\left(g^{-\frac{1}{2}}(\nabla' \nabla' y_m)g^{-\frac{1}{2}}\right)^2 \right) d\mathbf{x}' \quad (29)$$

and

$$E^{\text{bil}}(\mathbf{y}) := \frac{\alpha}{2} \int_{\Omega \setminus \Sigma} |\nabla' \nabla' \mathbf{y}|^2 d\mathbf{x}' - \alpha \sum_{i,j=1}^{2} \int_{\Omega \setminus \Sigma} \partial_{ij} \mathbf{y} \cdot (\partial_1 \mathbf{y} \times \partial_2 \mathbf{y}) Z_{ij} d\mathbf{x}'. \quad (30)$$

## 2.6. Forcing Term and Boundary Conditions

To complete the model, it remains to discuss how to incorporate external forces and boundary conditions.

In the presence of an external force $\mathbf{f}_s : \Omega_s \to \mathbb{R}^3$, the hyperelastic energy (7) must be modified to read

$$E_s(\mathbf{u}) = \int_{\Omega_s} W(\nabla \mathbf{u}(\mathbf{x})) d\mathbf{x} - \int_{\Omega_s} \mathbf{f}_s(\mathbf{x}) \cdot \mathbf{u}(\mathbf{x}) d\mathbf{x}.$$
In view of (11), the force applied to the three dimensional thin structure $\Omega_s = \Omega \times (-s/2, s/2)$ must satisfy

$$ \left| \lim_{s \to 0^+} s^{-3} \int_{-s/2}^{s/2} \mathbf{f}_s(\mathbf{x}', x_3) d\mathbf{x}_3 \right| < \infty, \quad \mathbf{x}' \in \Omega. $$

This allows us to define $\mathbf{f}(\mathbf{x}') := \lim_{s \to 0^+} s^{-3} \int_{-s/2}^{s/2} \mathbf{f}_s(\mathbf{x}', x_3) d\mathbf{x}_3$ which yields

$$ \lim_{s \to 0^+} s^{-3} \int_{\Omega_s} \mathbf{f}_s(\mathbf{x}) \cdot \mathbf{u}(\mathbf{x}) d\mathbf{x} = \int_{\Omega} \mathbf{f}(\mathbf{x}') \cdot \mathbf{y}(\mathbf{x}') d\mathbf{x}' $$

when using the modified Kirchhoff-Love ansatz (10) on the three dimensional deformation $\mathbf{u}$.

Therefore, when external forces are present, the term

$$ F(\mathbf{y}) := \int_{\Omega} \mathbf{f}(\mathbf{x}') \cdot \mathbf{y}(\mathbf{x}') d\mathbf{x}' \quad (31) $$

must be subtracted from the energies $E_s^{\text{pre}}(\mathbf{y})$, $E^{\text{pre}}(\mathbf{y})$, and $E^{\text{bil}}(\mathbf{y})$ (scaled by $s^2$ for the preasymptotic energy).

Different types of boundary conditions can be incorporated into the system. We say that *Dirichlet boundary conditions* are imposed on $\Gamma^D \subset \partial\Omega$ when

$$ \mathbf{y} = \boldsymbol{\varphi} \quad \text{and} \quad \nabla' \mathbf{y} = \Phi \quad \text{on} \quad \Gamma^D, \quad (32) $$

where $\boldsymbol{\varphi}: \Omega \to \mathbb{R}^3$ and $\Phi: \Omega \to \mathbb{R}^{3 \times 2}$ are sufficiently smooth, and with $\Phi$ satisfying the compatibility condition $\Phi^T \Phi = g$ a.e. in $\Omega$. In the case where $\boldsymbol{\varphi}(\mathbf{x}) = (\mathbf{x}, 0)^T$ and $\Phi = \nabla \boldsymbol{\varphi}$, the boundary conditions (32) are referred to as *clamped boundary conditions*. A *mixed boundary condition* occurs when only the value of the deformation on $\Gamma^M \subset \partial\Omega$ is specified, that is

$$ \mathbf{y} = \boldsymbol{\varphi} \quad \text{on} \quad \Gamma^M \quad (33) $$

for some smooth function $\boldsymbol{\varphi}: \Omega \to \mathbb{R}^3$.

All the energy models are defined on $[W^{1,\infty}(\Omega)]^3$ and thus also on continuous functions $[C^0(\bar{\Omega})]^3$. This means that *pointwise boundary conditions* can be enforced on the deformations as well. Given a collection of points $\Gamma^P := \{\mathbf{x}_1, \mathbf{x}_2, \dots, \mathbf{x}_P\}$ lying on $\partial\Omega$ and corresponding values $\{\boldsymbol{\varphi}_1, \boldsymbol{\varphi}_2, \dots, \boldsymbol{\varphi}_P\} \subset \mathbb{R}^3$, we set

$$ \mathbf{y}(\mathbf{x}_i) = \boldsymbol{\varphi}_i \quad \forall \mathbf{x}_i \in \Gamma^P. $$

The *free boundary case* refers to the case where $\Gamma^D \cup \Gamma^M \cup \Gamma^P = \emptyset$.

### 3. Space Discretization

We now describe numerical methods based on finite element methods to approximate minimizers of the preasymptotic energy $E_s^{\text{pre}}$ in (20) and the limiting energies $E^{\text{pre}}$ in (29) and $E^{\text{bil}}$ in (30). We point out two difficulties addressed below: (i) all the energies are defined on a subspace of $[H^2(\Omega \setminus \Sigma)]^3$ (recall that $\Sigma = \emptyset$ when no creases are present) and thus conforming discretizations can only be achieved by globally $C^0$ elements which are also $C^1$ on $\Omega \setminus \Sigma$; (ii) all but the preasymptotic energy require the deformations to satisfy the metric constraint (23). The latter is too rigid to be satisfied a.e. in $\Omega$ when using piecewise polynomial approximations.

Discontinuous Galerkin finite elements are put forward to circumvent the conformity requirement. In essence, discrete energies are designed by substituting the Hessian of the deformation by a suitable discrete Hessian. It is worth pointing out that discontinuous Galerkin elements are also preferred for their flexibility; in particular, they are more efficient for the approximation of the metric constraint [34, 35].

We first set up the notation related to discontinuous Galerkin finite element methods in Section 3.1. Then, in Section 3.2, we introduce a discrete Hessian operator, the crucial ingredient in all the numerical schemes discussed in Section 3.3.
To simplify the notation, from now on we will write **x** instead of **x'**, and accordingly for the differential operators (e.g. $\nabla'$ will be denoted $\nabla$). Moreover, the variable of integration will no longer be indicated.

### 3.1. Discontinuous Galerkin Finite Elements

For the description of the method we assume that the midplane $\Omega \subset \mathbb{R}^2$ is a polygonal domain, but extensions to general domains can be performed using standard techniques. Let $\{\mathcal{T}_h\}_{h>0}$ be a sequence of shape-regular but possibly graded partitions of $\Omega$ made of elements $T$, either triangles or quadrilaterals, of diameter at most $h$. Hanging nodes are thus allowed, and we assume that all the elements within each domain of influence have comparable diameters [55]. Let $\mathcal{E}_h := \mathcal{E}_h^0 \cup \mathcal{E}_h^b$ denote the set of edges, where $\mathcal{E}_h^0$ stands for the set of interior edges and $\mathcal{E}_h^b$ for the set of boundary edges. We also define $\mathcal{E}_h^D$ to be the set of Dirichlet boundary edges and $\mathcal{E}_h^M$ to be the set of mixed boundary edges.

In the presence of a crease $\Sigma$, we assume that the crease is well approximated by the subdivision $\mathcal{T}_h$ in the spirit of [30]. That is, the folding curve $\Sigma$ is approximated by the piecewise linear curve

$$ \Sigma_h := \cup_{l=1}^{L} e_l, \quad (34) $$

where $\mathcal{E}_h^\Sigma := \{e_l\}_{l=1}^L \subset \mathcal{E}_h^0$ is a collection of interior edges such that the endpoints of each edge $e_l$ belong to $\Sigma$.

It will be convenient later to group all the active edges into two sets: $\mathcal{E}_h^{\text{val}} := \mathcal{E}_h^0 \cup \mathcal{E}_h^D \cup \mathcal{E}_h^M$ for the deformations and $\mathcal{E}_h^{\text{grad}} := (\mathcal{E}_h^0 \setminus \mathcal{E}_h^\Sigma) \cup \mathcal{E}_h^D$ for the gradient of the deformations. These sets contain all the edges over which averages and jumps will be computed. Notice that in the presence of a crease, the crease edges are not included in the active edge set for the gradients to allow for discontinuous gradients. As we shall see, this is the only difference between the schemes with and without folding.

We define the diameter function $h$ on $\mathcal{T}_h \cup \mathcal{E}_h \cup \Gamma^P$ (recall that $\Gamma^P = \{\mathbf{x}_i\}_{i=1}^P$ are the locations where pointwise boundary conditions are prescribed) by

$$ h_{|T} := h_T := \text{diam}(T), \quad \forall T \in \mathcal{T}_h, \quad h_{|e} := h_e := \text{diam}(e), \quad \forall e \in \mathcal{E}_h, $$

$$ h_{|\mathbf{x}_i} := h_i := \frac{1}{\#(\omega_i)} \sum_{T \in \omega_i} h_T, \quad i = 1, 2, \dots, P, $$

where $\omega_i := \{T \in \mathcal{T}_h : \mathbf{x}_i \in \partial T\}$. The minimum mesh size is $h_{\min} := \min_{T \in \mathcal{T}_h} h_T$.

We assume that the boundary regions $\Gamma^D$ and $\Gamma^M$, where Dirichlet and mixed boundary conditions are enforced, are exactly captured by all the subdivisions. This means that for all $h > 0$ we have

$$ \Gamma^D = \{e: e \in \mathcal{E}_h^D\} \quad \text{and} \quad \Gamma^M = \{e: e \in \mathcal{E}_h^M\} $$

for some $\mathcal{E}_h^D, \mathcal{E}_h^M \subset \mathcal{E}_h^b$. We also assume that the locations $\Gamma^P = \{\mathbf{x}_i\}_{i=1}^P$, where pointwise boundary conditions are prescribed, are vertices of all the subdivisions. To keep a uniform notation, we introduce various skeletons of the mesh $\mathcal{T}_h$

$$ \Gamma_h^a = \cup \{e: e \in \mathcal{E}_h^a\} \quad \text{for } a \in \{0, D, M, \text{val, grad}\}. $$

For an integer $r \ge 0$, let $\mathbb{P}_r$ and $\mathbb{Q}_r$ denote polynomials of total degree at most $r$ and degree at most $r$ in each variable, respectively. We then define the finite element space

$$ \mathbb{V}_h^r := \{v_h \in L^2(\Omega) : v_{h|T} \circ F_T \in \mathbb{P}_r \quad \forall T \in \mathcal{T}_h\}, \quad (35) $$

where $\mathbb{P}_r$ is replaced by $\mathbb{Q}_r$ if the subdivision is made of quadrilaterals. Here $F_T : \hat{T} \to T$ is the map from the reference element $\hat{T}$ (unit simplex for triangulation and unit square for subdivisions
made of quadrilaterals) to the physical element $T \in \mathcal{T}_h$. In what follows, the deformation of the plate y will be approximated by a discrete function $y_h \in [\mathbb{V}_h^k]^3$ with $k \ge 2$.

The broken (elementwise) gradient of a scalar function $v_h \in \mathbb{V}_h^k$ is denoted $\nabla_h v_h$. We use a similar notation for other differential operators, for instance $D_h^2 v_h = \nabla_h \nabla_h v_h$ denotes the broken Hessian. Moreover, these operators are applied componentwise for vector-valued functions.

Let us now introduce the jump and average operators instrumental for discontinuous Galerkin methods to relax conformity conditions. For each $e \in \mathcal{E}_h^0$ let $\mathbf{n}_e$ be a unit normal to $e$ (the orientation is arbitrary but fixed once for all). With this notation, we set for $v_h \in \mathbb{V}_h^r$ and $e \in \mathcal{E}_h^0$

$$[v_h]|_e := v_h^{-} - v_h^{+}, \quad \{v_h}\!|_e := \frac{v_h^{+} + v_h^{-}}{2},$$

where $v_h^{\pm}(\mathbf{x}) := \lim_{\rho \to 0^+} v_h(\mathbf{x} \pm \rho \mathbf{n}_e)$ for any $\mathbf{x} \in e$. Note that if $e = \partial T_1 \cap \partial T_2$ with $T_1, T_2 \in \mathcal{T}_h$, then

$$[v_h]|_e \mathbf{n}_e = v_{h|T_1} \mathbf{n}_{\partial T_1} + v_{h|T_2} \mathbf{n}_{\partial T_2}$$

with $\mathbf{n}_{\partial T_i}$ the outward unit normal to $\partial T_i$, $i = 1, 2$. For edges on the boundary $e \in \mathcal{E}_h^b$, $\mathbf{n}_e$ is the outward unit normal vector to $\partial \Omega$ and $[v_h]|_e = \{v_h}\!|_e := v_h$. The jumps and averages of non-scalar functions are computed componentwise.

We introduce the mesh-dependent bilinear form $\langle \cdot, \cdot \rangle_{H_h^2(\Omega)}$ defined for any $v_h, w_h \in \mathbb{V}_h^k$ by

$$\langle v_h, w_h \rangle_{H_h^2(\Omega)} := (D_h^2 v_h, D_h^2 w_h)_{L^2(\Omega)} + (\mathbf{h}^{-1}[\nabla_h v_h], [\nabla_h w_h])_{L^2(\Gamma_h^{\text{grad}})} + (\mathbf{h}^{-3}[v_h], [w_h])_{L^2(\Gamma_h^{\text{val}})} + \sum_{\mathbf{x}_i \in \Gamma^P} h_i^{-2} v_h(\mathbf{x}_i) w_h(\mathbf{x}_i) \quad (36)$$

and we set

$$|v_h|_{H_h^2(\Omega)} := \langle v_h, v_h \rangle_{H_h^2(\Omega)}^{1/2}.$$

These definitions are extended componentwise to vector-valued functions.

## 3.2. Discrete Hessian

All the energies considered in this work involve the Hessian of the deformation. Conforming approximations come at the price of using costly and rigid $C^1$ finite element methods. Instead, we describe an LDG approach, which retains the simplicity of the discrete problem while not requiring $C^1$ elements.

The proposed LDG methods consist of replacing the continuous Hessian $D^2y$ by a suitable discrete (aka reconstructed) Hessian $H_h(y_h)$ that we introduce now. The discrete Hessian described here was originally proposed in [27] following ideas from, for example, [56, 57] but modified to possibly account for subdivision containing non-affine elements.

For $\mathbf{v}_h \in [\mathbb{V}_h^k]^3$, we recall that $D_h^2 \mathbf{v}_h : \Omega \to \mathbb{R}^{3 \times 2 \times 2}$ is the piecewise Hessian of $\mathbf{v}_h$. Because $\mathbf{v}_h$ does not belong to $[H^2(\Omega)]^3$, a better approximation of the Hessian must account for the jumps of $\mathbf{v}_h$ and $\nabla_h \mathbf{v}_h$ across elements. To achieve this, we use (local) lifting operators [58–60] $r_e : [L^2(e)]^2 \to [\mathbb{V}_h^{l_1}]^{2 \times 2}$ and $b_e : L^2(e) \to [\mathbb{V}_h^{l_2}]^{2 \times 2}$ to construct representations of these jumps in $L^2(\Omega)$. Given $l_1, l_2 \ge 0$, they are defined for $e \in \mathcal{E}_h$ by

$$r_e(\boldsymbol{\phi}) \in [\mathbb{V}_h^{l_1}]^{2 \times 2} : \int_{\Omega} r_e(\boldsymbol{\phi}) : \tau_h = \int_e \{\tau_h\} \mathbf{n}_e \cdot \boldsymbol{\phi} \quad \forall \tau_h \in [\mathbb{V}_h^{l_1}]^{2 \times 2},$$

and

$$b_e(\boldsymbol{\phi}) \in [\mathbb{V}_h^{l_2}]^{2 \times 2} : \int_{\Omega} b_e(\boldsymbol{\phi}) : \tau_h = \int_e \{\text{div} \, \tau_h\} \cdot \mathbf{n}_e \boldsymbol{\phi} \quad \forall \tau_h \in [\mathbb{V}_h^{l_2}]^{2 \times 2}.$$

Note that the support of the liftings $r_e$ and $b_e$ is the union of the two elements sharing $e$ as an edge (which reduces to a single element if $e$ is a boundary edge).
With these lifting operators, we can now define the discrete Hessian of a function $v_h \in \mathbb{V}_h^k$ as

$$ 
\begin{aligned}
H_h(v_h) &:= D_h^2(v_h) - R_h(|\nabla_h v_h|) + B_h([v_h]) \\
&:= D_h^2(v_h) - \sum_{e \in \mathcal{E}_h^{\text{grad}}} r_e(|\nabla_h v_h|) + \sum_{e \in \mathcal{E}_h^{\text{val}}} b_e([v_h]).
\end{aligned} 
\quad (37) $$

This definition is naturally extended to vector-valued functions in $[\mathbb{V}_h^k]^3$ by componentwise application. For the remainder of this work, we set $l_1 = l_2 = k$.

To motivate definition (37) of the discrete Hessian, we briefly sketch an argument showing a (weak) convergence property of $H_h^2(v_h)$. Let $\tau \in [C_0^\infty(\Omega)]^{2 \times 2}$ and let $\{v_h\}_{h>0} \subset \mathbb{V}_h^k$ be a sequence such that $|v_h|_{H_h^2(\Omega)} \le C$ for some constant $C$ independent of $h$ and such that $v_h \to v$ in $L^2(\Omega)$ as $h \to 0^+$ for some $v \in H^2(\Omega)$. On the one hand, notice that thanks to the regularity of $v$ and the convergence of $v_h$ to $v$, using two integrations by parts we have

$$ \int_{\Omega} D^2 v : \tau = \int_{\Omega} v \operatorname{div}(\operatorname{div}(\tau)) - \int_{\Omega} v_h \operatorname{div}(\operatorname{div}(\tau)) \quad \text{as } h \to 0^+. $$

On the other hand, using again two integrations by parts, we have for any $T \in \mathcal{T}_h$

$$ 
\begin{aligned}
\int_T v_h \operatorname{div}(\operatorname{div}(\tau)) &= \sum_{i,j=1}^2 \int_T v_h \partial_i (\partial_j \tau_{ij}) \\
&= \sum_{i,j=1}^2 \left[ \int_T \partial_j (\partial_i v_h) \tau_{ij} - \int_{\partial T} \partial_i v_h \tau_{ij} (\mathbf{n}_{\partial T})_j + \int_{\partial T} v_h \partial_j \tau_{ij} (\mathbf{n}_{\partial T})_i \right] \\
&= \int_T D^2 v_h : \tau - \int_{\partial T} \nabla v_h \cdot (\tau \mathbf{n}_{\partial T}) + \int_{\partial T} v_h \operatorname{div}(\tau) \cdot \mathbf{n}_{\partial T}.
\end{aligned} $$

Therefore, summing over the elements $T \in \mathcal{T}_h$ and using the fact that $\tau$ is smooth, namely that $[\tau] = \{\tau\} = \tau$ and similarly for $\operatorname{div}(\tau)$, we get

$$ \int_{\Omega} v_h \operatorname{div}(\operatorname{div}(\tau)) = \int_{\Omega} D_h^2 v_h : \tau - \sum_{e \in \mathcal{E}_h^{\text{grad}}} \int_{e} [\nabla_h v_h] \cdot \{\tau\} \mathbf{n}_e + \sum_{e \in \mathcal{E}_h^{\text{val}}} \int_{e} [v_h] \{\operatorname{div}(\tau)\} \cdot \mathbf{n}_e. $$

If $\tau$ were an admissible test function for the lifting operators, then the right hand side of the last relation would be $\int_{\Omega} H_h(v_h) : \tau$ and we would conclude that for all $\tau \in [C_0^\infty(\Omega)]^{2 \times 2}$

$$ \int_{\Omega} D^2 v : \tau \leftarrow \int_{\Omega} H_h(v_h) : \tau \quad \text{as } h \to 0^+, $$

i.e. $H_h(v_h)$ converges weakly to $D^2 v$ in $L^2(\Omega)^{2 \times 2}$. Because $\tau$ does not belong to $[\mathbb{V}_h^k]^{2 \times 2}$ in general, it remains to show that for the Lagrange interpolant $I_h : [C_0^\infty(\Omega)]^{2 \times 2} \to [\mathbb{V}_h^k \cap H_0^1(\Omega)]^{2 \times 2}$ we have

$$ \int_{\Omega} D_h^2 v_h : (\tau - I_h \tau) - \sum_{e \in \mathcal{E}_h^{\text{grad}}} \int_{e} [\nabla_h v_h] \cdot \{\tau - I_h \tau\} \mathbf{n}_e + \sum_{e \in \mathcal{E}_h^{\text{val}}} \int_{e} [v_h] \{\operatorname{div}(\tau - I_h \tau)\} \cdot \mathbf{n}_e \to 0 $$

and

$$ \int_{\Omega} H_h(v_h) : (\tau - I_h \tau) \to 0 $$

as $h \to 0^+$. This fact follows from the uniform boundedness assumptions $|v_h|_{H_h^2(\Omega)} \le C$ and standard interpolation estimates. We refer to [33, Lemma 2.4] for additional details when $\Sigma = \emptyset$ (see Appendix B for a proof in the free boundary case and Appendix C for the Dirichlet boundary conditions case). The case $\Sigma \ne \emptyset$ is derived similarly. Strong convergence properties are also available, for example, in [33].

We now turn our attention to the boundary conditions and start by noting that the Dirichlet/mixed boundary edges $e \in \mathcal{E}_h^D \cup \mathcal{E}_h^M$ are incorporated in the definition of the discrete Hessian. They do not influence the weak convergence property sketched above (because $\tau$ is compactly supported). However, the value of the deformation and its gradient imposed on $\Gamma^D \cup \Gamma^M$ must be
added separately for strong consistency [61]. To this end, we introduce the lifting of the boundary conditions (32) and (33)

$$L_h(\Phi_m, \varphi_m) := - \sum_{e \in \mathcal{E}_h^D} r_e(\Phi_m) + \sum_{e \in \mathcal{E}_h^D \cup \mathcal{E}_h^M} b_e(\varphi_m),$$

where we recall that $\Phi_m$ stands for the $m^{th}$ row of $\Phi$, $m = 1, 2, 3$. The discrete Hessian operator accounting for the boundary conditions is then chosen to be

$$H_h(\varphi_m, \Phi_m; y_{h,m}) := H_h(y_{h,m}) - L_h(\Phi_m, \varphi_m) \quad (38)$$

with the usual extension for vector-valued functions.

### 3.3. Discrete Energies

We are now in a position to describe the approximation of the different energies $E_s^{\text{pre}}(\mathbf{y})$ (preasymptotic), $E^{\text{pre}}(\mathbf{y})$ (prestrain), and $E^{\text{bil}}(\mathbf{y})$ (bilayer).

#### 3.3.1. Preasymptotic Energy

We start with the preasymptotic energy

$$E_s^{\text{pre}}(\mathbf{y}) = E^S(\mathbf{y}) + s^2 E^B(\mathbf{y})$$

and approximate the stretching and bending components separately. For the stretching part, we simply replace the derivatives appearing in (21) by their piecewiese counterparts and define for $\mathbf{y}_h \in [\mathbb{V}_h^k]^3$

$$E_h^S(\mathbf{y}_h) := \frac{\mu}{4} \int_{\Omega} |g^{-\frac{1}{2}} ((\nabla_h \mathbf{y}_h)^T \nabla_h \mathbf{y}_h - g) g^{-\frac{1}{2}}|^2 + \frac{\lambda}{8} \int_{\Omega} \text{tr} (g^{-\frac{1}{2}} ((\nabla_h \mathbf{y}_h)^T \nabla_h \mathbf{y}_h - g) g^{-\frac{1}{2}})^2. \quad (39)$$

For the bending part, we first replace as in [5] the second fundamental form in $E^B$ by the Hessian. This step is not justified but taken for computational convenience. We then take advantage of the discrete Hessian (38) and set for $\mathbf{y}_h \in [\mathbb{V}_h^k]^3$

$$E_h^B(\mathbf{y}_h) := \frac{\mu}{12} \sum_{m=1}^3 \int_{\Omega} |g^{-\frac{1}{2}} H_h(\varphi_m, \Phi_m; y_{h,m}) g^{-\frac{1}{2}}|^2 + \frac{\mu \lambda}{12(2\mu + \lambda)} \sum_{m=1}^3 \int_{\Omega} \text{tr} (g^{-\frac{1}{2}} H_h(\varphi_m, \Phi_m; y_{h,m}) g^{-\frac{1}{2}})^2; \quad (40)$$

compare with (22). In the preasymptotic case $\Sigma = \emptyset$ but for later use, we still define the approximation of the bending energy and subsequent quantities using $\Sigma$.

At this point neither $E_h^S$ nor $E_h^B$ enforce weak continuity condition nor are the boundary conditions (when $\Gamma^D \cup \Gamma^M \cup \Gamma^P \neq \emptyset$) imposed properly. This is the purpose of the stabilization term defined for $\mathbf{y}_h \in [\mathbb{V}_h^k]^3$ by

$$S_h(\mathbf{y}_h) := \frac{\gamma_0}{2} \left\| \mathbf{h}^{-\frac{3}{2}}[\mathbf{y}_h] \right\|_{L^2(\Gamma_h^0)}^2 + \frac{\gamma_0}{2} \left\| \mathbf{h}^{-\frac{3}{2}}(\mathbf{y}_h - \boldsymbol{\varphi}) \right\|_{L^2(\Gamma_h^{D \cup \Gamma_h^M})}^2 \\ + \frac{\gamma_1}{2} \left\| \mathbf{h}^{-\frac{1}{2}}[\nabla_h \mathbf{y}_h] \right\|_{L^2(\Gamma_h^{0 \setminus \Sigma_h})}^2 + \frac{\gamma_1}{2} \left\| \mathbf{h}^{-\frac{1}{2}}(\nabla_h \mathbf{y}_h - \boldsymbol{\Phi}) \right\|_{L^2(\Gamma_h^D)}^2 \\ + \frac{\gamma_2}{2} \sum_{\mathbf{x}_i \in \Gamma^P} h_i^{-2} |\mathbf{y}_h(\mathbf{x}_i) - \boldsymbol{\varphi}_i|^2, \quad (41)$$

where $\gamma_0, \gamma_1, \gamma_2 > 0$ are stabilization parameters. Note that a particularity of the LDG method is that these parameters are not required to be sufficiently large.
It remains to incorporate possible external forces - see (31) - to arrive at the following approximation of $E_s^{\text{pre}}$

$$E_{s,h}^{\text{pre}}(\mathbf{y}_h) := E_h^S(\mathbf{y}_h) + s^2 (E_h^B(\mathbf{y}_h) + S_h(\mathbf{y}_h) - F(\mathbf{y}_h)).$$

The discrete problem is thus to seek $\mathbf{y}_h^* \in [\mathbb{V}_h^k]^3$ such that

$$\begin{gathered} \mathbf{y}_h^* \in \arg\min_{\mathbf{y}_h \in [\mathbb{V}_h^k]^3} E_{s,h}^{\text{pre}}(\mathbf{y}_h). \quad (42) \\ \mathbf{y}_h \in [\mathbb{V}_h^k]^3 \end{gathered}$$

The $\Gamma$-convergence of $E_{s,h}^{\text{pre}}$ to $E_s^{\text{pre}}$ as $h \to 0^+$ is proved in [32] when using continuous (rather than discontinuous) piecewise polynomial approximations of the deformations and when $\Sigma = \emptyset$ (i.e. without folding).

### 3.3.2. Prestrain Energy

The approximation of the prestrain energy $E^{\text{pre}}$ in (25) directly follows from the above discussion since $E^{\text{pre}} = E^B$. Accounting for the potential external forces and boundary conditions, we thus set

$$E_h^{\text{pre}}(\mathbf{y}_h) := E_h^B(\mathbf{y}_h) + S_h(\mathbf{y}_h) - F(\mathbf{y}_h). \quad (43)$$

However, unlike the preasymptotic model, the deformations are required to satisfy the metric constraint (23). We have already pointed out that when using polynomial approximations, the metric cannot be satisfied everywhere in $\Omega$. We thus consider the following relaxation of the constraint

$$\mathcal{D}_h^{\text{aver}}(\mathbf{y}_h) := \sum_{T \in \mathcal{T}_h} \left| \int_T ((\nabla \mathbf{y}_h)^T \nabla \mathbf{y}_h - g) \right| \le \varepsilon \quad (44)$$

for some $\varepsilon > 0$ sufficiently small. This leads to the definition of the discrete admissible sets

$$\mathcal{A}_{h,\varepsilon}^{\text{aver}} := \{ \mathbf{y}_h \in [\mathbb{V}_h^k]^3 : \mathcal{D}_h^{\text{aver}}(\mathbf{y}_h) \le \varepsilon \}. \quad (45)$$

With these notations, the discrete prestrained minimization problem originally proposed in [27, 33] is to find $\mathbf{y}_h^* \in [\mathbb{V}_h^k]^3$ such that

$$\mathbf{y}_h^* \in \arg\min_{\mathbf{y}_h \in \mathcal{A}_{h,\varepsilon}^{\text{aver}}} E_h^{\text{pre}}(\mathbf{y}_h). \quad (46)$$

The $\Gamma$-convergence of $E_h^{\text{pre}}$ to $E^{\text{pre}}$ as $h \to 0^+$ is proved in [33] for the case $\Sigma = \emptyset$ (i.e. without folding). Moreover, the case $\Sigma \ne \emptyset$ is considered in [30] for the bending problem with isometry constraint (i.e. $g = I_2$) assuming that $\Sigma_h$ in (34) is an exact representation of the folding curve $\Sigma$. Finally, it is worth mentioning that extension to approximate curves obtained using isoparametric finite elements of degree $k \ge 2$ is investigated [36] for the linear bending problem, namely for the bending problem without the isometry constraint which is suitable for small displacements.

### 3.3.3. Bilayer Energy

The bilayer energy (27) is composed of two terms: the bending term $E^B$ (with $g = I_2$) approximated by $E_h^B$ and the nonlinear term

$$N(\mathbf{y}) := \alpha \sum_{i,j=1}^{2} \int_{\Omega} \partial_{ij} \mathbf{y} \cdot (\partial_1 \mathbf{y} \times \partial_2 \mathbf{y}) Z_{ij}.$$

Note that both terms are fourth order terms, requiring the nonlinearity to be discretized with care. To motivate the proposed scheme, we first note that

$$N(\mathbf{y}) \lesssim \|\mathbf{y}\|_{H^2(\Omega)} \|\nabla \mathbf{y}\|_{L^\infty(\Omega)}^2 \quad (47)$$
which is bounded for deformations $\mathbf{y} \in [H^2(\Omega)]^3$ satisfying the isometry constraint $\nabla \mathbf{y}^T \nabla \mathbf{y} = I_2$ (and thus $|\nabla \mathbf{y}| \in L^\infty(\Omega)$). In particular, the control on the nonlinear term is only possible thanks to the isometry constraint.

The relaxation (44) of the metric constraint does not provide enough control. Instead, we introduce an alternate discrete admissible set which guarantees a pointwise control at all the barycenter $x_T$ of the elements $T$ in the subdivision $\mathcal{T}_h$:

$$ \mathcal{D}_h^{\text{bary}}(\mathbf{y}_h) := \max_{T \in \mathcal{T}_h} |[(\nabla \mathbf{y}_h)^T \nabla \mathbf{y}_h - I_2](x_T)| \le \varepsilon \quad (48) $$

for $\varepsilon > 0$ sufficiently small; compare with (44). This leads to the definition of the discrete admissible set

$$ \mathcal{A}_{h,\varepsilon}^{\text{bary}} := \{ \mathbf{y}_h \in [\mathcal{V}_h^k]^3 : \mathcal{D}_h^{\text{bary}}(\mathbf{y}_h) \le \varepsilon \}. \quad (49) $$

The approximate constraint (48) implies (see [31])

$$ 1 - \varepsilon \le |\partial_i \mathbf{y}_h(x_T)| \le 1 + \varepsilon, \quad i = 1, 2, $$

so that $|\nabla_h \mathbf{y}_h(x_T)|$ is uniformly bounded for all $T \in \mathcal{T}_h$. Whence, in order to obtain a discrete version of (47), we approximate the nonlinear term $N(\mathbf{y}_h)$ using a 1-point quadrature with quadrature points localized at the barycenter of each element $T \in \mathcal{T}_h$

$$ N(\mathbf{y}) \approx \alpha \sum_{i,j=1}^{2} \sum_{T \in \mathcal{T}_h} |T| [\partial_{ij} \mathbf{y} \cdot (\partial_1 \mathbf{y} \times \partial_2 \mathbf{y}) Z_{ij}] (x_T). $$

In addition, it would be tempting to replace the second derivatives $\partial_{ij}\mathbf{y}(x_T)$ by the discrete hessian $H_h(\mathbf{y}_h)(x_T)$. However, estimates on the quadrature error would require regularity of $\mathbf{y}$ beyond $[H^2(\Omega)]^3$. To circumvent this issue, we introduce the $L^2(\Omega)$ projection $\overline{H}_h(\mathbf{y}_h)$ of $H_h(\mathbf{y}_h)$ onto the piecewise constant tensors and set

$$ N_h(\mathbf{y}_h) := \alpha \sum_{i,j=1}^{2} \sum_{T \in \mathcal{T}_h} |T| \left[ (\overline{H}_h(\mathbf{y}_h))_{ij} \cdot (\partial_1 \mathbf{y}_h \times \partial_2 \mathbf{y}_h) Z_{ij} \right] (x_T). $$

The discrete bilayer minimization problem is to seek $\mathbf{y}_h^* \in [\mathcal{V}_h^k]^3$ such that

$$ \mathbf{y}_h^* \in \underset{\mathbf{y}_h \in \mathcal{A}_{h,\varepsilon}^{\text{bary}}}{\operatorname{arg\,min}} E_h^{\text{bil}}(\mathbf{y}_h), \quad (50) $$

where for $\mathbf{y}_h \in [\mathcal{V}_h^k]^3$,

$$ E_h^{\text{bil}}(\mathbf{y}_h) := E_h^B(\mathbf{y}_h) + S_h(\mathbf{y}_h) - F(\mathbf{y}_h) - N_h(\mathbf{y}_h). \quad (51) $$

Note that the discrete minimization problem (50) was introduced in [31] following ideas from [29, 38]. Moreover, the proof that $E_h^{\text{bil}} \Gamma$-converges to $E^{\text{bil}}$ as $h \to 0^+$ can be found in [31].

## 4. Energy Minimization: Discrete Gradient Flows

Now that we have established the discrete minimization problems, it remains to discuss numerical procedures to construct approximate minimizers to each problem.

We recall that the considered energies are nonconvex and/or the minimization problem is subject to a nonconvex constraint. Hence, we use a gradient flow for its robustness and ability to approximately satisfy the type of nonconvex constraints encountered in our context. These benefits come at the expense of slow convergence towards minimizers.

Gradient flows require an initial condition, and their efficiency depends on both the energy of the initial condition and, for constrained problems, how well the initial condition satisfies the constraint. After introducing the (main) discrete gradient flow, we discuss preprocessing algorithms to construct suitable initial deformations. A bi-harmonic problem is advocated to enforce the boundary conditions (recall that the boundary conditions are part of the energies)
while an unconstrained gradient flow is put forward for reducing the metric defect of the initial deformation.

For the unconstrained preasymptotic minimization problem and in the metric preprocessing algorithm, an accelerated (Nesterov-type) algorithm considerably improves the number of iterations required to achieve a near minimum. This is the topic of the last part of this section.

## 4.1. Main Discrete Gradient Flows

Discrete gradient flows are used to minimize the discrete energies introduced in Section 3.3. Note that we have two types of minimization problems, namely an unconstrained minimization problem of the form

$$ \min_{\mathbf{y}_h \in [\mathbb{V}_h^k]^3} E_h(\mathbf{y}_h) $$

for the preasymptotic case and a constrained minimization problem of the form

$$ \min_{\mathbf{y}_h \in \mathbb{A}_{h, \varepsilon}^i} E_h(\mathbf{y}_h), \quad i \in \{\text{aver}, \text{bary}\}, $$

for all the other cases.

Gradient flows are chosen to minimize the different energies because of their robustness (energy decreasing property). In all the cases, we advocate an $H_h^2(\Omega)$ gradient flow based on the mesh-dependent inner product $(\cdot, \cdot)_{H_h^2(\Omega)}$ defined for any $\mathbf{v}_h, \mathbf{w}_h \in [\mathbb{V}_h^k]^3$ by

$$ (\mathbf{v}_h, \mathbf{w}_h)_{H_h^2(\Omega)} := \sigma(\mathbf{v}_h, \mathbf{w}_h)_{L^2(\Omega)} + \langle \mathbf{v}_h, \mathbf{w}_h \rangle_{H_h^2(\Omega)}, \quad (52) $$

where $\langle \mathbf{v}_h, \mathbf{w}_h \rangle_{H_h^2(\Omega)}$ is given by (36). Here $\sigma = 0$ if $\Gamma^D \neq \emptyset$ and $\sigma = 1$ otherwise to ensure that (52) is indeed an inner product on $[\mathbb{V}_h^k]^3$ when no Dirichlet boundary conditions are imposed. The gradient flow metric is thus

$$ \|\mathbf{v}_h\|_{H_h^2(\Omega)} := (\mathbf{v}_h, \mathbf{v}_h)_{H_h^2(\Omega)}^{1/2}, \quad \mathbf{v}_h \in [\mathbb{V}_h^k]^3. $$

We detail below the discrete gradient flows for the different energies. In order to simplify the discussion, we do not include the possible contribution from external forces (in other words we assume that $\mathbf{f} = 0$).

### 4.1.1. Preasymptotic Energy

We start with the unconstrained minimization problem (42) and recall that without external forces

$$ E_{s,h}^{\text{pre}}(\mathbf{y}_h) = E_h^S(\mathbf{y}_h) + s^2 (E_h^B(\mathbf{y}_h) + S_h(\mathbf{y}_h)). $$

A minimizing movements procedure is adopted to determine successive approximations reducing $E_{s,h}^{\text{pre}}$. The minimization process requires a pseudo time-step $\tau > 0$ and an initial guess $\mathbf{y}_h^0 \in [\mathbb{V}_h^k]^3$. Ideally, each iteration of the algorithm would read: given $\mathbf{y}_h^n \in [\mathbb{V}_h^k]^3$, find a deformation minimizing

$$ \mathbf{v}_h \mapsto \left( \frac{1}{2\tau} \|\mathbf{v}_h - \mathbf{y}_h^n\|_{H_h^2(\Omega)}^2 + E_{s,h}^{\text{pre}}(\mathbf{v}_h) \right) $$

over $[\mathbb{V}_h^k]^3$. Note that the Euler-Lagrange equation associated with the above minimization problem is nonlinear, and solving it would entail using an iterative method. Instead, we consider the following linearization

$$ \delta E_{s,h}^{\text{pre}}(\mathbf{y}_h^{n+1}, \mathbf{w}_h) \approx \delta E_h^S(\mathbf{y}_h^n, \mathbf{y}_h^{n+1}, \mathbf{w}_h) + s^2 (\delta E_h^B(\mathbf{y}_h^{n+1}, \mathbf{w}_h) + \delta S_h(\mathbf{y}_h^{n+1}, \mathbf{w}_h)). $$
Here $\delta E_h^B(\mathbf{v}_h, \mathbf{w}_h)$ is the Gâteau derivative of $E_h^B$ (see (40)) at $\mathbf{v}_h$ in the direction $\mathbf{w}_h$, $\delta E_h^S(\mathbf{y}_h^n; \mathbf{v}_h, \mathbf{w}_h)$ is the following linearization at $\mathbf{y}_h^n$ of the Gâteau derivative of $E_h^S$ (see (39)) at $\mathbf{v}_h$ in the direction $\mathbf{w}_h$

$$ \delta E_h^S (\mathbf{y}_h^n; \mathbf{v}_h, \mathbf{w}_h) := \frac{\mu}{2} \int_{\Omega} g^{-\frac{1}{2}} l(\mathbf{v}_h, \mathbf{w}_h) g^{-\frac{1}{2}} : g^{-\frac{1}{2}} m(\mathbf{y}_h^n) g^{-\frac{1}{2}} + \frac{\lambda}{4} \int_{\Omega} \mathrm{tr} (g^{-\frac{1}{2}} l(\mathbf{v}_h, \mathbf{w}_h) g^{-\frac{1}{2}}) \mathrm{tr} (g^{-\frac{1}{2}} m(\mathbf{y}_h^n) g^{-\frac{1}{2}}), $$

where

$$ l(\mathbf{v}_h, \mathbf{w}_h) := (\nabla_h \mathbf{v}_h)^T \nabla_h \mathbf{w}_h + (\nabla_h \mathbf{v}_h)^T \nabla_h \mathbf{w}_h \quad \text{and} \quad m(\mathbf{y}_h^n) := (\nabla_h \mathbf{y}_h^n)^T \nabla_h \mathbf{y}_h^n - g, $$

and $\delta S_h(\mathbf{v}_h, \mathbf{w}_h)$ is the Gâteau derivative of $S_h$ (see (41)) at $\mathbf{v}_h$ in the direction $\mathbf{w}_h$.

With these notations, we define $\mathbf{y}_h^{n+1} \in [\mathbb{V}_h^k]^3$ to be solution of

$$ \tau^{-1} (\mathbf{y}_h^{n+1} - \mathbf{y}_h^n, \mathbf{w}_h)_{H_h^2(\Omega)} + \delta E_h^S (\mathbf{y}_h^n, \mathbf{y}_h^{n+1}, \mathbf{w}_h) + s^2 (\delta E_h^B (\mathbf{y}_h^{n+1}, \mathbf{w}_h) + \delta S_h (\mathbf{y}_h^{n+1}, \mathbf{w}_h)) = 0 \quad (53) $$

for all $\mathbf{w}_h \in [\mathbb{V}_h^k]^3$. Note that (53) has a unique solution and the sequence of energies $\{E_{s,h}^{\mathrm{pre}}(\mathbf{y}_h^n)\}_{n \ge 0}$ is decreasing, i.e.

$$ \frac{1}{\tau} \| \mathbf{y}_h^{n+1} - \mathbf{y}_h^n \|_{H_h^2(\Omega)}^2 + E_{s,h}^{\mathrm{pre}}(\mathbf{y}_h^{n+1}) \le E_{s,h}^{\mathrm{pre}}(\mathbf{y}_h^n) $$

provided $\tau$ is sufficiently small (proportional to $h_{\min}$ and $E_{s,h}^{\mathrm{pre}}(\mathbf{y}_h^0)^{-1}$). As a consequence, we obtain $E_{s,h}^{\mathrm{pre}}(\mathbf{y}_h^{n+1}) < E_{s,h}^{\mathrm{pre}}(\mathbf{y}_h^n)$ if $\mathbf{y}_h^{n+1} \ne \mathbf{y}_h^n$. We refer to [32] for more details, see also [33].

### 4.1.2. Prestrain Energy

The prestrain energy $E_h^{\mathrm{pre}} = E_h^B + S_h$ appears simpler to reduce because it does not include the problematic non-quadratic stretching energy $E_h^S$. However, the metric condition appears as a constraint on the deformations, which have to belong to the admissible set $\mathcal{A}_{h,\varepsilon}^{\mathrm{aver}}$ (see (45)) for a prescribed $\varepsilon > 0$.

Note that the metric constraint $\nabla \mathbf{y}^T \nabla \mathbf{y} = g$ is relaxed in $\mathcal{A}_{h,\varepsilon}^{\mathrm{aver}}$ and only needs to be satisfied approximately. To achieve this, we compute increments in a (pseudo-)tangent space of the constraint and rely on the $H_h^2(\Omega)$ metric of the gradient flow to guarantee that all the deformations generated by the gradient flow indeed belong to $\mathcal{A}_{h,\varepsilon}^{\mathrm{aver}}$ for a specific choice of $\varepsilon$.

Given a deformation $\mathbf{y}_h^n \in [\mathbb{V}_h^k]^3$, we define the (pseudo-)tangent space at $\mathbf{y}_h^n$ as follows:

$$ \mathcal{W}_h^{\mathrm{aver}} (\mathbf{y}_h^n) := \left\{ \mathbf{w}_h \in [\mathbb{V}_h^k]^3 : \int_T \left( (\nabla_h \mathbf{w}_h)^T \nabla_h \mathbf{y}_h^n + (\nabla_h \mathbf{y}_h^n)^T \nabla_h \mathbf{w}_h \right) = 0 \quad \forall T \in \mathcal{T}_h \right\}. $$

Now, given a pseudo time-step $\tau > 0$ and an initial guess $\mathbf{y}_h^0 \in [\mathbb{V}_h^k]^3$, the discrete gradient flow algorithm for the constrained minimization problem (46) consists in computing successively $\mathbf{y}_h^{n+1} \in [\mathbb{V}_h^k]^3$ such that $\mathbf{y}_h^{n+1} - \mathbf{y}_h^n \in \mathcal{W}_h^{\mathrm{aver}}(\mathbf{y}_h^n)$ and

$$ \tau^{-1} (\mathbf{y}_h^{n+1} - \mathbf{y}_h^n, \mathbf{w}_h)_{H_h^2(\Omega)} + \delta E_h^{\mathrm{pre}} (\mathbf{y}_h^{n+1}, \mathbf{w}_h) = 0 \quad \forall \mathbf{w}_h \in \mathcal{W}_h^{\mathrm{aver}} (\mathbf{y}_h^n), $$

where

$$ \delta E_h^{\mathrm{pre}} (\mathbf{y}_h^{n+1}, \mathbf{w}_h) := \delta E_h^B (\mathbf{y}_h^{n+1}, \mathbf{w}_h) + \delta S_h (\mathbf{y}_h^{n+1}, \mathbf{w}_h). $$

In practice, the linear constraint encoded in $\mathcal{W}_h^{\mathrm{aver}}(\mathbf{y}_h^n)$ is enforced using piecewise constant Lagrange multipliers $\Lambda_h : \Omega \to [\mathbb{V}_h^0]^{2 \times 2}$ with $\Lambda_h^T = \Lambda_h$.

The resulting sequence of deformations is again energy decreasing, i.e.

$$ \frac{1}{\tau} \| \mathbf{y}_h^{n+1} - \mathbf{y}_h^n \|_{H_h^2(\Omega)}^2 + E_h^{\mathrm{pre}} (\mathbf{y}_h^{n+1}) \le E_h^{\mathrm{pre}} (\mathbf{y}_h^n) \quad (54) $$

irrespective of the choice of the pseudo timestep $\tau > 0$. Furthermore, because $\mathbf{y}_h^{n+1} - \mathbf{y}_h^n \in \mathcal{W}_h^{\mathrm{aver}}(\mathbf{y}_h^n)$, we can show that the metric defect $\mathcal{D}_h^{\mathrm{aver}}$ defined in (44) is uniformly controlled, namely

$$ \mathcal{D}_h^{\mathrm{aver}} (\mathbf{y}_h^n) \le \mathcal{D}_h^{\mathrm{aver}} (\mathbf{y}_h^0) + c_1 \tau (E_h^{\mathrm{pre}} (\mathbf{y}_h^0) + c_2), \quad n = 1, 2, \dots, \quad (55) $$
where the constant $c_1 > 0$ depends only on the constants in Poincaré-Friedrichs-type inequalities and $c_2 \ge 0$ depends on the input data (namely $\mu$, $g$, $\varphi$, and $\Phi$) as well as the stabilization parameters $\gamma_0, \gamma_1$. We refer to [33] for more details when $\Gamma^P = \emptyset$.

The control on the metric defect (55) implies that upon setting $\varepsilon = \mathcal{D}_h^{\text{aver}}(\mathbf{y}_h^0) + c_1\tau(E_h^{\text{pre}}(\mathbf{y}_h^0) + c_2)$, $\mathbf{y}_h^n \in \mathcal{A}_{h,\varepsilon}^{\text{aver}}$ for all $n = 1, 2, \dots$. In particular, the recursive algorithm produces a sequence of deformation in $\mathcal{A}_{h,\varepsilon}^{\text{aver}}$ with decreasing energy as desired.

### 4.1.3. Bilayer Energy

It remains to discuss the constrained minimization problem (50) associated with the bilayer energy

$$E_h^{\text{bil}} = E_h^B + S_h - N_h,$$

see (51) (recall that $\mathbf{f} = \mathbf{0}$). This case contains the two difficulties encountered above: the Gâteau derivative $\delta E_h^{\text{bil}}$ is not bilinear due to the cubic term $N_h$ and the minimization is subject to an approximate isometry constraint. Note that the latter is enforced at the cells' barycenters rather than in average, see (49).

Nevertheless, the algorithm considered in [31], which in turn is inspired from [38], combines the ideas of the previous two cases. We define for $\mathbf{y}_h^n \in [\mathbb{V}_h^k]^3$ the set

$$\mathbb{W}_h^{\text{bary}}(\mathbf{y}_h^n) := \left\{ \mathbf{w}_h \in [\mathbb{V}_h^k]^3 : \left[ (\nabla_h \mathbf{w}_h)^T \nabla_h \mathbf{y}_h^n + (\nabla_h \mathbf{y}_h^n)^T \nabla_h \mathbf{w}_h \right] (x_T) = 0 \quad \forall T \in \mathcal{T}_h \right\}$$

and seek $\mathbf{y}_h^{n+1} \in [\mathbb{V}_h^k]^3$ such that $\mathbf{y}_h^{n+1} - \mathbf{y}_h^n \in \mathbb{W}_h^{\text{bary}}(\mathbf{y}_h^n)$ and

$$\tau^{-1}(\mathbf{y}_h^{n+1} - \mathbf{y}_h^n, \mathbf{w}_h)_{H_h^2(\Omega)} + \delta E_h^B(\mathbf{y}_h^{n+1}, \mathbf{w}_h) + \delta S_h(\mathbf{y}_h^{n+1}, \mathbf{w}_h) = N_h(\mathbf{y}_h^n, \mathbf{w}_h) \quad (56)$$

for all $\mathbf{w}_h \in \mathbb{W}_h^{\text{bary}}(\mathbf{y}_h^n)$. Here

$$\begin{aligned} N_h(\mathbf{y}_h^n; \mathbf{w}_h) := & \alpha \sum_{i,j=1}^{2} \sum_{T \in \mathcal{T}_h} |T| \left[ (\overline{H}_h(\mathbf{w}_h))_{ij} \cdot (\partial_1 \mathbf{y}_h^n \times \partial_2 \mathbf{y}_h^n) Z_{ij} \right] (x_T) \\ & + \alpha \sum_{i,j=1}^{2} \sum_{T \in \mathcal{T}_h} |T| \left[ (\overline{H}_h(\mathbf{y}_h^n))_{ij} \cdot (\partial_1 \mathbf{w}_h \times \partial_2 \mathbf{y}_h^n) Z_{ij} \right] (x_T) \\ & + \alpha \sum_{i,j=1}^{2} \sum_{T \in \mathcal{T}_h} |T| \left[ (\overline{H}_h(\mathbf{y}_h^n))_{ij} \cdot (\partial_1 \mathbf{y}_h^n \times \partial_2 \mathbf{w}_h) Z_{ij} \right] (x_T) \end{aligned}$$

so that (56) corresponds to an explicit treatment of the cubic term $N_h$. As in the prestrain case, the linear constraint encoded in $\mathbb{W}_h^{\text{aver}}(\mathbf{y}_h^n)$ is enforced using piecewise constant Lagrange multipliers $\Lambda_h : \Omega \to [\mathbb{V}_h^0]^{2 \times 2}$ with $\Lambda_h^T = \Lambda_h$.

As it turns out, except for the mild requirement $\tau \lesssim |\ln(h_{\min})|^{-1}$, the explicit treatment of $N_h$ does not affect the convergence property (in [29] its treatment was implicit thereby requiring a fixed-point iteration at each step). In particular, the energy decay property (54) holds for $E_h^{\text{bil}}$ with a factor $\frac{1}{2\tau}$ instead of $\frac{1}{\tau}$ in front of the metric term and the isometry defect satisfies

$$\mathcal{D}_h^{\text{bary}}(\mathbf{y}_h^n) \le \mathcal{D}_h^{\text{bary}}(\mathbf{y}_h^0) + c_1 \tau |\log(h_{\min})| (E_h(\mathbf{y}_h^0) + c_2), \quad n = 1, 2, \dots, \quad (57)$$

where similarly to (55), the constant $c_1 > 0$ depends only on the constants in Poincaré-Friedrichs-type inequalities while $c_2 \ge 0$ depends on the input data (in this case $Z$, $\varphi$, and $\Phi$) as well as the stabilization parameters $\gamma_0, \gamma_1$. Refer to [31] for additional details in the case $\Gamma^P = \emptyset$. Again, this implies that all the iterations of the algorithms belong to $\mathcal{A}_{\varepsilon,h}^{\text{bary}}$ provided that $\varepsilon \ge \mathcal{D}_h^{\text{bary}}(\mathbf{y}_h^0) + c_1\tau|\log(h_{\min})|(E_h(\mathbf{y}_h^0) + c_2)$.
## 4.2. Preprocessing Step

In this section we introduce a preprocessing strategy which can be used to construct the initial deformation $\mathbf{y}_h^0$ required by the gradient flows detailed in Section 4. The estimates (55) and (57) for the metric defect satisfied by the constraint minimization problem suggest that the initial condition should satisfy approximately the boundary conditions (to reduce the energy) and have a small metric defect. Therefore, the preprocessing strategy consists of two main modules dedicated to each aspect. Note that in some cases, a suitable initial deformation is known a priori, as in the case of clamped boundary conditions and isometry constraint ($g = I_2$), which are satisfied by the identity map $\mathbf{y}_h^0(x_1, x_2) = (x_1, x_2, 0)^T$. However, in general, the initial deformation is not accessible.

### 4.2.1. Boundary Conditions Preprocessing

The BC preprocessing consists in solving a bi-Lapacian problem to obtain a deformation satisfying (approximately) the prescribed boundary conditions. We seek an approximation of the solution to the bi-Laplacian problem

$$ egin{cases} \Delta^2 \hat{\mathbf{y}} = \hat{\mathbf{f}} & \text{in } \Omega \\ \nabla \hat{\mathbf{y}} = \Phi & \text{on } \Gamma^D \\ \hat{\mathbf{y}} = \boldsymbol{\varphi} & \text{on } \Gamma^D \cup \Gamma^M \\ \hat{\mathbf{y}}(\mathbf{x}_i) = \boldsymbol{\varphi}_i & \text{for all } \mathbf{x}_i \in \Gamma^P, \end{cases} \qquad (58) $$

supplemented with the following natural boundary conditions

$$ \left\{ \begin{array}{ll} (D^2 \hat{\mathbf{y}}_m) \mathbf{n} = \mathbf{0} & \text{on } \partial\Omega \setminus \Gamma^D \\ \nabla(\Delta \hat{\mathbf{y}}_m) \cdot \mathbf{n} = 0 & \text{on } \partial\Omega \setminus (\Gamma^D \cup \Gamma^M) \end{array} \right. \qquad (59) $$

for $m = 1, 2, 3$. Note that $\hat{\mathbf{f}}$ is a fictitious forcing term. Its value is irrelevant but can be employed to obtain a non-planar deformation, namely a deformation $\hat{\mathbf{y}} = (y_1, y_2, y_3)$ with $y_3 \neq 0$, in particular when $\Gamma^D = \emptyset$. This happens to be critical because the main discrete gradient flow preserves planar configurations (when $\mathbf{f} = \mathbf{0}$) and can therefore not reach non-planar minimizers.

The approximation $\hat{\mathbf{y}}_h$ of $\hat{\mathbf{y}}$ is obtained using an LDG approach, that is

$$ \hat{\mathbf{y}}_h := \underset{\mathbf{y}_h \in [\mathbb{V}_h^k]^3}{\arg\min} E_h^{\text{pre}}(\mathbf{y}_h), $$

where $E_h^{\text{pre}}$ is given by (43) with $\mu = 6$, $\lambda = 0$, and $g = I_2$. In particular, the solution $\hat{\mathbf{y}}_h$ is obtained upon solving the associated Euler-Lagrange equation.

### 4.2.2. Metric Constraint Preprocessing

The metric preprocessing aims at minimizing the metric defect $\mathcal{D}_h^i$, $i \in \{\text{aver}, \text{bary}\}$. It consists of minimizing the simplified stretching energy

$$ E_h^{\text{str}}(\mathbf{y}_h) := \frac{1}{2} \int_{\Omega} |(\nabla_h \mathbf{y}_h)^T \nabla_h \mathbf{y}_h - g|^2 \qquad (60) $$

over $[\mathbb{V}_h^k]^3$. The stretching energy $E_h^{\text{str}}$ is not quadratic, but its reduction can be performed using an $H_h^2(\Omega)$ discrete gradient flow coupled with a linearization similar to the one for $E_{s,h}^{\text{pre}}$. Setting $\tilde{\mathbf{y}}_h^0 = \hat{\mathbf{y}}_h$, the deformation obtained by the BC preprocessing algorithm, recursively computes $\tilde{\mathbf{y}}_h^{n+1} \in [\mathbb{V}_h^k]^3$ as the solution to

$$ \tilde{\tau}^{-1}(\tilde{\mathbf{y}}_h^{n+1} - \tilde{\mathbf{y}}_h^n, \mathbf{w}_h)_{H_h^2(\Omega)} + \delta E_h^{\text{str}}(\tilde{\mathbf{y}}_h^n; \tilde{\mathbf{y}}_h^{n+1}, \mathbf{w}_h) = 0 \quad \forall \mathbf{w}_h \in [\mathbb{V}_h^k]^3, \qquad (61) $$
where $\tilde{\tau} > 0$ is a given (sufficiently small) pseudo time-step and where for $\mathbf{v}_h, \mathbf{w}_h \in [\mathbb{V}_h^k]^3$

$$ \delta E_h^{\text{str}}(\tilde{\mathbf{y}}_h^n; \mathbf{v}_h, \mathbf{w}_h) := \int_{\Omega} ((\nabla_h \mathbf{v}_h)^T \nabla_h \mathbf{w}_h + (\nabla_h \mathbf{w}_h)^T \nabla_h \mathbf{v}_h) : ((\nabla_h \tilde{\mathbf{y}}_h^n)^T \nabla_h \tilde{\mathbf{y}}_h^n - \mathbf{g}). $$

The gradient flow produces iterates with decreasing energy provided that $\tilde{\tau}$ is sufficiently small compared to $E_h^{\text{str}}(\tilde{\mathbf{y}}_h^0)^{-1}$ and $h_{\min}$, see [33]. Similar to the main discrete gradient flow, one draw-back of the metric constraint preprocessing is that planar configurations are local minimizers of (60) regardless of the target metric $g$, see [27]. A non-planar deformation can nevertheless be obtained using a non-planar initial deformation $\tilde{\mathbf{y}}_h^0$. Such an initial deformation can be generated, for instance, by solving (58)-(59) with $\hat{\mathbf{f}} \neq \mathbf{0}$. Moreover, the metric constraint preprocessing does not necessarily construct a deformation with finite energy $E_h^{\text{pre}}$. If needed, a bending term could be added to the energy (60) to generate iterates with a uniformly bounded $H_h^2(\Omega)$ semi-norm, see [33].

### 4.3. Accelerated Algorithm

It is well-documented that gradient flows exhibit slow convergence towards minimizers. Several accelerated algorithms have been introduced to improve the convergence rate of the standard gradient descent algorithm, but they are mainly available for convex problems.

Below we design an accelerated algorithm suited for the minimization of the preasymptotic energy and the stretching energy used in the preprocessing step. The proposed algorithm is inspired by Nesterov's original accelerated gradient method [62], see also [63, 64].

Typical results for the minimization of functions $f: \mathbb{R}^d \to \mathbb{R}$ is an $\mathcal{O}(n^{-2})$ convergence in $f$ when $f$ is a $C^1$ convex function with Lipschitz gradient. Here $n$ denotes the number of iterations. This is a significant improvement over the standard gradient method for which the convergence rate is only $\mathcal{O}(n^{-1})$.

The Fast Iterative Shrinkage/Thresholding Algorithm (FISTA) is an extension to when $f$ is the sum of two convex and lower-semicontinuous functions, one of class $C^1$ with Lipschitz gradient and the other potentially non-smooth. The convergence of the FISTA algorithm is also quadratic in the number of iterations, see [65] and [66].

Inspired by the FISTA algorithm (case $f_1 = 0$), but adapted to the functional setting, our accelerated algorithm for the minimization of the preasymptotic energy $E_{s,h}^{\text{pre}}$ and the stretching energy $E_h^{\text{str}}$ used in the metric constraint preprocessing reads. Given $\mathbf{y}_h^0 \in [\mathbb{V}_h^k]^3$, set $\mathbf{v}_h^0 = \mathbf{y}_h^0$. Then for $n = 0, 1, \dots$ do

(a) Find $\mathbf{y}_h^{n+1} \in [\mathbb{V}_h^k]^3$ satisfying (53) or (61) with $\mathbf{y}_h^n$ replaced by $\mathbf{v}_h^n$;

(b) Set $\mathbf{v}_h^{n+1} = \mathbf{y}_h^{n+1} + \eta_{n+1}(\mathbf{y}_h^{n+1} - \mathbf{y}_h^n)$, where $\eta_{n+1} = \frac{t_{n+1}-1}{t_{n+2}}$ and $\{t_n\}_{n \ge 1}$ satisfies

$$ t_1 = 1 \quad \text{and} \quad t_{n+1} = \sqrt{t_n^2 + \frac{1}{4} + \frac{1}{2}} \quad \text{for } n = 0, 1, \dots $$

As in the case of the main gradient flow, a linearization is performed at $\mathbf{v}_h^n$ to avoid a nonlinear system at each iteration.

Note that the accelerating algorithm does not have the energy decreasing property that the gradient descent algorithm does. This can be observed experimentally, see Figure 4 below. However, the algorithm does appear to converge when the appropriate step size is chosen.

Finding analytic estimates for the convergence rate of the proposed algorithm remains an open problem. In Section 5.1.2 we observe at least an $\mathcal{O}(n^{-2})$ convergence rate and illustrate the significant gain when the accelerated version of the algorithm is used.
### 4.4. Dynamics: Time Discretization

Several numerical experiments presented in Section 5 are dynamical in the sense that the data (boundary conditions and external forces) vary over the time interval $[0, T]$. The minimizers of the different energies $y$ are thus time dependent as well.

For a positive integer $M$, we let $\Delta t := T/M$ be the physical time-step (as opposed to the pseudo-time step used in the gradient flows) and consider the points $t_m := m\Delta t$ for $m = 0, 1, \dots, M$. We assume that the elastic relaxation is faster than the time relaxation, i.e. at each time step, the thin elastic structures have minimal energies. In particular, the deformations $y_h(t_m) \in [\mathbb{V}_h^k]_3$ at each time $t_m$ are discrete minimizers of the discrete energy defined using the input data evaluated at $t = t_m$.

The dynamic algorithm reads:

**Initialization** ($t = 0$): Obtain $\tilde{y}_h^{(0)}$ using the BC preprocessing algorithm with the prescribed conditions evaluated at $t = 0$; Set $y_h^{(0)} = \tilde{y}_h^{(0)}$;

**Dynamics** ($t \in (0, T]$), for $m = 1, 2, \dots, M$, do

(i) Obtain $\delta\tilde{y}_h^{(m)} := \tilde{y}_h^{(m)} - y_h^{(m-1)}$ using the BC preprocessing algorithm with the increment boundary conditions

$$\nabla \delta\tilde{y}_h^{(m)} = \Phi(t_m) - \Phi(t_{m-1}), \quad \text{on } \Gamma^D$$

$$\delta\tilde{y}_h^{(m)} = \varphi(t_m) - \varphi(t_{m-1}), \quad \text{on } \Gamma^D \cup \Gamma^M$$

$$\delta\tilde{y}_h^{(m)}(\mathbf{x}_i) = \varphi_i(t_m) - \varphi_i(t_{m-1}), \quad \text{for } i = 1, 2, \dots, P;$$

(ii) Starting from $\tilde{y}_h^{(m)}$, obtain $\tilde{y}_h^{(m)}$ using the metric preprocessing algorithm to satisfy approximately the metric constraint;

(iii) Starting from $\tilde{y}_h^{(m)}$, obtain $y_h^{(m)}$ using the (main) gradient flow to minimize the discrete energy with the data evaluated at $t = t_m$.

It is important to point out that the BC preprocessing procedure embedded in the dynamic algorithm above determines the increment $\delta\tilde{y}_h^{(m)} := \tilde{y}_h^{(m)} - y_h^{(m-1)}$ (rather than the preprocessed deformation $\tilde{y}_h^{(m)}$ directly) by solving (58)-(59) with incremental boundary conditions. This is done to take advantage of the previous step and, in particular, to avoid using a costly metric preprocessing algorithm starting from scratch at each time step.

## 5. Numerical Experiments

In this section, we illustrate the performance of the different algorithms. All the simulations are performed with $\mu = 6$ and $\lambda = 8$ (note that $\mu = 6$ and $\lambda = 0$ yield the factor $1/2$ for the bilayer case) and without any external force (i.e. $\mathbf{f} = \mathbf{0}$) except for the experiment in Section 5.3. The polynomial degree used for the approximated deformations $y_h$ is chosen to be $k = 2$, and, unless specified otherwise, we take $\gamma_0 = \gamma_1 = 1$ and $\gamma_2 = 10$ for the stabilization parameters. The main gradient flow ends when

$$\tau^{-1} |E_h(\mathbf{y}_h^{n+1}) - E_h(\mathbf{y}_h^n)| \le tol$$

for a prescribed tolerance $tol$ specified for each numerical experiment below. For the metric preprocessing we stop when

$$\tilde{\tau}^{-1} |\tilde{E}_h(\tilde{\mathbf{y}}_h^{n+1}) - \tilde{E}_h(\tilde{\mathbf{y}}_h^n)| \le tol \quad \text{or} \quad \mathcal{D}_h^i(\tilde{\mathbf{y}}_h^{n+1}) \le \tilde{\varepsilon}_0, \quad i \in \{\text{aver, bary}\},$$

for some tolerances $tol$ and $\tilde{\varepsilon}_0$ specified below. Finally, for the experiments in Sections 5.1 and 5.3, where the domain $\Omega$ consists of the unit disc, we use a quadratic mapping $F_T$ in (35) to better approximate the domain.
## 5.1. Preasymptotic

We first consider two experiments using the preasymptotic model. The first experiment illustrates the effect of the thickness of the plate on the final configuration, while the second demonstrates the advantages of the Nesterov-type acceleration discussed in Section 4.3. In both experiments the computational domain $\Omega$ is a disc of radius one, and we use $\tau = 0.01$ for the gradient flow pseudo time-step.

### 5.1.1. Disc with Oscillating Boundary

Our first experiment is inspired by [6, 8], in which a hydrogel disc of negative Gaussian curvature is observed to develop more oscillations along the boundary as its thickness is reduced. The material is prestrained according to the metric

$$g = J^T (\tilde{g} \circ \zeta) J,$$

where $J$ is the Jacobian matrix for the change of variables $(r, \theta) = \zeta(x_1, x_2)$ from polar to Cartesian coordinates and $\tilde{g}(r, \theta)$ is the first fundamental form of the following deformation

$$\tilde{y}(r, \theta) = (r \cos(\theta), r \sin(\theta), 0.2r^4 \sin(6\theta)). \qquad (62)$$

This deformation corresponds to a disk with six wrinkles.

For the discretization, our subdivision consists of 1280 quadrilaterals and a total of 34560 = 1280·9·3 degrees of freedom. The initial deformation $y_h^0$ is taken to be the continuous Lagrange interpolant of (62). We set $tol = 10^{-8}$ for the stopping criterion and run the experiment for $s^2 = 10^{-1}, 10^{-2}, 10^{-3}$, and 0.

The computed deformations minimizing the preasymptotic energy $E_s^{pre}$ are provided in Figure 3. In accordance with the laboratory experiments [6, 8], the number of wrinkles increases from 0 to 6 as $s$ decreases.

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Final configurations for the oscillating boundary experiment. Left to right and top to bottom: $s^2 = 10^{-1}, 10^{-2}, 10^{-3}$, and 0.</figcaption></figure>
### 5.1.2. Accelerated Algorithm

This example illustrates the benefits of using the accelerated algorithm introduced in Section 4.3 over our original gradient flow. This time the material is characterized by a “bubble” metric with positive Gaussian curvature, namely

$$g(x_1, x_2) = \begin{pmatrix} 1 + \alpha \frac{\pi^2}{4} \cos\left(\frac{\pi}{2}(1-r)\right)^2 \frac{x_1^2}{r^2} & \alpha \frac{\pi^2}{4} \cos\left(\frac{\pi}{2}(1-r)\right)^2 \frac{x_1 x_2}{r^2} \\ \alpha \frac{\pi^2}{4} \cos\left(\frac{\pi}{2}(1-r)\right)^2 \frac{x_1 x_2}{r^2} & 1 + \alpha \frac{\pi^2}{4} \cos\left(\frac{\pi}{2}(1-r)\right)^2 \frac{x_2^2}{r^2} \end{pmatrix}$$

with $r = \sqrt{x_1^2 + x_2^2}$ and $\alpha = 0.2$.

The mesh for this experiment consists of 320 quadrilaterals and a total of 8640 = 320 · 9 · 3 degrees of freedom. The initial deformation is taken to be the continuous Lagrange interpolant of the shallow paraboloid $-\frac{x_1^2}{16} + \frac{x_2^2}{16}$. This gives a slight initial bending that allows the simulation to find a non-flat minimizer. (Recall from Section 4.2 that the flat configuration is a local minimizer of the preasymptotic energy.) For the thickness of the material, we use the values $s^2 = 10^{-3}$, $10^{-4}$, $10^{-5}$, and 0 and we set $tol = 10^{-6}$ as stopping criterion.

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Evolution of $E_{s,h}^{\text{pre}}$ during the minimization when using the standard gradient flow and the accelerated algorithm. The energy are reported for the first few iterations to illustrate better the convergence of the acceleration.</figcaption></figure>

The discrete stretching energy $E_h^S$ and bending energy $s^2(E_h^B + S_h)$ along with the number of iterations needed to reach the tolerance $tol$ are reported in Table 1 when using the standard gradient flow and the accelerated algorithm. We observe that the accelerated algorithm finishes in much fewer iterations for all values of the thickness $s$. The approximate preasymptotic energy $E_{s,h}^{\text{pre}} = E_h^S + s^2(E_h^B + S_h)$ versus the algorithm iteration number $n$ is depicted in Figure 4 for the different values of $s$ considered. When using a standard gradient flow, we observe numerically
that the energy decay is like $\mathcal{O}(n^{-1})$, while the accelerated algorithm exhibits a decay similar to $\mathcal{O}(n^{-2})$. The latter is the expected convergence rate when using Nesterov-type accelerations.

<table><caption>**Table 1.** Final energy and number of iterations for different values of $s^2$, with and without acceleration.</caption><thead><tr><th rowspan="2">s<sup>2</sup></th><th colspan="3">Discrete Gradient Flow</th><th colspan="3">Accelerated Gradient Flow</th></tr><tr><th>E<sub>h</sub><sup>S</sup></th><th>E<sub>h</sub><sup>B</sup></th><th>Iterations</th><th>E<sub>h</sub><sup>S</sup></th><th>E<sub>h</sub><sup>B</sup></th><th>Iterations</th></tr></thead><tbody><tr><td>0</td><td>2.390e-4</td><td>0</td><td>17746</td><td>8.098e-6</td><td>0</td><td>1007</td></tr><tr><td>10<sup>-5</sup></td><td>1.525e-4</td><td>2.685e-4</td><td>26784</td><td>2.694e-4</td><td>3.990e-4</td><td>359</td></tr><tr><td>10<sup>-4</sup></td><td>4.479e-5</td><td>5.951e-4</td><td>37799</td><td>6.350e-6</td><td>4.095e-4</td><td>896</td></tr><tr><td>10<sup>-3</sup></td><td>5.595e-5</td><td>2.843e-3</td><td>39533</td><td>3.920e-5</td><td>2.741e-3</td><td>664</td></tr></tbody></table>

## 5.2. Bilayer and Folding

We consider here two experiments taken from [67], referred to as *diamond* and *bird* in the following. They illustrate the rigidity, robustness, and great variety of shapes achievable with this technology. In these examples, the computational domain consists of a collection of subdomains delimited by creases where folding is allowed. The bilayer material is designed to yield a piecewise constant spontaneous curvature tensor Z on each subdomain with $\alpha = 1$.

In both experiments we set $\tau = 0.1$ for the gradient flow pseudo time-step and $tol = 10^{-3}$ for the stopping criterion.

### 5.2.1. Diamond

The *diamond* computational domain is provided in Figure 5 (left). It consists of the square $(-1.5, 1.5) \times (-1.5, 1.5)$ rotated counter-clockwise by an angle of $\pi/4$. The two dashed red curves represent the creases. They are quadratic Bézier curves with the origin as control point. The point $\mathbf{z}_1$ is at distance $1/3$ from $\mathbf{x}_1$ on the segment from $\mathbf{x}_1$ to $\mathbf{x}_2$, and similarly for the three other points. We denote by $Z_i \in \mathbb{R}^{2 \times 2}$ the spontaneous curvature associated with each subdomain $i \in \{1, 2, 3\}$ and set

$$Z_1 = Z_3 = 0.6I_2 \quad \text{and} \quad Z_2 = -0.6I_2.$$

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Computational domain (left) and subdivision (right) for the *diamond* setting. The dashed red curves are creases across which folding is possible.</figcaption></figure>
The subdivision $\mathcal{T}_h$ is depicted in Figure 5 (right). It consists of 449 quadrilaterals with maximal mesh size 0.331678 (total of 13470 degrees of freedom, $12123 = 449 \cdot 9 \cdot 3$ for the deformation $\mathbf{y}_h$ and $1347 = 449 \cdot 3$ for the Lagrange multipliers used to enforce the linearized isometry constraint). Moreover, the initial deformation $\mathbf{y}_h^0$ is taken to be the identity map $\mathbf{y}_h^0(\Omega) = \Omega \times \{0\}$ for which $E_h^{\text{bil}} = 3.24$ and $\mathcal{D}_h^{\text{bary}} = 8.72497 \cdot 10^{-15}$. The deformations obtained at several steps of the gradient flow are provided in Figure 6, including the final (equilibrium) deformation reached in 305 iterations. For the final deformation, the energy and isometry defect are $E_h^{\text{bil}} = 0.679318$ and $\mathcal{D}_h^{\text{bary}} = 0.0334151$, respectively.

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Deformations for the diamond experiment. Left to right and top to bottom: deformation obtained after 20, 70, 120, 305 (top view), and 305 (side view) iterations.</figcaption></figure>

### 5.2.2. Bird

The bird geometry is depicted in Figure 7. It consists of 26 curves: 17 for the boundary of $\Omega$ (black plain curves), 7 creases (red dashed curves), and 1 extra curve on the tail (black dotted curve) which is only used for the construction of the subdivision (i.e. no folding is possible across this curve). All the curves are quadratic or cubic Bézier curves, i.e. obtained using one or two control points. The latter are chosen so that we obtain a geometry similar to [67]. In this example, the piecewise constant spontaneous curvature tensors $Z_1, \dots, Z_9$ are given on each subdomain by

$$Z_1 = Z_5 = \begin{pmatrix} 0 & 0 \\ 0 & 0.4 \end{pmatrix}, \quad Z_2 = Z_4 = \begin{pmatrix} 0 & 0 \\ 0 & -0.3 \end{pmatrix}, \quad Z_3 = 0.7 I_2,$$

$$Z_6 = \begin{pmatrix} -0.2 & 0 \\ 0 & 0 \end{pmatrix}, \quad Z_7 = \begin{pmatrix} 0.7 & 0 \\ 0 & 0 \end{pmatrix} \quad \text{and} \quad Z_8 = Z_9 = -0.7 I_2.$$
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Computational domain for the bird experiment.</figcaption></figure>

The mesh is constituted of 1468 quadrilaterals with maximal mesh size 0.690837 (total of 44040 degrees of freedom, $39636 = 1468 \cdot 9 \cdot 3$ for the deformation $\mathbf{y}_h$ and $4404 = 1468 \cdot 3$ for the Lagrange multipliers used to enforce the linearized isometry constraint), see Figure 8.

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Mesh for the bird example.</figcaption></figure>

The initial deformation $\mathbf{y}_h^0$ is again taken to be the identity map $\mathbf{y}_h^0(\Omega) = \Omega \times \{0\}$ for which $E_h^{\text{bil}} = 12.9388$ and $\mathcal{D}_h^{\text{bary}} = 3.0149 \cdot 10^{-14}$. The deformation obtained after 100, 200, 500 and 2383 (final) iterations are given in Figure 9. For the latter the energy is $E_h^{\text{bil}} = 4.2239$ while the isometry defect is $\mathcal{D}_h^{\text{bary}} = 0.0529014$.

It is worth pointing out that the equilibrium shape obtained is rather sensitive to the value of the spontaneous curvature tensors. To illustrate this, a coefficient $-0.4$ instead of $-0.3$ in $Z_2$ and $Z_4$ leads to an equilibrium deformation reported in Figure 10. After 2383 iteration of the gradient flow, the deformation is similar to the one obtained in Figure 9 (bottom-right), but the final deformation reached in 8085 iterations is quite different.
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Deformations for the bird experiment. Left to right and top to bottom: deformation obtained after 100, 200, 500, and 2383 iterations.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Deformations for the bird experiment when $Z_2 = Z_4 = [0,0;0,-0.4]$. (Left) Deformation obtained after 2383 iterations; (right) final deformation reached in 8085 iterations.</figcaption></figure>

### 5.3. Single Layer with Time-Dependent Force

In this experiment, the domain $\Omega$ is the unit disc endowed with a prestrain metric corresponding to a half-sphere, namely $g = I(y)$ with

$$ y(x_1, x_2) = \left( x_1, x_2, \sqrt{1 + \epsilon - x_1^2 - x_2^2} \right)^T, \quad \epsilon = 10^{-3}, \qquad (63) $$

where the small parameter $\epsilon$ is introduced to avoid singularities at the boundary $\partial\Omega$. We study the effect of a uniform axial force directed towards the center of the sphere. When the external force is sufficiently strong, we expect the sphere to crush like when the pressure inside a thin spherical reservoir is significantly smaller than the atmospheric pressure.

In this experiment, the external force is given by

$$ \mathbf{f}(x_1, x_2, t) = -t\mathbf{v}(x_1, x_2) = -\frac{t}{\sqrt{1+\epsilon}}\mathbf{y}(x_1, x_2), \qquad (64) $$

where $\mathbf{v}$ is the structure outward pointing normal and $t$ denotes the time. We prescribe a mixed boundary condition with deformation given by $\phi(x_1, x_2) = (x_1, x_2, \sqrt{\epsilon})^T$ on $\Gamma^M = \partial\Omega$, which is compatible with (63).

The subdivision of $\Omega$ consists of 992 quadrilaterals with hanging nodes and of maximal mesh size 0.194084 (total of 29760 degrees of freedom, 26874 = 992 · 9 · 3 for the deformation $y_h$ and 2976 = 992 · 3 for the Lagrange multipliers used to enforce the linearized metric constraint). It is finer near the boundary where more stretching is expected. The initial deformation $y_h^0$ is taken
to be the continuous Lagrange interpolant of $\mathbf{y}$ in (63), see Figure 11, for which $E_h^{\text{pre}} = 561.443$ and $\mathcal{D}_h^{\text{aver}} = 0.169394$. The other numerical parameters are $\tau = 0.2$, $\text{tol} = 10^{-3}$, and $\Delta t = 5$ so that $t_m = m\Delta t = 5m$.

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Initial deformation for the *half-sphere* experiment.</figcaption></figure>

For $m = 1, \dots, 44$, the equilibrium deformation is similar to the initial one given in Figure 11. The number of gradient flow iterations are 1594 for $m = 1$ and between 8 to 101 for $2 \le m \le 44$. When $m = 45$, $f(x_1, x_2, t_m) = -225\mathbf{v}(x_1, x_2)$ and the object cannot sustain the corresponding force anymore. It collapses and deforms all the way to a *half-sphere* in the opposite direction. In particular, the range of values for the third component of the deformation is [0.03094, 0.97704] at the first iteration and [-1.06951, 0.03254] at the last iteration of the gradient flow (81484). We refer to Figure 12 for the deformation obtained at different steps of the gradient flow and to Table 2 for the corresponding energy metric defect.

<table><caption>**Table 2.** Prestrain energy $E_h^{\text{pre}}$ and metric defect $\mathcal{D}_h^{\text{aver}}$ at different steps of the gradient flow for the *half-sphere* experiment when $t = 225$, see Figure 12 for the corresponding deformations.</caption><thead><tr><th>n</th><th>E<sub>h</sub><sup>pre</sup></th><th>\mathcal{D}<sub>h</sub><sup>aver</sup></th></tr></thead><tbody><tr><td>0</td><td>738.520</td><td>0.206666</td></tr><tr><td>400</td><td>711.633</td><td>0.245251</td></tr><tr><td>800</td><td>597.757</td><td>0.351605</td></tr><tr><td>1200</td><td>313.706</td><td>0.524551</td></tr><tr><td>81484</td><td>49.763</td><td>0.528641</td></tr></tbody></table>

### 5.4. Starshade Experiment

The starshade technology was developed as part of the NASA exoplanet exploration program. Starshades are occulters external to a telescope that shade the light from stars when imaging planets [68]. They fold to minimize the space they occupy when transported in rockets and are easily deployable when arrived at the destination.

The mathematical model corresponds to the single layer case with isometry constraint and time-dependent boundary conditions mimicking compression (closing) and decompression (opening). The geometry is taken from [69] and consists of a dodecagon, see Figure 13.

Inside the computational domain, there are 12 cubic Bézier curves as well as 6 straight lines (hexagon) across which folding is possible. The red dashed curves correspond to *valley folds* while the blue ones are *mountain folds*. The nodes of the dodecagon are given by
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Deformation for the *half-sphere* experiment when $t = 225$ in the forcing term (64). Left to right and top to bottom: initial deformation and deformation obtained after 400, 800, 1200, and 81484 steps of the gradient flow.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 13.** Domain for the *starshade* example.</figcaption></figure>

$$ \mathbf{x}_i = \left( R \cos\left(\frac{(i-1)\pi}{6}\right), R \sin\left(\frac{(i-1)\pi}{6}\right) \right), \quad i = 1, 2, \dots, 12, $$

with $R = 7$, while the nodes of the hexagon are given by

$$ \mathbf{z}_i = \left( r \cos\left(\frac{(i-1)\pi}{3} - \frac{\pi}{72}\right), r \sin\left(\frac{(i-1)\pi}{3} - \frac{\pi}{72}\right) \right), \quad i = 1, 2, \dots, 6, $$
with $r = 0.8$. For the curves, we use the following control points ($i = 1, 2, \dots, 6$):

$$ \mathbf{p}_i^V = \left( 3 \cos \left( \frac{(i-1)\pi}{3} - \frac{\pi}{8} \right), 3 \sin \left( \frac{(i-1)\pi}{3} - \frac{\pi}{8} \right) \right) $$

$$ \mathbf{q}_i^V = \left( 5 \cos \left( \frac{(i-1)\pi}{3} - \frac{\pi}{24} \right), 5 \sin \left( \frac{(i-1)\pi}{3} - \frac{\pi}{24} \right) \right) $$

$$ \mathbf{p}_i^M = \left( 3 \cos \left( \frac{(i-1)\pi}{3} + \frac{5\pi}{72} \right), 3 \sin \left( \frac{(i-1)\pi}{3} + \frac{5\pi}{72} \right) \right) $$

and

$$ \mathbf{q}_i^M = \left( 5 \cos \left( \frac{(i-1)\pi}{3} + \frac{\pi}{8} \right), 5 \sin \left( \frac{(i-1)\pi}{3} + \frac{\pi}{8} \right) \right) $$

The first valley curve is obtained using $\mathbf{z}_1$, $\mathbf{p}_1^V$, $\mathbf{q}_1^V$ and $\mathbf{x}_1$, the first mountain curve uses $\mathbf{z}_1$, $\mathbf{p}_1^M$, $\mathbf{q}_1^M$ and $\mathbf{x}_2$, the second valley curve uses $\mathbf{z}_2$, $\mathbf{p}_2^V$, $\mathbf{q}_2^V$ and $\mathbf{x}_3$, and so on.

For the boundary conditions, we compress the 6 *valley points* $\mathbf{x}_i$, $i = 1, 3, 5, 7, 9, 11$. More precisely, we consider a dynamical setting by prescribing for $i = 1, 3, 5, 7, 9, 11$, the following time-dependent pointwise boundary condition for the deformation $\mathbf{y}_h$

$$ \boldsymbol{\varphi}_i(t) = (1 - 2c_r t) \left( R \cos \left( \frac{(i-1)\pi}{6} \right), R \sin \left( \frac{(i-1)\pi}{6} \right), 0 \right)^T, \quad (65) $$

where $c_r = 0.25$ is a compression ratio. Using the strategy described in Section 4.4, we compute an approximate equilibrium deformation $\mathbf{y}_h^m$ at time $t = t_m$, $m = 0, 1, \dots, M$. Recall that this entails running the BC preprocessing for the variation $\delta \hat{\mathbf{y}}_h^m = \hat{\mathbf{y}}_h^m - \hat{\mathbf{y}}_h^{m-1}$ imposing the pointwise conditions

$$ \delta \boldsymbol{\varphi}_i = -2c_r \Delta t \left( R \cos \left( \frac{(i-1)\pi}{6} \right), R \sin \left( \frac{(i-1)\pi}{6} \right), 0 \right)^T, \quad i = 1, 3, 5, 7, 9, 11. $$

Note that at $t = 0$, $\boldsymbol{\varphi}_i = \mathbf{x}_i$ and in particular, the BC preprocessing with the above boundary conditions results in a flat deformation which is conserved by the metric preprocessing and gradient flow. In order to generate an out of plane deformation, the boundary conditions for the initialization step (see case $t = 0$ in Section 4.4) are modified to be (65) for the *valley points* and $\boldsymbol{\varphi}_i = (\mathbf{x}_i, 1.5)^T$ for the 6 *mountain points* $\mathbf{x}_i$, $i = 2, 4, 6, 8, 10, 12$.

The time-step for the time discretization is $\Delta t = 0.05$, and $M = 35$ steps are performed. The gradient flow parameters are $\tau = 0.05$ and $\text{tol} = 0.1$ while $\tilde{\tau} = 0.01$, $\tilde{\text{tol}} = 0.5$, and $\tilde{\varepsilon}_0 = 0.5$ are chosen for the metric preprocessing. The computational domain subdivision consists of 620 quadrilaterals of maximum mesh size 1.20925 (total of 18600 degrees of freedom, $16740 = 620 \cdot 9 \cdot 3$ for the deformation $\mathbf{y}_h$ and $1860 = 620 \cdot 3$ for the Lagrange multipliers used to enforce the linearized isometry constraint), see Figure 14. The penalty parameters are $\gamma_0 = \gamma_1 = \gamma_2 = 10$.
<figure><img src="image_15.png" /><figcaption>**Figure 14.** Mesh for the *starshade* example.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 15.** Left to right and top to bottom: *equilibrium* deformations for the *starshade* experiment at time $t_m$ for $m = 5, 10, 15, 20, 25, 30, 35$.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Figure 16.** Left to right: *equilibrium* deformations for the *starshade* experiment at time $t_m$ for $m = 10, 20, 30$ (side view).</figcaption></figure>

The equilibrium deformations obtained for several time $t_m = m\Delta t$ are reported in Figure 15, see Figure 16 for corresponding side views.

**Conflicts of interest**

The authors declare no competing financial interest.

**Dedication**

The manuscript was written through contributions of all authors. All authors have given approval
to the final version of the manuscript.

**References**

[1] Y. Forterre, J. Skotheim, J. Dumais, L. Mahadevan, "How the Venus flytrap Snaps", *Nature* **433** (2005), p. 421-425.

[2] R. Sachse, A. Westermeier, M. Mylo, J. Nadasdi, M. Bischoff, T. Speck, S. Poppinga, "Snapping mechanics of the Venus flytrap (*Dionaea muscipula*)", *Proc. Natl. Acad. Sci. USA* **117** (2020), no. 27, p. 16035-16042.

[3] A. Goriely, M. Ben Amar, "Differential growth and instability in elastic shells", *Phys. Rev. Lett.* **94** (2005), no. 19, article no. 198103.

[4] A. Yavari, "A geometric theory of growth mechanics", *J. Nonlinear Sci.* **20** (2010), no. 6, p. 781-830.

[5] P. Bella, R. Kohn, "Metric-Induced Wrinkling of a Thin Elastic Sheet", *J. Nonlinear Sci.* **24** (2014), p. 1147-1176.

[6] M. Lewicka, L. Mahadevan, "Geometry, analysis, and morphogenesis: Problems and prospects", *Bull. Am. Math. Soc.* **59** (2021), no. 3, p. 331-369.

[7] Y. Klein, E. Efrati, E. Sharon, "Shaping of elastic sheets by prescription of non-Euclidean metrics", *Science* **315** (2007), no. 5815, p. 1116-1120.

[8] Y. Klein, S. Venkataramani, E. Sharon, "An Experimental Study of Shape Transitions and Energy Scaling in Thin Non-Euclidean Plates", *Phys. Rev. Lett.* **106** (2011), no. 11, article no. 118303.

[9] E. Jager, E. Smela, O. Inganäs, "Microfabricating Conjugated Polymer Actuators", *Science* **290** (2000), no. 5496, p. 1540-1545.

[10] N. Bassik, B. Abebe, K. Laflin, D. Gracias, "Photolithographically patterned smart hydrogel based bilayer actuators", *Polymer* **51** (2010), no. 26, p. 6093-6098.

[11] A. Menges, S. Reichert, "Performative wood: Physically programming the responsive architecture of the HygroScope and HygroSkin projects", *Archit. Des.* **85** (2015), no. 5, p. 66-73.

[12] S. Schleicher, J. Lienhard, S. Poppinga, T. Speck, J. Knippers, "A methodology for tranferring principles of plant movements to elastic systems in architecture", *Comput. Aided Des.* **60** (2015), p. 105-117.

[13] H. Liu, P. Plucinsky, F. Feng, R. D. James, "Origami and Materials Science", *Philos. Trans. R. Soc. Lond., Ser. A* **379** (2021), no. 2201, article no. 20200113.

[14] G. Friesecke, R. D. James, S. Müller, "A Hierarchy of Plate Models Derived from Nonlinear Elasticity by Gamma-Convergence", *Arch. Ration. Mech. Anal.* **180** (2006), no. 2, p. 183-236.

[15] H. Le Dret, A. Raoult, "The Nonlinear Membrane Model as a Variational Limit of Nonlinear Three-Dimensional Elasticity", *J. Math. Pures Appl.* **73** (1995), p. 549-578.
[16] G. Friesecke, S. Müller, R. D. James, "Rigorous derivation of nonlinear plate theory and geometric rigidity", C. R. Math. Acad. Sci. Paris 334 (2002), no. 2, p. 173-178.

[17] G. Friesecke, R. D. James, S. Müller, "A theorem on geometric rigidity and the derivation of nonlinear plate theory from three-dimensional elasticity", Commun. Pure Appl. Math. 55 (2002), no. 11, p. 1461-1506.

[18] G. Friesecke, R. D. James, M. G. Mora, S. Müller, "Derivation of Nonlinear Bending Theory for Shells from Three Dimensional Nonlinear Elasticity by Gamma-Convergence", C. R. Math. Acad. Sci. Paris 336 (2003), no. 8, p. 697-702.

[19] G. Friesecke, R. D. James, S. Müller, "The Föppl-von Kármán Plate Theory as a Low Energy Γ-Limit of Nonlinear Elasticity", C. R. Math. Acad. Sci. Paris 335 (2002), no. 2, p. 201-206.

[20] M. Lewicka, M. Reza Pakzad, "Scaling laws for non-Euclidean plates and the $W^{2,2}$ isometric immersions of Riemannian metrics", ESAIM, Control Optim. Calc. Var. 17 (2011), no. 4, p. 1158-1173.

[21] G. Kirchhoff, "Über das Gleichgewicht und die Bewegung einer elastischen Scheibe", J. Reine Angew. Math. 40 (1850), p. 51-88.

[22] A. E. H. Love, A treatise on the mathematical theory of elasticity, fourth ed., Cambridge University Press, 1927.

[23] E. Efrati, E. Sharon, R. Kupferman, "Elastic theory of unconstrained non-Euclidean plates", J. Mech. Phys. Solids 57 (2009), no. 4, p. 762-775.

[24] K. Bhattacharya, M. Lewicka, M. Schäffner, "Plates with Incompatible Prestrain", Arch. Ration. Mech. Anal. 221 (2016), no. 1, p. 143-181.

[25] C. Maor, A. Shachar, "On the Role of Curvature in the Elastic Energy of Non-Euclidean Thin Bodies", J. Elasticity 134 (2019), p. 149-173.

[26] K. Böhnlein, S. Neukamm, D. Padilla-Garza, O. Sander, "A homogenized bending theory for prestrained plates", J. Nonlinear Sci. 33 (2022), no. 1, article no. 22.

[27] A. Bonito, D. Guignard, R. H. Nochetto, S. Yang, "LDG approximation of large deformations of prestrained plates", J. Comput. Phys. 448 (2022), article no. 110719.

[28] B. Schmidt, "Minimal energy configurations of strained multi-layers", Calc. Var. Partial Differ. Equ. 30 (2007), no. 4, p. 477-497.

[29] S. Bartels, A. Bonito, R. H. Nochetto, "Bilayer Plates: Model Reduction, Γ-Convergent Finite Element Approximation, and Discrete Gradient Flow", Commun. Pure Appl. Math. 70 (2017), no. 3, p. 547-589.

[30] S. Bartels, A. Bonito, P. Hornung, "Modeling and simulation of thin sheet folding", Interfaces Free Bound. 24 (2022), no. 4, p. 459-485.

[31] A. Bonito, R. H. Nochetto, S. Yang, "Γ-convergent LDG method for large bending deformations of bilayer plates", preprint arXiv:2301.03151, 2023.

[32] A. Bonito, D. Guignard, A. Morvant, "A Note on the Numerical Approximation of Thin Structures", in preparation, 2023.

[33] A. Bonito, D. Guignard, R. H. Nochetto, S. Yang, "Numerical analysis of the LDG method for large deformations of prestrained plates", IMA J. Numer. Anal. 43 (2023), no. 2, p. 627-662.

[34] A. Bonito, R. H. Nochetto, D. Ntogkas, "Discontinuous Galerkin Approach to Large Bending Deformation of a Bilayer Plate with Isometry Constraint", J. Comput. Phys. 423 (2020), article no. 109785.

[35] A. Bonito, R. H. Nochetto, D. Ntogkas, "DG approach to large bending deformations with isometry constraint", Math. Models Methods Appl. Sci. 31 (2021), no. 01, p. 133-175.

[36] S. Bartels, A. Bonito, P. Tscherner, "Error Estimates For A Linear Folding Model", IMA J. Numer. Anal. (2023).

[37] S. Bartels, "Finite element approximation of large bending isometries", Numer. Math. 124 (2013), no. 3, p. 415-440.

[38] S. Bartels, C. Palus, "Stable gradient flow discretizations for simulating bilayer plate bending with isometry and obstacle constraints", IMA J. Numer. Anal. 42 (2021), no. 3, p. 1903-1928.

[39] S. Bartels, "Finite element simulation of nonlinear bending models for thin elastic rods and plates", in Handbook of Numerical Analysis, vol. 21, Elsevier, 2020, p. 221-273.

[40] Q. Han, J.-X. Hong, Isometric embedding of Riemannian manifolds in Euclidean spaces, Mathematical Surveys and Monographs, vol. 13, American Mathematical Society, 2006.

[41] C. De Lellis, L. Székelyhidi Jr., "High dimensionality and h-principle in PDE", Bull. Am. Math. Soc. 54 (2017), p. 247-282.

[42] J. Nash, "$C^1$ Isometric Imbeddings", Ann. Math. 60 (1954), no. 3, p. 383-396.

[43] N. H. Kuiper, "On $C^1$-isometric imbeddings. I", Nederl. Akad. Wet., Proc., Ser. A 58 (1955), p. 545-556.

[44] N. H. Kuiper, "On $C^1$-isometric imbeddings. II", Nederl. Akad. Wet., Proc., Ser. A 58 (1955), p. 683-689.

[45] C. D. Modes, K. Bhattacharya, M. Warner, "Disclination-mediated thermo-optical response in nematic glass sheets", Phys. Rev. E 81 (2010), no. 6, article no. 060701.

[46] C. D. Modes, K. Bhattacharya, M. Warner, "Gaussian curvature from flat elastica sheets", Proc. R. Soc. Lond., Ser. A 467 (2010), no. 2128, p. 1121-1140.

[47] J. Kim, J. A. Hanna, R. C. Hayward, C. D. Santangelo, "Thermally responsive rolling of thin gel strips with discrete variations in swelling", Soft Matter 8 (2012), no. 8, p. 2375-2381.

[48] Z. L. Wu, M. Moshe, J. Greener, H. Therien-Aubin, Z. Nie, E. Sharon, E. Kumacheva, "Three-dimensional shape
transformations of hydrogel sheets induced by small-scale modulation of internal stresses", *Nat. Commun.* **4** (2013), article no. 1586.

[49] M. Gromov, *Partial Differential Relations*, Ergebnisse der Mathematik und ihrer Grenzgebiete. 3. Folge, vol. 13, Springer, 1986.

[50] R. Kupferman, J. P. Solomon, "A Riemannian approach to reduced plate, shell, and rod theories", *J. Funct. Anal.* **266** (2014), p. 2989-3039.

[51] E. Smela, O. Inganäs, Q. Pei, I. Lundström, "Electrochemical muscles: Micromachining fingers and corkscrews", *Adv. Mater.* **5** (1993), no. 9, p. 630-632.

[52] S. Alben, B. Balakrisnan, E. Smela, "Edge Effects Determine the Direction of Bilayer Bending", *Nano Lett.* **11** (2011), no. 6, p. 2280-2285.

[53] B. Balakrisnan, A. Nacev, E. Smela, "Design of bending multi-layer electroactive polymer actuators", *Smart Mater. Struct.* **24** (2015), no. 4, article no. 045032.

[54] D. Wood, C. Vailati, A. Menges, M. Rüggeberg, "Hygroscopically actuated wood elements for weather responsive and self-forming building parts - Facilitating upscaling and complex shape changes", *Constr. Build. Mater.* **165** (2018), p. 782-791.

[55] A. Bonito, R. H. Nochetto, "Quasi-optimal convergence rate of an adaptive discontinuous Galerkin method", *SIAM J. Numer. Anal.* **48** (2010), no. 2, p. 734-771.

[56] D. A. Di Pietro, A. Ern, *Mathematical Aspects of Discontinuous Galerkin Methods*, Mathématiques & Applications (Berlin), Springer, 2012.

[57] D. A. Di Pietro, A. Ern, "Discrete functional analysis tools for Discontinuous Galerkin methods with application to the incompressible Navier-Stokes equations", *Math. Comput.* **79** (2010), no. 271, p. 1303-1330.

[58] F. Bassi, S. Rebay, "A high-order accurate discontinuous finite element method for the numerical solution of the compressible Navier-Stokes equations", *J. Comput. Phys.* **131** (1997), no. 2, p. 267-279.

[59] D. N. Arnold, F. Brezzi, B. Cockburn, D. Marini, "Unified analysis of discontinuous Galerkin methods for elliptic problems", *SIAM J. Numer. Anal.* **39** (2002), no. 5, p. 1749-1779.

[60] F. Brezzi, G. Manzini, D. Marini, P. Pietra, A. Russo, "Discontinuous Galerkin approximations for elliptic problems", *Numer. Methods Partial Differ. Equations* **16** (2000), no. 4, p. 365-378.

[61] J. Nitsche, "Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind", *Abh. Math. Semin. Univ. Hamb.* **36** (1971), p. 9-15.

[62] Y. Nesterov, "A method for solving the convex programming problem with convergence rate $\mathcal{O}(1/k^2)$", *Dokl. Akad. Nauk SSSR* **269** (1983), no. 9, p. 543-547.

[63] Y. Nesterov, *Introductory Lectures on Convex Optimization: A Basic Course*, Applied Optimization, vol. 87, Springer, 2004.

[64] Y. Nesterov, *Lectures on convex optimization*, second ed., Springer Optimization and Its Applications, vol. 137, Springer, 2018.

[65] A. Beck, M. Teboulle, "A Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse Problems", *SIAM J. Imaging Sci.* **2** (2009), no. 1, p. 183-202.

[66] A. Chambolle, C. H. Dossal, "On the convergence of the iterates of "Fast Iterative Shrinkage/Thresholding Algorithm"", *J. Optim. Theory Appl.* **166** (2015), no. 3, p. 968-982.

[67] Y. Tahouni, T. Cheng, D. Wood, R. Sachse, R. Thierer, M. Bischoff, A. Menges, "Self-Shaping Curved Folding: A 4D-Printing Method for Fabrication of Self-Folding Curved Crease Structures", in *Symposium on Computational Fabrication*, SCF '20, ACM Press, 2020.

[68] NASA, "Starshade Technology Development", 2022, https://exoplanets.nasa.gov/exep/technology/starshade/.

[69] Jet Propulsion Laboratory - California Institute of Technology, "Space Origami: Make Your Own Starshade", 2021, https://www.jpl.nasa.gov/edu/learn/project/space-origami-make-your-own-starshade.