# Unified discrete multisymplectic Lagrangian formulation for hyperelastic solids and barotropic fluids

François Demoures, François Gay-Balmaz

► To cite this version:

François Demoures, François Gay-Balmaz. Unified discrete multisymplectic Lagrangian formulation for hyperelastic solids and barotropic fluids. Journal of Nonlinear Science, 2022, 32 (94), pp.94. 10.1007/s00332-022-09849-y. hal-03431531

HAL Id: hal-03431531

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

Submitted on 25 Jul 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Unified discrete multisymplectic Lagrangian formulation for hyperelastic solids and barotropic fluids

François Demoures<fnref n="1" />, François Gay-Balmaz²

October 4, 2021

## Abstract

We present a geometric variational discretization of nonlinear elasticity in 2D and 3D in the Lagrangian description. A main step in our construction is the definition of discrete deformation gradients and discrete Cauchy-Green deformation tensors, which allows for the development of a general discrete geometric setting for frame indifferent isotropic hyperelastic models. The resulting discrete framework is in perfect adequacy with the multisymplectic discretization of fluids proposed earlier by the authors. Thanks to the unified discrete setting, a geometric variational discretization can be developed for the coupled dynamics of a fluid impacting and flowing on the surface of an hyperelastic body. The variational treatment allows for a natural inclusion of incompressibility and impenetrability constraints via appropriate penalty terms. We test the resulting integrators in 2D and 3D with the case of a barotropic fluid flowing on incompressible rubber-like nonlinear models.

## 1 Introduction

We present a discrete variational formulation of hyperelasticity in the Lagrangian description. This formulation is associated to the multisymplectic geometric formulation of continuum mechanics which underlies both fluid dynamics and nonlinear elasticity. From the discrete variational setting, a structure preserving numerical integrator is derived for hyperelastic solids, which is spacetime multisymplectic, symplectic in time, preserves exactly the momenta associated to symmetries, and nearly preserves total energy. Thanks to its variational character, the numerical scheme can naturally incorporate incompressibility constraints. Being in perfect adequacy with the multisymplectic discretization of fluids proposed earlier in [9], the discrete framework can be extended to treat fluid-structure interaction dynamics by the inclusion of appropriate impenetrability constraints in the variational formulation.

<footnote n="1">¹EPFL, Doc & Postdoc Alumni. Av Druey 1, Lausanne 1018, Switzerland
francois.demoures@alumni.epfl.ch</footnote>

<footnote n="3">³Laboratoire de Météorologie Dynamique, École Normale Supérieure/CNRS, Paris, France.
francois.gay-balmaz@lmd.ens.fr</footnote>
The developments in this paper are founded on the multisymplectic geometric formulation of continuum mechanics, which is the spacetime version of the symplectic formulation of classical mechanics. This geometric setting allows to intrinsically formulate the multisymplectic form formula and the covariant Noether theorem, which are the field theoretic analogue to the symplecticity of the flow and to the conservation of momentum maps in classical mechanics. The more usual infinite dimensional dynamical system formulation of continuum mechanics is recovered by space and time splitting, which also has a discrete analogue.

Among the main concepts in this paper are the discrete deformation gradients, the discrete Cauchy-Green deformation tensors, and the discrete Jacobians. These objects, directly defined from the spacetime discrete configuration map of the body, allow a systematic treatment of various models of hyperelasticity via the appropriate definition of the discrete stored energy functions. The discrete variational geometric setting also systematically guides the choice of the appropriate degree of freedom of the discrete objects, for instance, at each spatial node there is one discrete Cauchy-Green tensor associated to each of the cells sharing this node. One advantage of the discrete variational framework is the possibility to include equality and inequality constraints via appropriate modification of the discrete action functional. This is developed in this paper for the treatment of fluid-structure interaction involving a barotropic fluid and an incompressible rubber-like solid. Such a coupling can be realized thanks to the unified geometric framework underlying both the solid and fluid components.

Geometric methods for discrete elasticity have been the subject of various developments. We can mention [2] concerning crystal elasticity discretization, where the discrete notions of stress and strain in lattices are introduced, and [35], [1], concerning geometric discretization of elasticity from the point of view of discrete exterior calculus.

Several points of view have been developed to treat fluid-structure interaction problems, see [20], [18], [14] for the comparison of various fluid-structure interaction methods. We can mention the *immersed boundary method* (see [27], [28]) where the interactions between the solid boundaries and the fluid are taken into account by means of local body forces, and the *fictitious domain method* (see [12]) where Lagrange multipliers are used instead of calculating the interaction forces. Other point of views to couple fluids and elastic bodies can be obtained, for example, via the *arbitrary Lagrangian Eulerian method*, see [13], techniques based on the *inverse motion function*, see [4] and [15], or the *smoothed particle hydrodynamics method*, see [17], and references therein.

Our approach is based on the geometric formulation of continuum mechanics seen as a particular classical field theory, see [11] and [23]. The discrete formulation we develop follows the variational multisymplectic discretization initially developed in [22, 19]. Such an approach has been exploited for instance towards its application to Lie group valued field theories for the treatment of geometrically exact (Cosserat) rods [8, 5] and towards its application to nonsmooth mechanics [6, 7].

The organization of the paper is as follows. Section 2 summarizes briefly the Lagrangian variational setting for elasticity on the continuous side. The elements of the constitutive continuous theory of elasticity necessary for our development are described
in rather general terms. Some aspects of the geometric formulation of field theory are also given, including the Noether theorem. In sections 3 and 4 we describe, in 2D and 3D, the discrete variational setting for elasticity. We define the discrete deformation gradients, the discrete Cauchy-Green deformation tensors, as well as the discrete Jacobians, in a parallel way with the continuous setting. Based on these definitions, a discrete Lagrangian is constructed for frame indifferent isotropic materials. It is then shown how this approach can be naturally coupled with the variational discretization of barotropic fluids thanks to the introduction of appropriate penalty functions associated to the impenetrability constraints. Illustrative numerical examples are given in 2D and 3D for a barotropic fluid flowing into an hyperelastic container described with the St. Venant-Kirchhoff or the Mooney-Rivlin model.

## 2 Elasticity and field theory

In this section we briefly review the variational geometric setting for elasticity in the Lagrangian (or material) description by focusing on isotropic hyperelastic materials. The description is given in terms of the (right) Cauchy-Green tensor, which plays a central role in this paper. The formulation is then recasted in a multisymplectic variational setting, which allows to formulate intrinsically the Hamilton principle, the multisymplectic property of the solutions, and the covariant Noether theorem associated to symmetries. This gives the geometric framework to be discretized in a structure preserving way later.

### 2.1 Description of hyperelastic bodies

We assume that the elastic body is represented by a compact manifold $\mathcal{B} \subset \mathbb{R}^3$ with piecewise smooth boundary, and moves in the ambient space $\mathcal{M} = \mathbb{R}^3$. The motion of the body is described by a time dependent configuration map $\varphi : \mathbb{R} \times \mathcal{B} \to \mathcal{M}$ which indicates the location $m = \varphi(t, X) \in \mathcal{M}$ at time $t$ of the material point $X \in \mathcal{B}$. The deformation gradient is denoted $\mathbf{F}(t, X) = \nabla\varphi(t, X)$, given in coordinates by $\mathbf{F}^a_i = \varphi^a_{,i}$, with $X^i, i = 1, 2, 3$ the Cartesian coordinates on $\mathcal{B}$ and $m^a, a = 1, 2, 3$ the Cartesian coordinates on $\mathcal{M}$. We assume that the configuration map is regular enough so that all the computations below are valid.

A material is called hyperelastic if there is a stored energy function $\mathcal{W}$ depending on the points $X \in \mathcal{B}$ and the gradient deformation $\mathbf{F}(t, X)$ such that the first Piola-Kirchhoff stress tensor of the body in the configuration $\varphi(t, X)$ is $\mathbf{P} = \rho_0 \frac{\partial \mathcal{W}}{\partial \mathbf{F}}$, where $\rho_0$ is the mass density of the body in the Lagrangian description. Given Riemannian metrics $\mathbf{G}$ and $\mathbf{g}$ on $\mathcal{B}$ and $\mathcal{M}$, the stored energy function of an hyperelastic material with configuration map $\varphi : \mathcal{B} \to \mathcal{M}$ (at some fixed time) takes the general form

$$ \mathcal{W}(X, \mathbf{F}(X), \mathbf{G}(X), \mathbf{g}(\varphi(X))). $$

This expression of the stored energy function explicitly written in terms of Riemannian metrics allows to naturally reformulate the axiom of material frame indifference and the
isotropic property as invariances of W with respect to spatial and material diffeomorphisms, i.e., diffeomorphisms of M and B, [21].

From the axiom of material frame indifference the stored energy function depends on the deformation gradient only through the Cauchy-Green deformation tensor defined by

$$ \mathbf{C} = \varphi^*\mathbf{g}, \quad \mathbf{C}_{ij}(X) = \mathbf{g}_{ab}(\varphi(X))\varphi^{a}_{,i}(X)\varphi^{b}_{,j}(X), \quad (1) $$

i.e., we have $W(X, \mathbf{F}(X), \mathbf{G}(X), \mathbf{g}(\varphi(X))) = W(X, \mathbf{G}(X), \mathbf{C}(X))$, for some function W, where $\varphi^*$ denotes the pull-back of a tensor field by $\varphi$.

The material (or Lagrangian) strain tensor E and the second Piola-Kirchhoff stress tensor S are defined respectively by

$$ \mathbf{E}(\mathbf{C}) = \frac{1}{2}(\mathbf{C} - \mathbf{G}), \quad \mathbf{E}_{ij} = \frac{1}{2}(\mathbf{C}_{ij} - \mathbf{G}_{ij}), \quad (2) $$

$$ \mathbf{S}(\mathbf{C}) = 2\rho_0 \frac{\partial W}{\partial \mathbf{C}}, \quad \mathbf{S}^{ij} = 2\rho_0 \frac{\partial W}{\partial \mathbf{C}_{ij}}. \quad (3) $$

The elasticity tensor on the reference configuration is defined by

$$ \mathbf{C} = \frac{\partial \mathbf{S}}{\partial \mathbf{C}}, \quad \mathbf{C}^{ijkl} = \frac{\partial \mathbf{S}^{ij}}{\partial \mathbf{C}_{kl}} = 2\rho_0 \frac{\partial^2 W}{\partial \mathbf{C}_{ij} \partial \mathbf{C}_{kl}}. \quad (4) $$

For simplicity, we will consider the Euclidean case, where $\mathbf{G}_{ij} = \delta_{ij}$ and $\mathbf{g}_{ab} = \delta_{ab}$ in Cartesian coordinates, in which case we can write

$$ \mathbf{C} = \mathbf{F}^T\mathbf{F} \quad \text{and} \quad \mathbf{C}_{ij} = \varphi^{a}_{,i}\varphi^{a}_{,j} = \langle \varphi_{,i}, \varphi_{,j} \rangle, $$

where $\langle \cdot, \cdot \rangle$ denotes the Euclidean inner product on $\mathbb{R}^n$.

## 2.2 Euler-Lagrange equations for an hyperelastic body

The Lagrangian density of an hyperelastic body evaluated on a configuration map $\varphi(t, X)$ is

$$ \begin{aligned} \mathcal{L}(\varphi, \dot{\varphi}, \nabla\varphi) &= L(\varphi, \dot{\varphi}, \nabla\varphi)d^3X \wedge dt \\ &= \left[ \frac{1}{2}\rho_0 \langle \dot{\varphi}, \dot{\varphi} \rangle - \rho_0 W(\mathbf{G}, \mathbf{C}) - \rho_0 \Pi(\varphi) \right] d^3X \wedge dt, \end{aligned} \quad (5) $$

with $\Pi$ a potential energy term such as the gravitational potential. Note that we didn't write explicitly the dependence on X and t. Hamilton's principle reads

$$ \delta \int_{0}^{T} \int_{B} L(\varphi, \dot{\varphi}, \nabla\varphi) d^3X \wedge dt = 0, $$

for variations $\delta\varphi$ of the body configuration map with $\delta\varphi = 0$ at $t = 0, T$. It yields the Euler-Lagrange equations

$$ \frac{d}{dt} \frac{\partial L}{\partial \dot{\varphi}} - \frac{\partial L}{\partial \varphi} = - \mathrm{DIV} \frac{\partial L}{\partial \mathbf{F}} \quad (6) $$
together with the natural boundary conditions

$$ \frac{\partial L}{\partial \mathbf{F}_i^a} n_i \delta \varphi^a = 0 \quad \text{on} \quad \partial \mathcal{B} \qquad (7) $$

for the allowed variations $\delta\varphi$ at the boundary, with $\mathbf{F} = \nabla\varphi$. For the Lagrangian (5) we compute

$$ \frac{\partial L}{\partial \mathbf{F}_i^a} = -2\rho_0 \mathbf{F}_j^a \frac{\partial W}{\partial C_{ji}} = -\mathbf{F}_j^a \mathbf{S}^{ji} \quad \text{i.e.} \quad \frac{\partial L}{\partial \mathbf{F}} = -\mathbf{F} \mathbf{S}, $$

hence the Euler-Lagrange equations (6) for the Lagrangian (5) gives the equations of motion for an hyperelastic body as

$$ \rho_0 \ddot{\varphi} = \text{DIV}(\mathbf{F}\mathbf{S}) - \rho_0 \frac{\partial \Pi}{\partial \varphi}. \qquad (8) $$

In this paper we focus exclusively on the variational formulation and the equations in the material description, see [10] for the variational formulation in the Eulerian and convective descriptions and its relation with material frame indifference and material covariance discussed above.

## 2.3 Isotropic hyperelastic materials

A frame indifferent hyperelastic material is isotropic if and only if the stored energy function can be written as $W(X, \mathbf{G}(X), \mathbf{C}(X)) = \Phi(X, \lambda_1(X), \lambda_2(X), \lambda_3(X))$, where $\Phi$ is a symmetric function of the principal stretches $\lambda_1, \lambda_2, \lambda_3$. This is equivalent to state that $W$ is a function of the invariants $I_1, I_2, I_3$ of $\mathbf{C}$ defined by$^2$

$$ \begin{aligned} I_1(\mathbf{C}) &= \text{Tr}(\mathbf{C}), & I_2(\mathbf{C}) &= \frac{1}{2} \left[ (I_1(\mathbf{C}))^2 - \text{Tr}(\mathbf{C}^2) \right] & I_3(\mathbf{C}) &= \det(\mathbf{C}). \end{aligned} \qquad (9) \\ = \det(\mathbf{C})\text{Tr}(\mathbf{C}^{-1}), $$

The derivative of the invariants with respect to $\mathbf{C}$ are given by

$$ \frac{\partial I_1}{\partial \mathbf{C}} = \mathbf{G}^{-1}, \quad \frac{\partial I_2}{\partial \mathbf{C}} = I_2 \mathbf{C}^{-1} - I_3 \mathbf{C}^{-2}, \quad \frac{\partial I_3}{\partial \mathbf{C}} = I_3 \mathbf{C}^{-1}, \qquad (10) $$

from which the second Piola-Kirchhoff stress tensor (3) reads

$$ \mathbf{S} = 2\rho_0 \left[ \frac{\partial W}{\partial I_1} \mathbf{G}^{-1} + \left( \frac{\partial W}{\partial I_2} I_2 + \frac{\partial W}{\partial I_3} I_3 \right) \mathbf{C}^{-1} - \frac{\partial W}{\partial I_2} I_3 \mathbf{C}^{-2} \right]. \qquad (11) $$

We refer to [21, §3.5] for more details.

<footnote>$^1$Given the polar decomposition $\mathbf{F} = \mathbf{R}\mathbf{U}$ of the gradient deformation, the principal stretches $\lambda_1, \lambda_2, \lambda_3$ are the eigenvalues of $\mathbf{U}$. In particular $\lambda_1^2, \lambda_2^2, \lambda_3^2$ are the eigenvalues of $\mathbf{C} = \mathbf{U}^2$.</footnote>

<footnote>$^2$For simplicity we didn't write explicitly the dependence of the invariants on $\mathbf{G}$ and assumed $\mathbf{G}_{ij} = \delta_{ij}$.</footnote>
The Lagrangian density of a frame indifferent isotropic hyperelastic material is therefore

$$ \mathcal{L}(\varphi, \dot{\varphi}, \nabla \varphi) = \left[ \frac{1}{2} \rho_0 (\langle \dot{\varphi}, \dot{\varphi} \rangle - \rho_0 W(I_1, I_2, I_3) - \rho_0 \Pi(\varphi)) \right] d^3 X \wedge dt, \quad (12) $$

while the equations of motion take the form

$$ \rho_0 \ddot{\varphi} = \text{DIV} \left[ 2\rho_0 \mathbf{F} \left( \frac{\partial W}{\partial I_1} \mathbf{G}^{-1} + \left( \frac{\partial W}{\partial I_2} I_2 + \frac{\partial W}{\partial I_3} I_3 \right) \mathbf{C}^{-1} - \frac{\partial W}{\partial I_2} I_3 \mathbf{C}^{-2} \right) \right] - \rho_0 \frac{\partial \Pi}{\partial \varphi}. \quad (13) $$

Incompressible materials can be treated by inserting the constraint $I_3 = 1$ in the Hamilton principle as

$$ \delta \int_0^T \int_B (L(\varphi, \dot{\varphi}, \nabla \varphi) + \mu(I_3 - 1)) d^3 X \wedge dt = 0, $$

where $\mu$ is a Lagrange multiplier.

**Example 2.1** A well-known example of isotropic hyperelastic model is the *Moonley-Rivlin* model ([25], [29, 30, 31]), which is suitable for the description of certain incompressible rubber-like materials. Its stored energy function is defined by the simple expression

$$ W(I_1, I_2) = C_1(I_1 - 3) + C_2(I_2 - 3), \quad (14) $$

with positive material constants $C_1, C_2$.

## 2.4 Multisymplectic variational continuum mechanics

We briefly recall here the geometric variational framework of classical field theory, as it applies to continuum mechanics (fluid and elasticity), see [23] and [19]. In §3 and §4 we shall carry out a structure preserving discretization of this setting which allows the identification of the notion of discrete multisymplecticity, discrete momentum map, and discrete Noether theorems.

**Configuration bundle, jet bundle, and Lagrangian density.** The geometric formulation of classical field theories starts with the identification of the configuration bundle of the theory, denoted $\pi_{y,X}: \mathcal{Y} \to \mathcal{X}$. The configuration fields $\varphi$ of the theory are sections of this fiber bundle, i.e., they are smooth maps $\varphi: \mathcal{X} \to \mathcal{Y}$ such that $\pi_{y,X} \circ \varphi = \text{id}_{\mathcal{X}}$, where $\text{id}_{\mathcal{X}}$ denotes the identity map on $\mathcal{X}$.

The first jet bundle of the configuration bundle $\pi_{y,X}: \mathcal{Y} \to \mathcal{X}$ is the field theoretic analogue of the tangent bundle of classical mechanics. It is defined as the fiber bundle $\pi_{J^1y,y}: J^1\mathcal{Y} \to \mathcal{Y}$ whose fiber at $y \in \mathcal{Y}$ consists of linear maps $\gamma: T_x\mathcal{X} \to T_y\mathcal{Y}$ satisfying $T\pi_{y,X} \circ \gamma = \text{id}_{T_x\mathcal{X}}$, where $x = \pi_{y,X}(y)$. The derivative of a field $\varphi$ can be regarded as a section of the fiber bundle $\pi_{J^1y,X} = \pi_{y,X} \circ \pi_{J^1y,y}: J^1\mathcal{Y} \to \mathcal{X}$, by considering its first jet extension $x \in \mathcal{X} \to j^1\varphi(x) = T_x\varphi \in J^1_{\varphi(x)}\mathcal{Y}$, with $T_x\varphi: T_x\mathcal{X} \to T_{\varphi(x)}\mathcal{Y}$ the tangent map of $\varphi$.
A Lagrangian density is a smooth bundle map $\mathcal{L}: J^1\mathcal{Y} \to \Lambda^{n+1}\mathcal{X}$ over $\mathcal{X}$, where $\Lambda^{n+1}\mathcal{X} \to \mathcal{X}$ is the vector bundle of $(n+1)$-forms on $\mathcal{X}$, with $\dim \mathcal{X} = n+1$. The associated action functional is given by

$$S(\varphi) := \int_{\mathcal{X}} \mathcal{L}(j^1\varphi(x)), \qquad (15)$$

and the Euler-Lagrange equations follow from the stationary condition for $\delta S(\varphi) = 0$ for appropriate boundary conditions.

For continuum mechanics, the configuration bundle is the trivial fiber bundle

$$\mathcal{Y} = \mathcal{M} \times \mathcal{X} \to \mathcal{X}, \quad \text{with} \quad \mathcal{X} = \mathbb{R} \times \mathcal{B},$$

where $\mathcal{B}$ is the reference configuration of the elastic body and $\mathcal{M}$ is the ambient space. We have the variables $x = (t, X) \in \mathcal{X}$ and $y = (x, m) = (t, X, m) \in \mathcal{Y}$. A section of this bundle is a map $\varphi: \mathcal{X} \to \mathcal{X} \times \mathcal{M}$, whose first component is $id_{\mathcal{X}}$. It is canonically identified with a map $\varphi: \mathcal{X} = \mathbb{R} \times \mathcal{B} \to \mathcal{M}$ referred to as the elastic body configuration in §2.1. The first jet extension of $\varphi$ is $j^1\varphi(t, X) = (\varphi(t, X), \dot{\varphi}(t, X), \nabla\varphi(t, X))$ and the Lagrangian density reads

$$\mathcal{L}(\varphi, \dot{\varphi}, \nabla\varphi) = L(\varphi, \dot{\varphi}, \nabla\varphi)d^3X \wedge dt.$$

**Multisymplectic form formula and Noether theorem.** A solution $\varphi$ of the Euler-Lagrange equation associated to $\mathcal{L}: J^1\mathcal{Y} \to \Lambda^{n+1}\mathcal{X}$ satisfies the multisymplectic form formula given by

$$\int_{\partial U} (j^1 \varphi)^* \dot{\mathbf{i}}_{j^1 V} \mathbf{i}_{j^1 W} \Omega_{\mathcal{L}} = 0, \qquad (16)$$

for all open subset $U \subset \mathcal{X}$ with piecewise smooth boundary, and for all vector fields $V, W$ solutions of the first variation of the Euler-Lagrange equations at $\varphi$, see [22]. In this formula, $\Omega_{\mathcal{L}}$ denotes the Cartan $(n+2)$-form on $J^1\mathcal{Y}$ associated to $\mathcal{L}$ and $j^1V, j^1W$ are the first jet extension of $V, W$ to $J^1\mathcal{Y}$, see [11]. Property (16) is an extension to field theory of the symplectic property of the solution of the Euler-Lagrange equations of classical mechanics.

Consider the action of a Lie group $G$ on $\mathcal{Y}$, assume that this action covers a diffeomorphism of $\mathcal{X}$, and assume that the Lagrangian density $\mathcal{L}$ is $G$-equivariant. Then a solution $\varphi$ of the Euler-Lagrange equation satisfies the covariant Noether theorem

$$\int_{\partial U} (j^1 \varphi)^* J^\mathcal{L}(\xi) = 0, \qquad (17)$$

for all open subset $U \subset \mathcal{X}$ with piecewise smooth boundary and for all $\xi \in \mathfrak{g}$, the Lie algebra of $G$, with $J^\mathcal{L}: J^1\mathcal{Y} \to \mathfrak{g}^* \otimes \Lambda^n J^1\mathcal{Y}$ the covariant momentum map associated to $\mathcal{L}$ and to the Lie group action, see [11] for details.

The properties (16) and (17) admit discrete versions that will be considered later.
## 3 2D discrete fluid - elastic body interactions

In this section we develop a multisymplectic variational discretization of two dimensional hyperelasticity. A main step in our approach is the definition of the discrete Cauchy-Green tensor associated to a discrete configuration map. We then show how this approach can be naturally coupled with the multisymplectic discretization of barotropic fluids. This is then used for the development of a variational discretization of the coupled dynamics of a barotropic fluid flowing on a hyperelastic body.

### 3.1 Multisymplectic discretization and the discrete Cauchy-Green tensor

**Discrete Lagrangian setting.** Consider the configuration bundle of continuum mechanics $\mathcal{Y} = \mathcal{X} \times \mathcal{M} \to \mathcal{X} = \mathbb{R} \times \mathcal{B}$. In this section we assume that $\mathcal{B}$ is a domain in $\mathbb{R}^2$ and take $\mathcal{M} = \mathbb{R}^2$. To discretize the geometric setting, one first considers a discrete parameter space $\mathcal{U}_d$ and a discrete base-space configuration, which is a one-to-one map

$$ \phi_{\mathcal{X}_d}: \mathcal{U}_d \to \phi_{\mathcal{X}_d}(\mathcal{U}_d) = \mathcal{X}_d \subset \mathcal{X} $$

whose image $\mathcal{X}_d$ is the discrete spacetime. The discrete fields are the sections of the discrete configuration bundle $\pi_d: \mathcal{Y}_d = \mathcal{X}_d \times \mathcal{M} \to \mathcal{X}_d$, which can be identified with maps $\varphi_d: \mathcal{X}_d \to \mathcal{M}$. The discrete field and the discrete base-space configuration can be simultaneously described by introducing the discrete configuration, which is a map $\phi_d: \mathcal{U}_d \to \mathcal{Y}$. From $\phi_d$ we obtain the discrete base-space configuration and the discrete field as $\phi_{\mathcal{X}_d} = \pi_d \circ \phi_d$ and $\varphi_d = \phi_d \circ \phi_{\mathcal{X}_d}^{-1}$, see Fig. 1.

<figure><img src="image_2.png" /><figcaption>Figure 1: Discrete configuration and discrete configuration bundle</figcaption></figure>

We choose the discrete parameter space $\mathcal{U}_d = \{0, ..., j, ..., N\} \times \mathbb{B}_d$. Here $\{0, ..., j, ..., N\}$ describes an increasing sequence of time and $\mathbb{B}_d$ encodes the nodes and simplexes of the discretization of $\mathcal{B}$. Assuming for simplicity that $\mathcal{B}$ is a rectangle, we consider $\mathbb{B}_d = \{0, ..., A\} \times \{0, ..., B\}$ and denote $(j, a, b) \in \mathcal{U}_d$ the element of the discrete parameter space. The latter determines a set of parallelepipeds, denoted $\mathbb{B}_{a,b}^j$, defined by the following eight pairs of indices

$$ \mathbb{B}_{a,b}^j = \{(j, a, b), (j+1, a, b), (j, a+1, b), (j, a, b+1), (j, a+1, b+1), \\ (j+1, a+1, b), (j+1, a, b+1), (j+1, a+1, b+1)\}, \quad (18) $$
with $j = 0, ..., N-1$, $a = 0, ..., A-1$, $b = 0, ..., B-1$, see Fig. 2. We assume that the discrete base-space configuration is of the form

$$ \phi_{\chi_d}: \mathcal{U}_d \ni (j, a, b) \mapsto s_{a,b}^j = (t^j, s_a^j, s_b^j) \in \mathcal{X}_d \quad (19) $$

and denote by $\varphi_{a,b}^j := \varphi_d(s_{a,b}^j)$ the value of the discrete field at $s_{a,b}^j$.

<figure><img src="image_3.png" /><figcaption>Figure 2: Discrete spacetime domain $\mathcal{U}_d$.</figcaption></figure>

Denoting by $\mathcal{U}_d^{\mathcal{D}}$ the set of all parallelepipeds in the discrete parameter space $\mathcal{U}_d$, we write $\chi_d^{\mathcal{D}} := \phi_{\chi_d}(\mathcal{U}_d^{\mathcal{D}})$ and recall that the discrete version of the first jet bundle is given by

$$ J^1 \mathcal{Y}_d := \chi_d^{\mathcal{D}} \times \underbrace{\mathcal{M} \times \dots \times \mathcal{M}}_{8 \text{ times}} \to \chi_d^{\mathcal{D}}. \quad (20) $$

Given a discrete field $\varphi_d$, its first jet extension is the section of (20) defined by

$$ j^1 \varphi_d(\mathcal{D}_{a,b}^j) = (\varphi_{a,b}^j, \varphi_{a,b}^{j+1}, \varphi_{a+1,b}^j, \varphi_{a+1,b}^{j+1}, \varphi_{a,b+1}^j, \varphi_{a,b+1}^{j+1}, \varphi_{a+1,b+1}^j, \varphi_{a+1,b+1}^{j+1}), \quad (21) $$

which associates to each parallelepiped, the values of the field at its nodes. Finally, we recall that a discrete Lagrangian is a map

$$ \mathcal{L}_d: J^1 \mathcal{Y}_d \to \mathbb{R}, $$

from which the discrete action functional is constructed as

$$ S_d(\varphi_d) = \sum_{\mathcal{D} \in \chi_d^{\mathcal{D}}} \mathcal{L}_d(j^1 \varphi_d(\mathcal{D})). $$

**Discrete Cauchy-Green deformation tensor.** Given a discrete base-space configuration $\phi_{\chi_d}$ and a discrete field $\varphi_d$, the following four vectors $\mathbf{F}_{\ell;a,b}^j \in \mathbb{R}^2$, $\ell = 1, 2, 3, 4$ are defined at each node $(j, a, b) \in \mathcal{U}_d$, see Fig. 3 on the right:

$$ \mathbf{F}_{1;a,b}^j = \frac{\varphi_{a+1,b}^j - \varphi_{a,b}^j}{|s_{a+1}^j - s_a^j|} \quad \text{and} \quad \mathbf{F}_{2;a,b}^j = \frac{\varphi_{a,b+1}^j - \varphi_{a,b}^j}{|s_{b+1}^j - s_b^j|} \quad (22) $$

$$ \mathbf{F}_{3;a,b}^j = \frac{\varphi_{a-1,b}^j - \varphi_{a,b}^j}{|s_a^j - s_{a-1}^j|} = -\mathbf{F}_{1;a-1,b}^j \quad \text{and} \quad \mathbf{F}_{4;a,b}^j = \frac{\varphi_{a,b-1}^j - \varphi_{a,b}^j}{|s_b^j - s_{b-1}^j|} = -\mathbf{F}_{2;a,b-1}^j. $$
<figure><img src="image_4.png" /><figcaption>Figure 3: Discrete field $\phi_d = \varphi_d \circ \phi_{\chi_d}$ evaluated on $\mathcal{D}_{a,b}^j, \mathcal{D}_{a,b}^j, \mathcal{D}_{a,b}^j, \mathcal{D}_{a,b}^j$ at time $t^j$ and the associated vectors $\mathbf{F}_{\ell; a, b}, \ell = 1, 2, 3, 4$.</figcaption></figure>

**Definition 3.1** The discrete deformation gradients of a discrete field $\varphi_d$ at the parallelepiped $\mathcal{D}_{a,b}^j$ are the four $2 \times 2$ matrices $\mathbf{F}_{\ell}(\mathcal{D}_{a,b}^j), \ell = 1, 2, 3, 4$, defined by

$$ \begin{aligned} \mathbf{F}_1(\mathcal{D}_{a,b}^j) &= [\mathbf{F}_{1;a,b}^j \quad \mathbf{F}_{2;a,b}^j], & \mathbf{F}_2(\mathcal{D}_{a,b}^j) &= [\mathbf{F}_{2;a+1,b}^j \quad \mathbf{F}_{3;a+1,b}^j], \\ \mathbf{F}_3(\mathcal{D}_{a,b}^j) &= [\mathbf{F}_{4;a,b+1}^j \quad \mathbf{F}_{1;a,b+1}^j], & \mathbf{F}_4(\mathcal{D}_{a,b}^j) &= [\mathbf{F}_{3;a+1,b+1}^j \quad \mathbf{F}_{4;a+1,b+1}^j]. \end{aligned} \tag{23} $$

Note that the discrete deformation gradients are defined at each of the four nodes of $\mathcal{D}_{a,b}^j$ that are associated to time $t^j$, see Fig. 3 on the left for the ordering $\ell = 1, 2, 3, 4$ with respect to the nodes. Also, for each $\mathcal{D}$, the discrete deformation gradients $\mathbf{F}_{\ell}(\mathcal{D})$, $\ell = 1, 2, 3, 4$, depend on the discrete field $\varphi_d$ only through its first jet extension $j^1\varphi_d(\mathcal{D})$ at the given $\mathcal{D}$, see (21).

It is assumed that the discrete field $\varphi_d$ is such that the determinant of the discrete deformation gradient are positive. Using Definition 3.1, the discrete Cauchy-Green deformation tensors can be defined exactly as in the continuous case.

**Definition 3.2** The discrete Cauchy-Green deformation tensors of a discrete field $\varphi_d$ at the parallelepiped $\mathcal{D}_{a,b}^j$ are the four $2 \times 2$ symmetric matrices

$$ \mathbf{C}_{\ell}(\mathcal{D}_{a,b}^j) = \mathbf{F}_{\ell}(\mathcal{D}_{a,b}^j)^{\top} \mathbf{F}_{\ell}(\mathcal{D}_{a,b}^j). \tag{24} $$

With the assumption on the discrete field mentioned above, the discrete deformation gradients and the discrete Cauchy-Green deformation tensors are maps

$$ \mathbf{F}_{\ell}: \mathcal{X}_{d}^{\mathcal{D}} \to GL(2), \quad \mathbf{C}_{\ell}: \mathcal{X}_{d}^{\mathcal{D}} \to \mathrm{Sym}_{+}(2), \quad \ell = 1, 2, 3, 4, $$

defined on $\mathcal{X}_{d}^{\mathcal{D}}$, with $\mathrm{Sym}_{+}(2)$ the set of symmetric positive definite $2 \times 2$ matrices. By using the definitions (23) and (24), the discrete Cauchy-Green deformation tensors can be explicitly written in terms of the vectors defined in (22) as follows, which is a discrete counterpart of the continuous relation $\mathbf{C}_{ij} = \langle \varphi_{,i}, \varphi_{,j} \rangle$.
**Proposition 3.3** Given a discrete field $\varphi_d: \mathcal{X}_d \to \mathcal{M}$, the discrete Cauchy-Green deformation tensors are explicitly given as

$$ \mathbf{C}_1(\mathbb{D}_{a,b}^j) = \begin{bmatrix} \langle \mathbf{F}_{1;a,b}^j, \mathbf{F}_{1;a,b}^j \rangle & \langle \mathbf{F}_{1;a,b}^j, \mathbf{F}_{2;a,b}^j \rangle \\ \langle \mathbf{F}_{2;a,b}^j, \mathbf{F}_{1;a,b}^j \rangle & \langle \mathbf{F}_{2;a,b}^j, \mathbf{F}_{2;a,b}^j \rangle \end{bmatrix}, \qquad (25) $$

$$ \mathbf{C}_2(\mathbb{D}_{a,b}^j) = \begin{bmatrix} \langle \mathbf{F}_{3;a+1,b}^j, \mathbf{F}_{3;a+1,b}^j \rangle & \langle \mathbf{F}_{3;a+1,b}^j, \mathbf{F}_{2;a+1,b}^j \rangle \\ \langle \mathbf{F}_{2;a+1,b}^j, \mathbf{F}_{3;a+1,b}^j \rangle & \langle \mathbf{F}_{2;a+1,b}^j, \mathbf{F}_{2;a+1,b}^j \rangle \end{bmatrix}, \qquad (26) $$

$$ \mathbf{C}_3(\mathbb{D}_{a,b}^j) = \begin{bmatrix} \langle \mathbf{F}_{1;a,b+1}^j, \mathbf{F}_{1;a,b+1}^j \rangle & \langle \mathbf{F}_{1;a,b+1}^j, \mathbf{F}_{4;a,b+1}^j \rangle \\ \langle \mathbf{F}_{4;a,b+1}^j, \mathbf{F}_{1;a,b+1}^j \rangle & \langle \mathbf{F}_{4;a,b+1}^j, \mathbf{F}_{4;a,b+1}^j \rangle \end{bmatrix}, \qquad (26) $$

$$ \mathbf{C}_4(\mathbb{D}_{a,b}^j) = \begin{bmatrix} \langle \mathbf{F}_{3;a+1,b+1}^j, \mathbf{F}_{3;a+1,b+1}^j \rangle & \langle \mathbf{F}_{3;a+1,b+1}^j, \mathbf{F}_{4;a+1,b+1}^j \rangle \\ \langle \mathbf{F}_{4;a+1,b+1}^j, \mathbf{F}_{3;a+1,b+1}^j \rangle & \langle \mathbf{F}_{4;a+1,b+1}^j, \mathbf{F}_{4;a+1,b+1}^j \rangle \end{bmatrix}, $$

in terms of the vector $\mathbf{F}_{\ell;a,b}^j$ defined in (22).

**Discrete Jacobian.** We recall from [9] that the discrete Jacobian of a discrete field $\varphi_d$ at the parallelepiped $\mathbb{D}_{a,b}^j$ are the four numbers $J_\ell(\mathbb{D}_{a,b}^j)$, $\ell = 1, 2, 3, 4$, defined by

$$ \begin{aligned} J_1(\mathbb{D}_{a,b}^j) &= |\mathbf{F}_{1;a,b}^j \times \mathbf{F}_{2;a,b}^j|, & J_2(\mathbb{D}_{a,b}^j) &= |\mathbf{F}_{2;a+1,b}^j \times \mathbf{F}_{3;a+1,b}^j|, \\ J_3(\mathbb{D}_{a,b}^j) &= |\mathbf{F}_{4;a,b+1}^j \times \mathbf{F}_{1;a,b+1}^j|, & J_4(\mathbb{D}_{a,b}^j) &= |\mathbf{F}_{3;a+1,b+1}^j \times \mathbf{F}_{4;a+1,b+1}^j|. \end{aligned} \qquad (27) $$

The discrete Jacobian can be written in terms of the discrete deformation gradient and the discrete Cauchy-Green tensors exactly as in the continuous case, namely, we have the following relations.

**Proposition 3.4** The discrete Jacobian is related to the discrete deformation gradients and discrete Cauchy-Green tensors as follows

$$ \det(\mathbf{F}_{\ell}(\mathbb{D}_{a,b}^j)) = J_{\ell}(\mathbb{D}_{a,b}^j), \quad \det(\mathbf{C}_{\ell}(\mathbb{D}_{a,b}^j)) = J_{\ell}(\mathbb{D}_{a,b}^j)^2, \qquad (28) $$

for each $\mathbb{D}_{a,b}^j \in \mathcal{X}_d^\mathbb{D}$ and all $\ell = 1, 2, 3, 4$.

## 3.2 Discrete frame indifferent linear hyperelastic model

In the 2D case, we shall focus on the St. Venant-Kirchhoff model

$$ \rho_0 W(\mathbf{C}) = \lambda \operatorname{Tr}(\mathbf{E})^2 + 2\mu \operatorname{Tr}(\mathbf{E}^2) = \mathbf{E}^\mathrm{T} \mathbf{C} \mathbf{E}, \qquad (29) $$

where $\lambda$ and $\mu$ are the first and second Lamé coefficients and $\mathbf{E}$ is the material strain tensor defined in (2) with $\mathbf{G}_{ij} = \delta_{ij}$. In the last term we rewrote the stored energy
function in a vectorial form by introducing the vector E and the symmetric matrix C defined by

$$ E = (\mathbf{E}_{11} \ \mathbf{E}_{22} \ 2\mathbf{E}_{12})^{\mathrm{T}} \in \mathbb{R}^3 \quad \text{and} \quad \mathbf{C} = \frac{E}{1 - \nu^2} \begin{bmatrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & \frac{1-\nu}{2} \end{bmatrix} $$

with $\nu = \frac{\lambda}{\lambda+2\mu}$ the Poisson ratio and $E = \frac{4\mu(\lambda+\mu)}{\lambda+2\mu}$ the Young modulus in 2D.

**Discrete stored energy function.** The discrete stored energy function $W_d: J^1\mathcal{Y}_d \to \mathbb{R}$ is defined in terms of the continuous one as follows

$$ \rho_0 W_d(j^1 \varphi_d(\mathbf{u})) := \frac{\rho_0}{4} \sum_{\ell=1}^{4} W(\mathbf{C}_\ell(\mathbf{u})) = \frac{1}{4} \sum_{\ell=1}^{4} \mathbf{E}_\ell(\mathbf{u})^\mathrm{T} \mathbf{C} \mathbf{E}_\ell(\mathbf{u}), \quad (30) $$

with $\mathbf{E}_\ell(\mathbf{u}) = (\mathbf{E}_{\ell 11}(\mathbf{u}) \ \mathbf{E}_{\ell 22}(\mathbf{u}) \ 2\mathbf{E}_{\ell 12}(\mathbf{u}))^\mathrm{T}$. For instance, for $\ell = 1$, we have

$$ \mathbf{E}_1(\mathbf{u}_{a,b}^j) = \left( \frac{1}{2} (\langle \mathbf{F}_{1;a,b}, \mathbf{F}_{1;a,b} \rangle - 1), \frac{1}{2} (\langle \mathbf{F}_{2;a,b}, \mathbf{F}_{2;a,b} \rangle - 1), \langle \mathbf{F}_{1;a,b}, \mathbf{F}_{2;a,b} \rangle \right). $$

**Discrete second Piola-Kirchhoff stress.** The discrete second Piola-Kirchhoff stress at time $t^j$ and spatial position $\ell$, with the ordering $\ell = 1$ to $\ell = 4$, respectively associated to the nodes $(j, a, b)$, $(j, a+1, b)$, $(j, a, b+1)$, $(j, a+1, b+1)$, is defined as

$$ \mathbf{S}_\ell(\mathbf{u}_{a,b}^j) = 2\rho_0 \frac{\partial W}{\partial \mathbf{C}}(\mathbf{C}_\ell(\mathbf{u}_{a,b}^j)), \quad (31) $$

for all $\mathbf{u}_{a,b}^j \in \mathcal{X}_d^{\mathbf{u}}$.

### 3.3 2D Barotropic fluid - elastic body interactions

We describe here the discrete variational setting for the coupled dynamics of a barotropic fluid flowing on a hyperelastic body. The Lagrangian density of a barotropic fluid is of the form

$$ \begin{aligned} \mathcal{L}^f(\varphi, \dot{\varphi}, \nabla\varphi) &= L^f(\varphi, \dot{\varphi}, \nabla\varphi)d^2X \wedge dt \\ &= \left[ \frac{1}{2}\rho_0^f \langle \dot{\varphi}, \dot{\varphi} \rangle - \rho_0^f W^f(\rho_0^f, J) - \rho_0^f \Pi(\varphi) \right] d^2X \wedge dt, \end{aligned} \quad (32) $$

with $\rho_0^f$ the fluid mass density and $J = \det \mathbf{F} = \sqrt{\det \mathbf{C}}$, where we assume $\mathbf{G}_{ij} = \delta_{ij}$ and $g_{ab} = \delta_{ab}$. Here $W^f(\rho_0^f, J) = w(\rho_0^f/J)$, with $w(\rho)$ the specific internal energy barotropic fluid. We recall that Hamilton's principle yields the barotropic fluid equations in the Lagrangian (or material) description as

$$ \rho_0^f \ddot{\varphi} = -\frac{\partial}{\partial x^i} (P_W J \mathbf{F}^{-1})^i - \rho_0^f \frac{\partial \Pi}{\partial \varphi}, \quad (33) $$

with $P_W(\rho_0^f, J) = -\rho_0^f \frac{\partial W^f}{\partial J}(\rho_0^f, J)$ the pressure in the material description.
**Discrete Lagrangians for the fluid and the elastic body.** For the elastic body, we consider discrete Lagrangians of the form

$$ \mathcal{L}_d^e(j^1\varphi_d(\mathbb{D})) = \text{vol}(\mathbb{D}) (\rho_0^e K_d(j^1\varphi_d(\mathbb{D})) - \rho_0^e W_d^e(j^1\varphi_d(\mathbb{D})) - \rho_0^e \Pi_d(j^1\varphi_d(\mathbb{D}))) \quad (34) $$

where vol($\mathbb{D}$) is the volume of the parallelepiped $\mathbb{D} \in \mathcal{X}_d^{\mathbb{D}}$. For simplicity, we have assumed that $\rho_0^e$ is constant in space. The discrete kinetic energy $K_d: J^1\mathcal{Y}_d \to \mathbb{R}$ is defined by

$$ K_d(j^1\varphi_d(\mathbb{D}_{a,b}^j)) := \frac{1}{4} \sum_{\alpha=a}^{a+1} \sum_{\beta=b}^{b+1} \frac{1}{2} |\mathbf{v}_{\alpha,\beta}^j|^2 \quad (35) $$

with $\mathbf{v}_{\alpha,\beta}^j = (\varphi_{\alpha,\beta}^{j+1} - \varphi_{\alpha,\beta}^j)/\Delta t^j$. The discrete stored energy function $W_d^e: J^1\mathcal{Y}_d \to \mathbb{R}$ is defined by

$$ W_d^e(j^1\varphi_d(\mathbb{D}_{a,b}^j)) := \frac{1}{4} \sum_{\ell=1}^{4} W^e(\mathbf{C}_\ell(\mathbb{D}_{a,b}^j)), \quad (36) $$

with $W^e$ the stored energy function of the hyperelastic body. Finally, the discrete potential energy $\Pi_d: J^1\mathcal{Y}_d \to \mathbb{R}$ is defined by

$$ \Pi_d(j^1\varphi_d(\mathbb{D}_{a,b}^j)) := \frac{1}{4} \sum_{\alpha=a}^{a+1} \sum_{\beta=b}^{b+1} \Pi(\varphi_{\alpha,\beta}^j), \quad (37) $$

where $\Pi$ is the potential energy of the continuous model. We shall focus on the gravitation potential $\Pi(\varphi) = \langle \mathbf{g}, \varphi \rangle$, with $\mathbf{g}$ the gravitational acceleration vector.

For the barotropic fluid, we shall consider discrete versions of (32) given in a similar way with the elasticity case (34) by

$$ \mathcal{L}_d^f(j^1\varphi_d(\mathbb{D})) = \text{vol}(\mathbb{D}) (\rho_0^f K_d(j^1\varphi_d(\mathbb{D})) - \rho_0^f W_d^f(\rho_0^f, j^1\varphi_d(\mathbb{D})) - \rho_0^f \Pi_d(j^1\varphi_d(\mathbb{D}))) \quad (38) $$

where now $W_d^f: J^1\mathcal{Y}_d \to \mathbb{R}$ is the discrete internal energy of the fluid given by

$$ W_d^f(\rho_0^f, j^1\varphi_d(\mathbb{D}_{c,d}^j)) := \frac{1}{4} \sum_{\ell=1}^{4} W^f(\rho_0^f, J_\ell(\mathbb{D}_{c,d}^j)), \quad (39) $$

with $W^f(\rho_0^f, J)$ the internal energy of the barotropic fluid and $J_\ell(\mathbb{D}_{c,d}^j)$ the discrete Jacobians associated to $\mathbb{D}_{c,d}^j$. We assumed that $\rho_0^f$ is constant in space for simplicity.

**Discrete action and variations.** To simplify the exposition, we assume that the discrete base space configuration for the elastic body is fixed and given by $\phi_{\mathcal{X}_d}(j, a, b) = (j\Delta t, a\Delta s_1, b\Delta s_2)$, for given $\Delta t, \Delta s_1, \Delta s_2$, similarly for the fluid. In this case, we have $\text{vol}(\mathbb{D}) = \Delta t\Delta s_1\Delta s_2$ in the discrete Lagrangian and the mass of each 2D cell in $\phi_{\mathcal{X}_d}(\mathbb{B}_d)$ is $M^e = \rho_0^e\Delta s_1\Delta s_2$.
We consider the discrete action functional

$$
\begin{align}
S_d(\varphi_d) &= \sum_{\mathbb{D} \in \mathcal{X}_d^{\mathbb{D}}} \mathcal{L}_d^{\mathbf{e}}(j^1 \varphi_d(\mathbb{D})) + \mathcal{L}_d^{\mathbf{f}}(j^1 \varphi_d(\mathbb{D})) \nonumber \\
&= \sum_{j=0}^{N-1} \left( \sum_{a=0}^{A-1} \sum_{b=0}^{B-1} \mathcal{L}_d^{\mathbf{e}}(j^1 \varphi_d(\mathbb{D}_{a,b}^j)) + \sum_{c=0}^{C-1} \sum_{d=0}^{D-1} \mathcal{L}_d^{\mathbf{f}}(j^1 \varphi_d(\mathbb{D}_{c,d}^j)) \right), \tag{40}
\end{align}
$$

where $\mathcal{L}_d^e$ and $\mathcal{L}_d^f$ are the discrete Lagrangians associated to the elastic body and to the fluid. For notational simplicity we have used the same notation $\varphi_d$ for the both the discrete configuration of the fluid and the elastic body.

Let us denote by $D_k$, $k = 1, ..., 8$, the partial derivative of the discrete Lagrangians with respect to the $k$-th component of $j^1\varphi_d(\mathbb{D})$, in the order given in (21). From the discrete Hamilton principle $\delta S_d(\varphi_d) = 0$, we get the discrete Euler-Lagrange equations in the general form

$$
\begin{equation}
\begin{split}
& D_1 \mathcal{L}_{a,b}^j + D_2 \mathcal{L}_{a,b}^{j-1} + D_3 \mathcal{L}_{a-1,b}^j + D_4 \mathcal{L}_{a-1,b}^{j-1} \\
& \qquad + D_5 \mathcal{L}_{a,b-1}^j + D_6 \mathcal{L}_{a,b-1}^{j-1} + D_7 \mathcal{L}_{a-1,b-1}^j + D_8 \mathcal{L}_{a-1,b-1}^{j-1} = 0
\end{split}
\tag{41}
\end{equation}
$$

for both $\mathcal{L}_d^e$ and $\mathcal{L}_d^f$, where we used the abbreviate notation $\mathcal{L}_{a,b}^j := \mathcal{L}(j^1\varphi_d(\mathbb{D}_{a,b}^j))$. By using the expressions of the discrete Lagrangians given in (34) and (38), the general form (69) yields

$$
\begin{equation}
\begin{split}
\rho_0^n \frac{\mathbf{v}_{a,b}^{j-1} - \mathbf{v}_{a,b}^j}{\Delta t} - \frac{1}{4} \sum_{\ell=1}^4 \Biggl[ & \mathfrak{D}_1^\ell (\mathbf{F}_\ell(\mathbb{D}_{a,b}^j)\mathbf{S}_\ell(\mathbb{D}_{a,b}^j)) + \mathfrak{D}_3^\ell (\mathbf{F}_\ell(\mathbb{D}_{a-1,b}^j)\mathbf{S}_\ell(\mathbb{D}_{a-1,b}^j)) \\
& + \mathfrak{D}_5^\ell (\mathbf{F}_\ell(\mathbb{D}_{a,b-1}^j)\mathbf{S}_\ell(\mathbb{D}_{a,b-1}^j)) + \mathfrak{D}_7^\ell (\mathbf{F}_\ell(\mathbb{D}_{a-1,b-1}^j)\mathbf{S}_\ell(\mathbb{D}_{a-1,b-1}^j)) \Biggr] \\
& - \rho_0^n \frac{\partial \Pi}{\partial \varphi}(\varphi_{a,b}^j) = 0,
\end{split}
\tag{42}
\end{equation}
$$

$\mathbf{n} = \mathbf{e}, \mathbf{f}$, see Appendix A.1 for the derivation, as well as for the definition of the operators $\mathfrak{D}_k$, $k = 1, 3, 5, 7$. Note that these are the discrete equations of motion corresponding to (8). In (42), $\mathbf{S}_\ell$ is the discrete second Piola-Kirchhoff tensor of the model considered, see (31). For the St. Venant-Kirchhoff model one has $\mathbf{S} = \lambda \mathrm{Tr}(\mathbf{E})\mathbf{l} + 2\mu \mathbf{E}$ while for the barotropic fluid $\mathbf{S} = -P_W \mathbf{J} \mathbf{C}^{-1}$.

Implicit versions of these Euler-Lagrange equations can be derived by considering the appropriate modifications in the definition of the discrete deformation gradients, see [9, §3.2.1]. We also refer to that paper for the discrete boundary conditions emerging from the discrete Hamilton principle.

**Discrete multisymplectic form formula.** By evaluating the differential $dS_d$ of the discrete action functional, restricted to an arbitrary subdomain $U'_d$ of $U_d$, on a solution $\varphi_d$ of the discrete Euler-Lagrange equation and by taking the exterior derivative of the
resulting expression along the first variations of this solution, one obtains a discrete analogue of the multisymplectic form formula (16), see [22] for a detailed treatment. This formula extends to spacetime discretization, the symplectic property of variational time integrators [24]. It also encodes a discrete version of the reciprocity theorem of continuum mechanics, as well as discrete time symplecticity of the solution, see [19]. We refer to [9] for application to the case of the barotropic fluid, the case of the elastic body being similar. In the present situation, the solution of the discrete Euler-Lagrange equations (42) satisfies the discrete multisymplectic form formula on each subdomain $U'_d$ of the discrete spacetime associated to either the fluid motion or the elastic body motion.

**Symmetries and discrete Noether theorems.** Consider the action of the special Euclidean group $SE(2) \ni (R_\theta, u)$ on a discrete field given by $\varphi_d \mapsto R_\theta\varphi_d + u$ with $R_\theta$ the rotation of angle $\theta$. The discrete covariant momentum maps $J_{L_d}^P : J^1Y_d \to \mathfrak{sc}(2)^*$, $\mathbf{p} = 1, \dots, 8$, associated to the $SE(2)$ action and to a discrete Lagrangian $\mathcal{L}_d : J^1Y_d \to \mathbb{R}$, are found as

$$J_{L_d}^P (j^1 \varphi_d (\mathbb{D}_{a,b}^j)) = (\varphi^{(\mathbf{p})} \times D_{\mathbf{p}} \mathcal{L}_{a,b}^j, D_{\mathbf{p}} \mathcal{L}_{a,b}^j) \in \mathfrak{sc}(2)^*. \quad (43)$$

In (43), $\varphi^{(\mathbf{p})}$ is the $\mathbf{p}$-th component of $j^1\varphi_d$, see (21), $D_{\mathbf{p}}$ denotes the partial derivative with respect to the $\mathbf{p}$-th argument of $\mathcal{L}_d$, and $\mathcal{L}_{a,b}^j := \mathcal{L}_d(j^1\varphi_d(\mathbb{D}_{a,b}^j))$. From the discrete covariant momentum maps (43), the discrete classical momentum map is given by

$$\mathbf{J}_d(\varphi^j, \varphi^{j+1}) = \sum_{a=0}^{A-1} \sum_{b=0}^{B-1} (\mathcal{J}_{L_d}^2 + \mathcal{J}_{L_d}^4 + \mathcal{J}_{L_d}^6 + \mathcal{J}_{L_d}^8) (j^1 \varphi_d(\mathbb{D}_{a,b}^j)) \\ = - \sum_{a=0}^{A-1} \sum_{b=0}^{B-1} (\mathcal{J}_{L_d}^1 + \mathcal{J}_{L_d}^3 + \mathcal{J}_{L_d}^5 + \mathcal{J}_{L_d}^7) (j^1 \varphi_d(\mathbb{D}_{a,b}^j)) \in \mathfrak{sc}(2)^*$$

with $\varphi^j$ the collection of all positions at time $t^j$, see [5]. For the discrete Lagrangians (34) and (38), we get the expressions

$$\mathbf{J}_d(\varphi^j, \varphi^{j+1}) = \sum_{a=0}^{A-1} \sum_{b=0}^{B-1} \mathbf{J}^e(j^1 \varphi_d(\mathbb{D}_{a,b}^j)) + \sum_{c=0}^{C-1} \sum_{d=0}^{D-1} \mathbf{J}^f(j^1 \varphi_d(\mathbb{D}_{c,d}^j)) \quad (44)$$

with

$$\mathbf{J}^e(j^1 \varphi_d(\mathbb{D}_{a,b}^j)) = \sum_{\alpha=a}^{a+1} \sum_{\beta=b}^{b+1} \left( \varphi_{\alpha,\beta}^j \times \frac{M^e}{4} v_{\alpha,\beta}^j, \frac{M^e}{4} v_{\alpha,\beta}^j \right) \quad (45)$$

similarly for the fluid.

In absence of gravitational potential, the discrete Lagrangians (34) and (38) are $SE(2)$ invariant. Indeed, when $\varphi_d \mapsto R_\theta\varphi_d + u$ the associated discrete deformation gradients (23) transform as $\mathbf{F}_\ell \mapsto R_\theta\mathbf{F}_\ell$, so that the discrete Cauchy-Green deformation tensors (24) and the discrete Jacobians (27) are $SE(2)$-invariant. The discrete covariant Noether theorem for the discrete covariant momentum maps $J_{L_d}^P$, given as

$$\sum_{\mathbb{D} \in \mathcal{U}'_d} \sum_{\mathbf{p}; \mathbb{D}^{(\mathbf{p})} \in \partial \mathcal{U}'_d} J_{L_d}^{\mathbf{p}}(\mathbb{D}) = 0,$$
thus holds on arbitrary subdomain $U'_d$ of the discrete spacetime associated to either the fluid motion or the elastic body motion. Furthermore, the classical discrete Noether theorem holds for the discrete classical momentum map (73) as

$$ \mathbf{J}_d(\varphi^j, \varphi^{j+1}) = \mathbf{J}_d(\varphi^{j-1}, \varphi^j), $$

for all $t^j$.

When the gravitational potential is taken into account, only the component of the discrete momentum maps associated to horizontal translations satisfies the discrete Noether theorems above.

**Incompressibility and impenetrability conditions.** The constraints will be included in the discrete variational formulation by the addition of appropriate penalty terms.

The first constraint is the incompressibility of the elastic body which is imposed via the condition $J_\ell(\mathbb{D}) - 1 = 0$, for all $\mathbb{D} \in \mathcal{X}_d^{\mathbb{D}}$ and all $\ell = 1, 2, 3, 4$, which mimics the properties of rubber. The associated quadratic penalty function is given by

$$ \Phi_{\text{in}}(j^1 \varphi_d(\mathbb{D})) := \frac{1}{4} \sum_{\ell=1}^{4} \frac{r}{2} (J_{\ell}(\mathbb{D}) - 1)^2, \quad (46) $$

where $r$ is the penalty parameter.

The second constraint is the impenetrability condition between the fluid node $\varphi_{c,d}^j$ flowing above the elastic body nodes $\varphi_{a,b}^j, \varphi_{a+1,b}^j$, which is defined as follows

$$ \Psi_{\text{im}}(\varphi_{a,b}^j, \varphi_{a+1,b}^j, \varphi_{c,d}^j) = \langle (\varphi_{c,d}^j - \varphi_{a+1,b}^j), \mathbf{R}_{\pi/2}(\varphi_{a+1,b}^j - \varphi_{a,b}^j) \rangle \geq 0, \quad (47) $$

where $\mathbf{R}_{\pi/2}$ is a rotation of $+\pi/2$. In this expression $\mathbf{R}_{\pi/2}(\varphi_{a+1,b}^j - \varphi_{a,b}^j)$ represents the outward pointing normal vector to the body. See, e.g., [3] for the use of such constraint functions in impact problems. The function $\Psi_{\text{im}}$ in (47) is evaluated on all the boundary nodes $\varphi_{c,d}^j$ of the fluid and on the two boundary nodes of the elastic body that are the closest to $\varphi_{c,d}^j$. The associated quadratic penalty term is given by

$$ \Phi_{\text{im}}(\varphi_{a,b}^j, \varphi_{a+1,b}^j, \varphi_{c,d}^j) = \frac{1}{2} K \Psi_{\text{im}}(\varphi_{a,b}^j, \varphi_{a+1,b}^j, \varphi_{c,d}^j)^2, \quad (48) $$

with $K \in ]0, \infty[$ if $\Phi_{\text{im}}(\varphi_{a,b}^j, \varphi_{a+1,b}^j, \varphi_{c,d}^j) \leq 0$ and $K = 0$ if $\Phi_{\text{im}}(\varphi_{a,b}^j, \varphi_{a+1,b}^j, \varphi_{c,d}^j) > 0$. The derivatives of $\Phi_{\text{im}}$ with respect to the nodes $\varphi_{a,b}^j, \varphi_{a+1,b}^j, \varphi_{c,d}^j$ are given by

$$ \begin{aligned} D_1 \Phi_{\text{im}}(\varphi_{a,b}^j, \varphi_{a+1,b}^j, \varphi_{c,d}^j) &= -K \Psi_{\text{im}}(\varphi_{a,b}^j, \varphi_{a+1,b}^j, \varphi_{c,d}^j) \mathbf{R}_{\pi/2}^{\top}(\varphi_{c,d}^j - \varphi_{a+1,b}^j) \\ D_2 \Phi_{\text{im}}(\varphi_{a,b}^j, \varphi_{a+1,b}^j, \varphi_{c,d}^j) &= K \Psi_{\text{im}}(\varphi_{a,b}^j, \varphi_{a+1,b}^j, \varphi_{c,d}^j) \mathbf{R}_{\pi/2}(\varphi_{a,b}^j - \varphi_{c,d}^j) \quad (49) \\ D_3 \Phi_{\text{im}}(\varphi_{a,b}^j, \varphi_{a+1,b}^j, \varphi_{c,d}^j) &= K \Psi_{\text{im}}(\varphi_{a,b}^j, \varphi_{a+1,b}^j, \varphi_{c,d}^j) \mathbf{R}_{\pi/2}(\varphi_{a+1,b}^j - \varphi_{a,b}^j), \end{aligned} $$

which express the directions of the reactions forces.
The resulting discrete action functional is found as

$$ \tilde{S}_d(\varphi_d) = S_d(\varphi_d) - \sum_{j=0}^{N-1} \sum_{a=0}^{A-1} \sum_{j=0}^{B-1} \Delta t \Delta s_1 \Delta s_2 \Phi_{\text{in}}(j^1 \varphi_d(\Box^{j}_{a,b})) \\ - \sum_{j=0}^{N-1} \sum_{(c,d) \in \mathbb{B}^f} \Phi_{\text{im}}(\varphi_{a,b}^j, \varphi_{a+1,b}^j, \varphi_{c,d}^j), \qquad (50) $$

where $S_d$ is associated to the discrete Lagrangians of the fluid and the elastic body, see (40). In the third term the second sum is taken over the boundary nodes of the fluid and, given such a boundary node $\varphi_{c,d}^j$, the nodes $\varphi_{a,b}^j$, $\varphi_{a+1,b}^j$ are chosen as explained above. For the boundary nodes $\varphi_{c,d}^j$ that are not in contact, the corresponding term vanishes by definition of $\Phi_{\text{im}}$. Away from the interface between the barotropic fluid and the elastic body, the boundary nodes are subject to the ambient pressure. This boundary condition directly arises from the variations of $S_d$ as explained below. Parts of the boundary of the body can also be prescribed, see below.

In order to avoid the possible instabilities related to the presence of a large penalty parameter, one could use the methods of augmented Lagrangians, in which both a Lagrange multiplier term and a penalty term are used, see [32]. From a practical point of view, see [26], this approach is however much more expensive as it needs to determine the Lagrange multiplier at each contact point.

## 3.4 Barotropic fluid flowing over a St. Venant-Kirchhoff hyperelastic container

We consider an incompressible elastic body with stored energy given by the St. Venant-Kirchhoff model (29) and a compressible barotropic fluid described by the Tait equation

$$ W_d^f(\rho_0, j^1 \varphi_d(\Box)) = \frac{1}{4} \sum_{\ell=1}^{4} \left[ \frac{A}{\gamma-1} \left( \frac{J_{\ell}(\Box)}{\rho_0} \right)^{1-\gamma} + B \left( \frac{J_{\ell}(\Box)}{\rho_0} \right) \right]. \qquad (51) $$

We shall use the expression (51) for the treatment of an isentropic perfect fluid, where the value of the constant $B \neq 0$ does not affect the dynamics, while it allows to naturally impose the boundary condition $P|_{\partial B} = B$, with $P = A (\rho_0 / j)^\gamma$ and $B$ the external pressure. This is crucial for the discretization, since it allows to find the appropriate discretization of the boundary condition directly from the boundary terms of the discrete variational principle, see [9].

The barotropic fluid is subject to gravity and is flowing without friction inside an elastic container. In this case, see Fig. 4, the impenetrability constraints take the following form along the three edges

$$ \begin{aligned} \Psi_{\text{im1}}(\varphi_{a,b}^j, \varphi_{a+1,b}^j, \varphi_{c,d}^j) &= \langle (\varphi_{c,d}^j - \varphi_{a+1,b}^j), \text{R}_{\pi/2}(\varphi_{a+1,b}^j - \varphi_{a,b}^j) \rangle \geq 0 \\ \Psi_{\text{im2}}(\varphi_{a,b}^j, \varphi_{a,b+1}^j, \varphi_{c,d}^j) &= \langle (\varphi_{c,d}^j - \varphi_{a,b}^j), \text{R}_{\pi/2}(\varphi_{a,b}^j - \varphi_{a,b+1}^j) \rangle \geq 0 \\ \Psi_{\text{im3}}(\varphi_{a,b}^j, \varphi_{a,b+1}^j, \varphi_{c,d}^j) &= \langle (\varphi_{c,d}^j - \varphi_{a,b+1}^j), \text{R}_{\pi/2}(\varphi_{a,b+1}^j - \varphi_{a,b}^j) \rangle \geq 0. \end{aligned} \qquad (52) $$
<figure><img src="image_5.png" /><figcaption>Figure 4: Elastic body and fluid configurations respectively in nodes $\varphi_{a,b}^j$, $\varphi_{a+1,b}^j$, $\varphi_{a,b+1}^j$, and $\varphi_{c,d}^j$ at time $t^j$.</figcaption></figure>

The integrator is found by computing the criticality condition for the discrete action functional (50) with incompressibility and impenetrability penalty terms. We simulate the following situation:

1. The barotropic fluid has the properties $\rho_0 = 997 \text{ kg/m}^3$, $\gamma = 6$, $A = \tilde{A}\rho_0^{-\gamma}$ with $\tilde{A} = 3.041 \times 10^4 \text{ Pa}$, and $B = 3.0397 \times 10^4 \text{ Pa}$. The size of the discrete reference configuration at time $t^0$ is $0.5\text{m} \times 0.15\text{m}$, with space-steps $\Delta s_1 = 0.025$, $\Delta s_2 = 0.015\text{m}$.

2. The hyperelastic St. Venant-Kirchhoff model has the properties $\rho_0 = 945 \text{ kg/m}^3$, Poisson ratio $\nu = 0.4999$, and Young modulus $E = 4.5 \times 10^6$. The incompressibility parameter is $r = 10^4$ and the space-steps are $\Delta s_1 = 0.02$, $\Delta s_2 = 0.0125\text{m}$.

3. The two corners at the bottom of the hyperelastic body are fixed.

4. The time-step is $\Delta t = 10^{-4}$ and the test is carried out for $0.4\text{s}$. After this time, the impact between the fluid and the obstacle requires taking into account the physical phenomena of breaking waves.

The results are reproduced in Fig. 5 and Fig. 6, which also show the evolution of the total momentum maps (fluid+solid), the relative total energy behavior (fluid+solid), and the norm of the resultant of forces acting on the right and on the left of the fluid. We used the notation $\mathbf{J}_d = (J_1, J_2, J_3) \in \mathfrak{se}(2)^*$ for the angular, horizontal, and vertical component, respectively. Note that only the component $J_2$ corresponds to a symmetry of the discrete Lagrangian in presence of gravity. Its conservation is however broken, even before the fluid impacts the elastic body, see Fig. 5, due to the fixing of the two corners at the bottom of the hyperelastic body. Its exact conservation in Fig. 6 is due to the symmetric character of the initial configuration. The graph of the resultant force also illustrates the times of impact of the fluid on the right and left parts of the hyperelastic container.

The forces acting on the right at the spacetime nodes $(j, C-1, d)$, $d = 1, \dots, D-1$, of the fluid associated to pressures $P_2(\mathbb{D}_{C-1,d}^j)$ & $P_4(\mathbb{D}_{C-1,d-1}^j)$ are given by

$$ \frac{P_2(\mathbb{D}_{C-1,d}^j)}{4} (\varphi_{C,d}^j - \varphi_{C,d+1}^j) \times \frac{\mathbf{n}_2(\mathbb{D}_{C-1,d}^j)}{|\mathbf{n}_2(\mathbb{D}_{C-1,d}^j)|} + \frac{P_4(\mathbb{D}_{C-1,d-1}^j)}{4} (\varphi_{C,d-1}^j - \varphi_{C,d}^j) \times \frac{\mathbf{n}_4(\mathbb{D}_{C-1,d-1}^j)}{|\mathbf{n}_4(\mathbb{D}_{C-1,d-1}^j)|}, $$

see Appendix A.2 in [9] for details. Similarly, the forces acting on the left at the spacetime nodes $(j, 0, d)$, $d = 1, \dots, D-1$, of the fluid associated to pressures $P_1(\mathbb{D}_{0,d}^j)$ & $P_3(\mathbb{D}_{0,d-1}^j)$
are

$$ \frac{P_1(\mathcal{W}_{0,d}^j)}{4} (\varphi_{0,d+1}^j - \varphi_{0,d}^j) \times \frac{\mathbf{n}_1(\mathcal{W}_{0,d}^j)}{|\mathbf{n}_1(\mathcal{W}_{0,d}^j)|} + \frac{P_3(\mathcal{W}_{0,d-1}^j)}{4} (\varphi_{0,d}^j - \varphi_{0,d-1}^j) \times \frac{\mathbf{n}_3(\mathcal{W}_{0,d-1}^j)}{|\mathbf{n}_3(\mathcal{W}_{0,d-1}^j)|} $$

<figure><img src="image_6.png" /><figcaption>Figure 5: Top to bottom on the left: Configuration after 0.01s, 0.1s, 0.15s, 0.35s, 0.4s. Top to bottom on the right: Evolution of momentum maps, relative energy, and norm of the resultant of forces acting on the right and on the left of the fluid, associated respectively to pressures $P_2(\mathcal{W}_{C-1,d}^j)$ & $P_4(\mathcal{W}_{C-1,d-1}^j)$ and $P_1(\mathcal{W}_{0,d}^j)$ & $P_3(\mathcal{W}_{0,d-1}^j)$ during 0.4s.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 6: *Top to bottom on the left:* Configuration after 0.01s, 0.15s, 0.25s, 0.35s, 0.4s. *Top to bottom on the right:* Evolution of momentum maps, relative energy, and norm of the resultant of forces acting on the right and on the left of the fluid, associated respectively to pressures $P_2(\tilde{\omega}_{C-1,d}^j)$ & $P_4(\tilde{\omega}_{C-1,d-1}^j)$ and $P_1(\tilde{\omega}_{0,d}^j)$ & $P_3(\tilde{\omega}_{0,d-1}^j)$ during 0.5s.</figcaption></figure>

## 3.5 Convergence tests for the hyperelastic model

We consider a St. Venant-Kirchhoff hyperelastic cantilever with mass density $\rho_0 = 945 \text{ kg/m}^3$, Poisson ratio $\nu = 0.4999$, and Young modulus $E = 2.5 \times 10^6$. The model is incompressible, with coefficient $r = 10^4$ in (46). The size of the rectangular reference cantilever configuration at time $t^0$ is $L_1 \times L_2 = 0.8 \text{ m} \times 0.2 \text{ m}$, see Fig. 7.
We consider the explicit integrator described in (42), for the elastic body only, and study the convergence with respect to $\Delta t$ and $\Delta s_i$, $i = 1, 2$.

1. Given a fixed mesh, with values $\Delta s_1 = \Delta s_2 = 0.025\text{m}$, we vary the time-steps as $\Delta t \in \{2 \times 10^{-4}, 10^{-4}, 5 \times 10^{-5}, 2.5 \times 10^{-5}\}$. We compute the $L^2$-errors in the position $\varphi_d$ at time $t^N = 0.5\text{s}$, by comparing $\varphi_d$ with an “exact solution” obtained with the time-step $\Delta t_{\text{ref}} = 6.25 \times 10^{-6}\text{s}$. That is, for each value of $\Delta t$ we calculate

$$ \|\varphi_d - \varphi_{\text{ref}}\|_{L^2} = \left( \sum_a \sum_b \|\varphi_{a,b}^N - \varphi_{\text{ref};a,b}^N\|^2 \right)^{1/2} . \quad (53) $$

This yields the following convergence with respect to $\Delta t$

<table><thead><tr><th>Δt</th><th>2 × 10<sup>-4</sup></th><th>10<sup>-4</sup></th><th>5 × 10<sup>-5</sup></th><th>2.5 × 10<sup>-5</sup></th></tr></thead><tbody><tr><th>||<i>φ</i><sub>d</sub> − <i>φ</i><sub>ref</sub>||<sub>L<sup>2</sup></sub></th><td>1.3 × 10<sup>-3</sup></td><td>6.47 × 10<sup>-4</sup></td><td>3.02 × 10<sup>-4</sup></td><td>1.29 × 10<sup>-4</sup></td></tr><tr><th>rate</th><td></td><td>1.007</td><td>1.1</td><td>1.23</td></tr></tbody></table>

2. Given a fixed time-step $\Delta t = 5 \times 10^{-5}\text{s}$, we vary the space-steps as $\Delta s_1 = \Delta s_2 \in \{0.1, 0.05, 0.025, 0.0125\}$. The “exact solution” is chosen with $\Delta s_{1;\text{ref}} = \Delta s_{2;\text{ref}} = 0.00625\text{m}$. We compute the $L^2$-errors in the position $\varphi_d$ at time $t^N = 0.1\text{s}$. We get the following convergence with respect to $\Delta s_1 = \Delta s_2$.

<table><thead><tr><th>Δs<sub>1</sub> = Δs<sub>2</sub></th><th>0.1</th><th>0.05</th><th>0.025</th><th>0.0125</th></tr></thead><tbody><tr><th>||φ<sub>d</sub> - φ<sub>ref</sub>||<sub>L<sup>2</sup></sub></th><td>0.0481</td><td>0.0151</td><td>0.0044</td><td>0.0012</td></tr><tr><th>rate</th><td></td><td>1.6715</td><td>1.779</td><td>1.8745</td></tr></tbody></table>

The results indicate that the method is second-order accurate with respect to $\Delta s_1 = \Delta s_2$.

<figure><img src="image_8.png" /><figcaption>Figure 7: *From left to right:* Cantilever configuration, relative energy, and momentum maps after 0.2s, when $\Delta t = 5 \times 10^{-5}$ and $\Delta s_1 = \Delta s_2 = 6.25 \times 10^{-3}$.</figcaption></figure>

# 4 3D Discrete fluid - elastic body interactions

In this section we show how the setting developed previously can be extended to the three dimensional case. The main step is to find an appropriate definition of the 3D discrete deformation gradients.
## 4.1 3D discrete setting and the discrete Cauchy-Green tensor

**Discrete Lagrangian setting.** With $\mathcal{B}$ a domain in $\mathbb{R}^3$ with piecewise smooth boundary and $\mathcal{M} = \mathbb{R}^3$, the discrete configuration bundle is defined as in the 2D case, see the beginning of §3.1 and Fig. 1. The discrete parameter space is now $\mathcal{U}_d = \{0, ..., j, ..., N\} \times \mathbb{B}_d$, where we consider $\mathbb{B}_d = \{0, ..., A\} \times \{0, ..., B\} \times \{0, ..., C\}$ and we denote $(j, a, b, c) \in \mathcal{U}_d$. It determines a set $\mathcal{U}_d^{\mathcal{D}}$ of parallelepipeds denoted $\mathcal{D}_{a,b,c}^j$, defined by 16 pairs of indices, see Fig. 8 for the eight pairs of indices in $\mathcal{D}_{a,b,c}^j$ at time $t^j$. We assume that the discrete base-space configuration is of the form

$$ \phi_{\chi_d} : \mathcal{U}_d \ni (j, a, b, c) \mapsto s_{a,b,c}^j = (t^j, s_a^j, s_b^j, s_c^j) \in \chi_d \quad (54) $$

and denote by $\varphi_{a,b,c}^j := \varphi_d(s_{a,b,c}^j)$ the value of the discrete field at $s_{a,b,c}^j$. The discrete first jet bundle and discrete first jet extensions are defined as in the 2D case earlier. In particular, we write the first jet extension as follows

$$ \begin{aligned} j^1 \varphi_d(\mathcal{D}_{a,b,c}^j) = & (\varphi_{a,b,c}^j, \varphi_{a,b,c}^{j+1}, \varphi_{a+1,b,c}^j, \varphi_{a+1,b,c}^{j+1}, \varphi_{a,b+1,c}^j, \varphi_{a,b+1,c}^{j+1}, \varphi_{a,b,c+1}^j, \varphi_{a,b,c+1}^{j+1}, \\ & \varphi_{a+1,b+1,c}^j, \varphi_{a+1,b+1,c}^{j+1}, \varphi_{a,b+1,c+1}^j, \varphi_{a,b+1,c+1}^{j+1}, \varphi_{a+1,b,c+1}^j, \varphi_{a+1,b,c+1}^{j+1}, \varphi_{a+1,b+1,c+1}^j, \varphi_{a+1,b+1,c+1}^{j+1}) \end{aligned} \quad (55) $$

<figure><img src="image_9.png" /><figcaption>Figure 8: Discrete field $\phi_d = \varphi_d \circ \phi_{\chi_d}$ evaluated on $\mathcal{D}_{a,b,c}^j$ at time $t^j$.</figcaption></figure>

**Discrete Cauchy-Green deformation tensor in 3D.** Given a discrete base-space configuration $\phi_{\chi_d}$ and a discrete field $\varphi_d$, the following six vectors $\mathbf{F}_{\ell;a,b,c}^j \in \mathbb{R}^3$, $\ell = 1, ..., 6$ are defined at each node $(j, a, b, c) \in \mathcal{U}_d$, see Fig. 8 on the right:

$$ \begin{aligned} \mathbf{F}_{1;a,b,c}^j &= \frac{\varphi_{a+1,b,c}^j - \varphi_{a,b,c}^j}{|s_{a+1}^j - s_a^j|}, & \mathbf{F}_{2;a,b,c}^j &= \frac{\varphi_{a,b+1,c}^j - \varphi_{a,b,c}^j}{|s_{b+1}^j - s_b^j|}, & \mathbf{F}_{3;a,b,c}^j &= \frac{\varphi_{a,b,c+1}^j - \varphi_{a,b,c}^j}{|s_{c+1}^j - s_c^j|}, \\ \mathbf{F}_{4;a,b,c}^j &= \frac{\varphi_{a-1,b,c}^j - \varphi_{a,b,c}^j}{|s_a^j - s_{a-1}^j|}, & \mathbf{F}_{5;a,b,c}^j &= \frac{\varphi_{a,b-1,c}^j - \varphi_{a,b,c}^j}{|s_b^j - s_{b-1}^j|}, & \mathbf{F}_{6;a,b,c}^j &= \frac{\varphi_{a,b,c-1}^j - \varphi_{a,b,c}^j}{|s_c^j - s_{c-1}^j|}. \end{aligned} $$
**Definition 4.1** The discrete deformation gradients of a discrete field $\varphi_d$ at the parallelepiped $\mathbb{D}_{a,b}^j$ are the eight $3 \times 3$ matrices $\mathbf{F}_\ell(\mathbb{D}_{a,b,c}^j)$, $\ell = 1, \dots, 8$, defined by

$$
\begin{align}
\mathbf{F}_1(\mathbb{D}_{a,b,c}^j) &= [\mathbf{F}_{1;a,b,c}^j \quad \mathbf{F}_{2;a,b,c}^j \quad \mathbf{F}_{3;a,b,c}^j], \nonumber \\
\mathbf{F}_2(\mathbb{D}_{a,b,c}^j) &= [\mathbf{F}_{2;a+1,b,c}^j \quad \mathbf{F}_{4;a+1,b,c}^j \quad \mathbf{F}_{3;a+1,b,c}^j], \nonumber \\
\mathbf{F}_3(\mathbb{D}_{a,b,c}^j) &= [\mathbf{F}_{5;a,b+1,c}^j \quad \mathbf{F}_{1;a,b+1,c}^j \quad \mathbf{F}_{3;a,b+1,c}^j], \nonumber \\
\mathbf{F}_4(\mathbb{D}_{a,b,c}^j) &= [\mathbf{F}_{2;a,b,c+1}^j \quad \mathbf{F}_{1;a,b,c+1}^j \quad \mathbf{F}_{6;a,b,c+1}^j], \tag{56} \\
\mathbf{F}_5(\mathbb{D}_{a,b,c}^j) &= [\mathbf{F}_{4;a+1,b+1,c}^j \quad \mathbf{F}_{5;a+1,b+1,c}^j \quad \mathbf{F}_{3;a+1,b+1,c}^j], \nonumber \\
\mathbf{F}_6(\mathbb{D}_{a,b,c}^j) &= [\mathbf{F}_{1;a,b+1,c+1}^j \quad \mathbf{F}_{5;a,b+1,c+1}^j \quad \mathbf{F}_{6;a,b+1,c+1}^j], \nonumber \\
\mathbf{F}_7(\mathbb{D}_{a,b,c}^j) &= [\mathbf{F}_{4;a+1,b,c+1}^j \quad \mathbf{F}_{2;a+1,b,c+1}^j \quad \mathbf{F}_{6;a+1,b,c+1}^j], \nonumber \\
\mathbf{F}_8(\mathbb{D}_{a,b,c}^j) &= [\mathbf{F}_{5;a+1,b+1,c+1}^j \quad \mathbf{F}_{4;a+1,b+1,c+1}^j \quad \mathbf{F}_{6;a+1,b,c+1}^j]. \nonumber
\end{align}
$$

We note that the discrete deformation gradients are defined at each of the eight nodes of $\mathbb{D}_{a,b,c}^j$ that are associated to time $t^j$. We assume that the determinant of these matrices is always positive. The ordering $\ell = 1$ to $\ell = 8$ is respectively associated to the nodes $(a, b, c)$, $(a+1, b, c)$, $(a, b+1, c)$, $(a, b, c+1)$, $(a+1, b+1, c)$, $(a, b+1, c+1)$, $(a+1, b, c+1)$, $(a+1, b+1, c+1)$ see Fig. 8 on the left. Note that with this concept of discrete deformation gradient, one can consider the right and left polar decompositions $\mathbf{F} = \mathbf{RU} = \mathbf{VR}$ exactly as in the continuous case.

**Definition 4.2** The 3D discrete Cauchy-Green deformation tensors of a discrete field $\varphi_d$ at the parallelepiped $\mathbb{D}_{a,b,c}^j$ are the eight $3 \times 3$ symmetric matrices

$$ \mathbf{C}_\ell(\mathbb{D}_{a,b,c}^j) = \mathbf{F}_\ell(\mathbb{D}_{a,b,c}^j)^\mathrm{T} \mathbf{F}_\ell(\mathbb{D}_{a,b,c}^j), \quad \ell = 1, \dots, 8. \qquad (57) $$

As in Proposition 3.3, we have following result.

**Proposition 4.3** Given a discrete field $\varphi_d: \mathcal{X}_d \to \mathcal{M}$, the 3D discrete Cauchy-Green deformation tensor associated to $\mathbb{D}_{a,b,c}^j$ at the node $(a, b, c)$ at time $t^j$, for $\ell = 1$, is

$$
\mathbf{C}_1(\mathbb{D}_{a,b,c}^j) = \begin{bmatrix} \langle \mathbf{F}_{1;a,b,c}^j, \mathbf{F}_{1;a,b,c}^j \rangle & \langle \mathbf{F}_{1;a,b,c}^j, \mathbf{F}_{2;a,b,c}^j \rangle & \langle \mathbf{F}_{1;a,b,c}^j, \mathbf{F}_{3;a,b,c}^j \rangle \\ \langle \mathbf{F}_{2;a,b,c}^j, \mathbf{F}_{1;a,b,c}^j \rangle & \langle \mathbf{F}_{2;a,b,c}^j, \mathbf{F}_{2;a,b,c}^j \rangle & \langle \mathbf{F}_{2;a,b,c}^j, \mathbf{F}_{3;a,b,c}^j \rangle \\ \langle \mathbf{F}_{3;a,b,c}^j, \mathbf{F}_{1;a,b,c}^j \rangle & \langle \mathbf{F}_{3;a,b,c}^j, \mathbf{F}_{2;a,b,c}^j \rangle & \langle \mathbf{F}_{3;a,b,c}^j, \mathbf{F}_{3;a,b,c}^j \rangle \end{bmatrix}. \quad (58)
$$

See Appendix A.3 for the discrete Cauchy-Green deformation tensor at the nodes $(a+1, b, c)$, $(a, b+1, c)$, $(a, b, c+1)$, $(a+1, b+1, c)$, $(a, b+1, c+1)$, $(a+1, b, c+1)$, $(a+1, b+1, c+1)$.
**Discrete Jacobian.** The 3D discrete Jacobian of $\varphi_d$ associated to $\mathbb{D}_{a,b,c}^j$ at the node $(j, a, b, c)$ is defined by

$$J_1(\mathbb{D}_{a,b,c}^j) = (\mathbf{F}_{1;a,b,c}^j \times \mathbf{F}_{2;a,b,c}^j) \cdot \mathbf{F}_{3;a,b,c}^j \qquad (59)$$

We refer to Appendix A.2 for the expression of the 3D discrete Jacobians at the other nodes. The relations with the 3D discrete deformation gradients and discrete Cauchy-Green tensors hold exactly as in the 2D case in Proposition 3.4.

## 4.2 Discrete frame indifferent and isotropic nonlinear hyperelastic models

**Invariants of the discrete Cauchy-Green deformation tensor.** Our definition of the discrete Cauchy-Green deformation tensor as a symmetric matrix allows the definition of its invariants exactly as in the continuous case. Note that in the discrete case each node has several associated discrete Cauchy-Green deformation tensors, namely, one for each parallelepiped sharing this node. This follows directly from the discrete field theoretic Lagrangian setting that we follow. We have the following expressions.

**Definition 4.4** The invariants of the symmetric matrix $C_\ell(\mathbb{D}_{a,b,c}^j)$, $\ell = 1, \dots, 8$, are defined by

$$ \begin{aligned} I_1(\mathbf{C}_\ell(\mathbb{D}_{a,b,c}^j)) &= \mathrm{Tr}(\mathbf{C}_\ell(\mathbb{D}_{a,b,c}^j)) \\ I_2(\mathbf{C}_\ell(\mathbb{D}_{a,b,c}^j)) &= \frac{1}{2} \left[ I_1(\mathbf{C}_\ell(\mathbb{D}_{a,b,c}^j))^2 - \mathrm{Tr}(\mathbf{C}_\ell(\mathbb{D}_{a,b,c}^j)^2) \right] \\ &= \det(\mathbf{C}_\ell(\mathbb{D}_{a,b,c}^j)) \mathrm{Tr}(\mathbf{C}_\ell(\mathbb{D}_{a,b,c}^j)^{-1}) \\ I_3(\mathbf{C}_\ell(\mathbb{D}_{a,b,c}^j)) &= \det(\mathbf{C}_\ell(\mathbb{D}_{a,b,c}^j)). \end{aligned} \qquad (60) $$

For each fixed $\mathbb{D}_{a,b,c}^j$, these invariants depend only of the first jet extension $j^1\varphi_d(\mathbb{D}_{a,b,c}^j)$ of the discrete field at the given $\mathbb{D}_{a,b,c}^j$. Recall that these invariants are related to the coefficients in the characteristic polynomial $P(\nu)$ of the $3 \times 3$ symmetric matrix $C_\ell(\mathbb{D}_{a,b,c}^j)$ as follows

$$P(\nu) = \nu^3 - I_1(\mathbf{C}_\ell(\mathbb{D}_{a,b,c}^j))\nu^2 + I_2(\mathbf{C}_\ell(\mathbb{D}_{a,b,c}^j))\nu - I_3(\mathbf{C}_\ell(\mathbb{D}_{a,b,c}^j)).$$

In particular, the eigenvalues $\nu_1, \nu_2, \nu_3$ are defined at each node with respect to each of the eight neighbouring parallelepipeds. For each $\mathbb{D}_{a,b,c}^j \in \mathcal{X}_d^j$ and $\ell = 1, \dots, 8$, we have

$$ \begin{aligned} I_1(\mathbf{C}_\ell(\mathbb{D}_{a,b,c}^j)) &= \nu_1 + \nu_2 + \nu_3, & I_2(\mathbf{C}_\ell(\mathbb{D}_{a,b,c}^j)) &= \nu_1\nu_2 + \nu_1\nu_3 + \nu_2\nu_3, \\ I_3(\mathbf{C}_\ell(\mathbb{D}_{a,b,c}^j)) &= \nu_1\nu_2\nu_3, \end{aligned} $$

and $\nu_1 = \lambda_1^2, \nu_2 = \lambda_2^2, \nu_3 = \lambda_3^2$, with $\lambda_1, \lambda_2, \lambda_3$ the discrete principal stretches. At each node, there are thus several classes of discrete principal stretches, namely, one for each parallelepiped sharing this nodes. We can also recall the following standard result, which holds here in the discrete case and can be regarded as a discrete analogue of the characterization of the stored energy function for frame indifferent isotropic hyperelastic materials.
**Proposition 4.5** The following are equivalent:

(i) A scalar function f of $\mathbf{C}_\ell(\mathbb{D}_{a,b,c}^j)$ is invariant under orthogonal transformations.

(ii) f is a function of the invariants of $\mathbf{C}_\ell(\mathbb{D}_{a,b,c}^j)$.

(iii) f is a symmetric function of the principal stretches.

**3D Discrete stored energy function.** Given the stored energy function $W(I_1, I_2, I_3)$ of a frame indifferent isotropic hyperelastic material, the associated discrete stored energy function is defined as $W_d: J^1\mathcal{Y}_d \to \mathbb{R}$ with

$$W_d(j^1\varphi_d(\mathbb{D})) = \frac{1}{8} \sum_{\ell=1}^{8} W(I_1(\mathbf{C}_\ell(\mathbb{D})), I_2(\mathbf{C}_\ell(\mathbb{D})), I_3(\mathbf{C}_\ell(\mathbb{D}))). \quad (61)$$

For each $\mathbb{D} \in \mathcal{X}_d^{\mathbb{D}}$ the value of $W_d$ is found as the averaged value of the stored energy function evaluated on the discrete Cauchy-Green tensor at the nodes of $\mathbb{D}$ at time $t^j$.

**3D Discrete second Piola-Kirchhoff stress.** It is defined at each node of each $\mathbb{D} \in \mathcal{X}_d^{\mathbb{D}}$ as

$$\mathbf{S}_\ell(\mathbb{D}) = 2\rho_0 \frac{\partial W}{\partial \mathbf{C}}(\mathbf{C}_\ell(\mathbb{D})) = \mathbf{S}(\mathbf{C}_\ell(\mathbb{D})), \quad (62)$$

$\ell = 1, ..., 8$. In particular, for an isotropic material, it takes the form

$$\mathbf{S}_\ell(\mathbb{D}) = 2\rho_0 \left[ \frac{\partial W}{\partial I_1} \mathbf{G}^{-1} + \left( \frac{\partial W}{\partial I_2} I_2 + \frac{\partial W}{\partial I_3} I_3 \right) (\mathbf{C}_\ell(\mathbb{D}))^{-1} - \frac{\partial W}{\partial I_2} I_3 (\mathbf{C}_\ell(\mathbb{D}))^{-2} \right]. \quad (63)$$

**Example 4.6 (Nonlinear Mooney-Rivlin incompressible model)** For this model, see (14), the discrete stored energy becomes

$$W_d(j^1\varphi_d(\mathbb{D}_{a,b,c}^j)) = \frac{1}{8} \sum_{\ell=1}^{8} \left[ C_1(I_1(\mathbf{C}_\ell(\mathbb{D}_{a,b,c}^j)) - 3) + C_2(I_2(\mathbf{C}_\ell(\mathbb{D}_{a,b,c}^j)) - 3) \right] \quad (64)$$

and the 3D discrete second Piola-Kirchhoff stress is

$$\mathbf{S}_\ell(\mathbb{D}) = 2\rho_0 \left( C_1 \mathbf{G}^{-1} + C_2 I_2(\mathbf{C}_\ell(\mathbb{D})) \mathbf{C}_\ell(\mathbb{D})^{-1} - C_2 I_3(\mathbf{C}_\ell(\mathbb{D})) \mathbf{C}_\ell(\mathbb{D})^{-2} \right).$$

Incompressibility is imposed via a Lagrange multiplier exactly as in the 2D case earlier.
## 4.3 3D Barotropic fluid - elastic body interactions

**Discrete Lagrangians for the fluid and the elastic body.** The discrete Lagrangian $\mathcal{L}_d^e$ of the elastic body is defined exactly as in the 2D case in (34) with the obvious modifications, namely we take the discrete kinetic energy

$$ K_d(j^1 \varphi_d(\mathbb{D}_{a,b,c}^j)) := \frac{1}{8} \sum_{\alpha=a}^{a+1} \sum_{\beta=b}^{b+1} \sum_{\gamma=c}^{c+1} \frac{1}{2} |\mathbf{v}_{\alpha,\beta,\gamma}^j|^2, \quad (65) $$

where $\mathbf{v}_{\alpha,\beta,\gamma}^j = (\varphi_{a,b,c}^{j+1} - \varphi_{a,b,c}^j)/\Delta t$, we take the discrete stored energy function

$$ W_d^e(j^1 \varphi_d(\mathbb{D}_{a,b,c}^j)) := \frac{1}{8} \sum_{\ell=1}^{8} W^e(\mathbf{C}_\ell(\mathbb{D}_{a,b,c}^j)) \quad (66) $$

and the gravitational potential energy density

$$ \Pi_d(\mathbb{D}_{a,b,c}^j) := \frac{1}{8} \sum_{\alpha=a}^{a+1} \sum_{\beta=b}^{b+1} \sum_{\gamma=c}^{c+1} \Pi(\varphi_{\alpha,\beta,\gamma}^j). \quad (67) $$

We proceed similarly for the discrete Lagrangian $\mathcal{L}_d^f$ of the fluid, see (38).

**Discrete action and discrete Euler-Lagrange equations.** As earlier, we assume $\phi_{\mathcal{X}_d}(j, a, b, c) = (j\Delta t, a\Delta s_1, b\Delta s_2, c\Delta s_3)$ and the mass of each 3D cell in $\Phi_{\mathcal{X}_d}(\mathbb{B}_d)$ is $M^e = \rho_0^e \Delta s_1 \Delta s_2 \Delta s_3$. The discrete action functional reads

$$ S_d(\varphi_d) = \sum_{j=0}^{N-1} \left( \sum_{a=0}^{A-1} \sum_{b=0}^{B-1} \sum_{c=0}^{C-1} \mathcal{L}_d^e(j^1 \varphi_d(\mathbb{D}_{a,b,c}^j)) + \sum_{d=0}^{D-1} \sum_{e=0}^{E-1} \sum_{f=0}^{F-1} \mathcal{L}_d^f(j^1 \varphi_d(\mathbb{D}_{d,e,f}^j)) \right). \quad (68) $$

Let us denote by $D_k, k = 1, \dots, 16$, the partial derivative of the discrete Lagrangians with respect to the k-th component of $j^1 \varphi_d(\mathbb{D})$, in the order given in (55). From the discrete Hamilton principle $\delta S_d(\varphi_d) = 0$, we get the discrete Euler-Lagrange equations in the general form

$$ \begin{aligned} & D_1 \mathcal{L}_{a,b,c}^j + D_2 \mathcal{L}_{a,b,c}^{j-1} + D_3 \mathcal{L}_{a-1,b,c}^j + D_4 \mathcal{L}_{a-1,b,c}^{j-1} + D_5 \mathcal{L}_{a,b-1,c}^j + D_6 \mathcal{L}_{a,b-1,c}^{j-1} + D_7 \mathcal{L}_{a,b,c-1}^j \\ & \qquad + D_8 \mathcal{L}_{a,b,c-1}^{j-1} + D_9 \mathcal{L}_{a-1,b-1,c}^j + D_{10} \mathcal{L}_{a-1,b-1,c}^{j-1} + D_{11} \mathcal{L}_{a,b-1,c-1}^j + D_{12} \mathcal{L}_{a,b-1,c-1}^{j-1} \\ & \qquad + D_{13} \mathcal{L}_{a-1,b,c-1}^j + D_{14} \mathcal{L}_{a-1,b,c-1}^{j-1} + D_{15} \mathcal{L}_{a-1,b-1,c-1}^j + D_{16} \mathcal{L}_{a-1,b-1,c-1}^{j-1} = 0 \end{aligned} \quad (69) $$

for both $\mathcal{L}_d^e$ and $\mathcal{L}_d^f$, where we used the abbreviate notation $\mathcal{L}_{a,b,c}^j := \mathcal{L}(j^1 \varphi_d(\mathbb{D}_{a,b,c}^j))$. For the discrete Lagrangians based on the discretizations given in (65)-(67), the discrete Euler-Lagrange equations (69) yield

$$ \begin{aligned} & \rho_0^n \frac{\mathbf{v}_{a,b,c}^{j-1} - \mathbf{v}_{a,b,c}^j}{\Delta t} - \frac{1}{8} \sum_{\ell=1}^{4} \left[ \mathfrak{D}_1^\ell (\mathbf{F}_\ell(\mathbb{D}_{a,b,c}^j)\mathbf{S}_\ell(\mathbb{D}_{a,b,c}^j)) + \mathfrak{D}_3^\ell (\mathbf{F}_\ell(\mathbb{D}_{a-1,b,c}^j)\mathbf{S}_\ell(\mathbb{D}_{a-1,b,c}^j)) \right. \\ & \qquad \left. + \mathfrak{D}_5^\ell (\mathbf{F}_\ell(\mathbb{D}_{a,b-1,c}^j)\mathbf{S}_\ell(\mathbb{D}_{a,b-1,c}^j)) + \dots + \mathfrak{D}_{15}^\ell (\mathbf{F}_\ell(\mathbb{D}_{a-1,b-1,c-1}^j)\mathbf{S}_\ell(\mathbb{D}_{a-1,b-1,c-1}^j)) \right] \\ & \qquad - \rho_0^n \frac{\partial \Pi}{\partial \varphi}(\varphi_{a,b}^j) = 0, \end{aligned} \quad (70) $$
n = e, f. The derivation of (70) is similar to the 2D case, see Appendix A.1. These equations give the dynamics at the interior of the elastic body and the fluid.

**Impenetrability conditions.** The impenetrability conditions between the fluid nodes $\varphi_{d,e,f}^j, \varphi_{d+1,e,f}^j, \varphi_{d,e+1,f}^j, \varphi_{d+1,e+1,f}^j$ and the elastic body nodes are defined as follows, see Fig. 9,

$$
\begin{align}
\Psi_{\text{im1}}(\varphi_{a,b,c}^j, \varphi_{a-1,b,c}^j, \varphi_{a,b-1,c}^j, \varphi_{d,e,f}^j) & = \langle (\varphi_{d,e,f}^j - \varphi_{a,b,c}^j), (\varphi_{a-1,b,c}^j - \varphi_{a,b,c}^j) \times (\varphi_{a,b-1,c}^j - \varphi_{a,b,c}^j) \rangle \ge 0 \nonumber \\
\Psi_{\text{im2}}(\varphi_{a,b,c}^j, \varphi_{a+1,b,c}^j, \varphi_{a,b-1,c}^j, \varphi_{d+1,e,f}^j) & = \langle (\varphi_{d+1,e,f}^j - \varphi_{a,b,c}^j), (\varphi_{a,b-1,c}^j - \varphi_{a,b,c}^j) \times (\varphi_{a+1,b,c}^j - \varphi_{a,b,c}^j) \rangle \ge 0 \tag{71} \\
\Psi_{\text{im3}}(\varphi_{a,b,c}^j, \varphi_{a-1,b,c}^j, \varphi_{a,b+1,c}^j, \varphi_{d,e+1,f}^j) & = \langle (\varphi_{d,e+1,f}^j - \varphi_{a,b,c}^j), (\varphi_{a,b+1,c}^j - \varphi_{a,b,c}^j) \times (\varphi_{a-1,b,c}^j - \varphi_{a,b,c}^j) \rangle \ge 0 \nonumber \\
\Psi_{\text{im4}}(\varphi_{a,b,c}^j, \varphi_{a+1,b,c}^j, \varphi_{a,b+1,c}^j, \varphi_{d+1,e+1,f}^j) & = \langle (\varphi_{d+1,e+1,f}^j - \varphi_{a,b,c}^j), (\varphi_{a+1,b,c}^j - \varphi_{a,b,c}^j) \times (\varphi_{a,b+1,c}^j - \varphi_{a,b,c}^j) \rangle \ge 0. \nonumber
\end{align}
$$

The expression $(\varphi_{a-1,b,c}^j - \varphi_{a,b,c}^j) \times (\varphi_{a,b-1,c}^j - \varphi_{a,b,c}^j)$ represents the outward pointing normal vector field to the body, similarly in the other three conditions.

<figure><img src="image_10.png" /><figcaption>Figure 9: Fluid flowing above the hyperelastic body.</figcaption></figure>

The associated penalty terms $\Phi_{\text{im}i}$, $i = 1, ..., 4$ are defined as in (48). The directions of the reaction forces are given by the derivative with respect to the positions of the nodes, i.e.,

$$
\begin{align*}
D_1 \Psi_{\text{im1}}(\varphi_{a,b,c}^j, \varphi_{a-1,b,c}^j, \varphi_{a,b-1,c}^j, \varphi_{d,e,f}^j), \quad D_2 \Psi_{\text{im1}}(\varphi_{a,b,c}^j, \varphi_{a-1,b,c}^j, \varphi_{a,b-1,c}^j, \varphi_{d,e,f}^j), \\
D_3 \Psi_{\text{im1}}(\varphi_{a,b,c}^j, \varphi_{a-1,b,c}^j, \varphi_{a,b-1,c}^j, \varphi_{d,e,f}^j), \quad D_4 \Psi_{\text{im1}}(\varphi_{a,b,c}^j, \varphi_{a-1,b,c}^j, \varphi_{a,b-1,c}^j, \varphi_{d,e,f}^j),
\end{align*}
$$

similarly for $\Psi_{\text{im2}}$, $\Psi_{\text{im3}}$, $\Psi_{\text{im4}}$, see Appendix A.4 for the expressions of these derivatives.
The resulting discrete action functional is

$$
\begin{align}
\tilde{S}_d(\varphi_d) &= S(\varphi_d) - \sum_{j=0}^{N-1} \sum_{a=0}^{A-1} \sum_{b=0}^{B-1} \sum_{c=0}^{C-1} \Delta t \Delta s_1 \Delta s_2 \Delta s_3 \Phi_{\text{in}}(j^1 \varphi_d(\mathbb{D}_{a,b,c}^j)) \nonumber \\
&\quad - \sum_{j=0}^{N-1} \sum_{(d,e,f) \in \mathbb{B}^f} \left[ \begin{aligned}
&\Phi_{\text{im}_1}(\varphi_{a,b,c}^j, \varphi_{a-1,b,c}^j, \varphi_{a,b-1,c}^j, \varphi_{d,e,f}^j) + \Phi_{\text{im}_2}(\varphi_{a,b,c}^j, \varphi_{a-1,b,c}^j, \varphi_{a,b-1,c}^j, \varphi_{d,e,f}^j) \\
&+ \Phi_{\text{im}_3}(\varphi_{a,b-1,c}^j, \varphi_{a-1,b-1,c}^j, \varphi_{a,b,c}^j, \varphi_{d,e,f}^j) + \Phi_{\text{im}_4}(\varphi_{a-1,b-1,c}^j, \varphi_{a,b-1,c}^j, \varphi_{a-1,b,c}^j, \varphi_{d,e,f}^j)
\end{aligned} \right], \tag{72}
\end{align}
$$

where $S_d$ is given in (68). The second term in (72) is associated to the incompressibility of the elastic body. In the third term, the second sum is taken over the boundary nodes of the fluid and, given such a boundary node $\varphi_{d,e,f}^j$, the nodes of the elastic body in the four terms $i = 1, \dots, 4$ are chosen in accordance with (71) for each functions $\Phi_{\text{im}_i}$, $i = 1, \dots, 4$.

**Symmetries and discrete Noether theorems.** The discrete classical momentum map associated to the action of the special Euclidean group $SE(3)$ takes the form

$$
\mathbf{J}_d(\varphi^j, \varphi^{j+1}) = \sum_{a=0}^{A-1} \sum_{b=0}^{B-1} \sum_{c=0}^{C-1} \mathbf{J}^e(j^1 \varphi_d(\mathbb{D}_{a,b,c}^j)) + \sum_{d=0}^{D-1} \sum_{e=0}^{E-1} \sum_{f=0}^{F-1} \mathbf{J}^f(j^1 \varphi_d(\mathbb{D}_{d,e,f}^j)) \quad (73)
$$

with

$$
\mathbf{J}^e(j^1 \varphi_d(\mathbb{D}_{a,b,c}^j)) = \sum_{\alpha=a}^{a+1} \sum_{\beta=b}^{b+1} \sum_{\gamma=c}^{c+1} \left( \varphi_{\alpha,\beta}^j \times \frac{M^e}{8} v_{\alpha,\beta,\gamma}^j, \frac{M^e}{8} v_{\alpha,\beta,\gamma}^j \right) \quad (74)
$$

similarly for the fluid.

4.4 Barotropic fluid flowing over a Mooney-Rivlin hyperelastic support

We consider an incompressible Mooney-Rivlin hyperelastic body, see (64), and a barotropic
fluid described by the Tait equation with discrete internal energy term

$$
W_d^f(\rho_0, j^1 \varphi_d(\mathbb{D})) = \frac{1}{8} \sum_{\ell=1}^{8} \left[ \frac{A}{\gamma-1} \left( \frac{J_\ell(\mathbb{D})}{\rho_0} \right)^{1-\gamma} + B \left( \frac{J_\ell(\mathbb{D})}{\rho_0} \right) \right]. \quad (75)
$$

The integrator is found by computing the criticality condition for the discrete action
functional which includes the incompressibility and impenetrability penalty terms. We
simulate the following situation, see Fig. 10:

1. The barotropic fluid has the properties $\rho_0 = 997 \text{ kg/m}^2$, $\gamma = 6$, $A = \tilde{A}\rho_0^{-\gamma}$ with $\tilde{A} = 3.041 \times 10^4 \text{ Pa}$, and $B = 3.0397 \times 10^4 \text{ Pa}$. The size of the discrete reference configuration at time $t^0$ is 0.3m×0.3m×0.05m, with space-steps $\Delta s_1 = \Delta s_2 = 0.05$, $\Delta s_3 = 0.025 \text{m}$.
2. The Mooney-Rivlin model has the properties $\rho_0 = 945 \text{kg/m}^3$, $C_1 = 1.848$, and $C_2 = 0.264$ such that $C_1/C_2 = 7$. The model is nearly incompressible, with coefficient $r = 10^4$ in the corresponding penalty term. The size of the discrete reference configuration at time $t^0$ is $1\text{m} \times 1\text{m} \times 0.3\text{m}$ with space-steps $\Delta s_1 = \Delta s_2 = 0.0625\text{m}$, $\Delta s_3 = 0.1\text{m}$.

3. The four upper edges of the Mooney-Rivlin incompressible model are fixed.

4. The time-step is $\Delta t = 2 \times 10^{-4}$ and the test is carried out for $0.4s$.

The results are reproduced in Fig. 10, where we used the notation $\mathbf{J}_d = (J_1, J_2, J_3, J_4, J_5, J_6) \in \mathfrak{S}\mathfrak{C}(3)^*$. Due to the presence of the gravitation potential, the discrete Lagrangian is only invariant under the subgroup $\mathcal{SE}(2)$ of Euclidean transformations of the horizontal plane. This is why only the components $J_3, J_4, J_5$ of the discrete momentum map $\mathbf{J}_d$ are almost preserved. As in our 2D test, the fixing of the four upper edges prevents the exact conservation of these momentum maps.
<figure><img src="image_11.png" /><figcaption>Figure 10: *Top to bottom on the left:* Configuration after 0.002s, 0.1s, 0.4s, 0.5s, 0.6s. *Top to bottom on the right:* Evolution of total (fluid+solid) momentum maps, total (fluid+solid) energy, and relative total energy during 0.5s.</figcaption></figure>
## 4.5 Convergence tests for the isotropic nonlinear hyperelastic model

Consider the Mooney-Rivlin incompressible model described above, with properties $\rho_0 = 945 \text{ kg/m}^3$, $C_1 = 1.848$, and $C_2 = 0.264$ with $C_1/C_2 = 7$. The model is incompressible, with coefficient $r = 10^4$ in the penalty term. The size of the discrete reference configuration at time $t^0$ is $0.4 \times 0.4 \times 0.2 \text{ m}$. As earlier, the four upper edges of the Mooney-Rivlin incompressible model are fixed, see Fig. 11. The motion is only induced by gravity forces.

We consider the explicit integrator described in (70), for the elastic body only, and study the convergence with respect to $\Delta t$ and $\Delta s_i$, $i = 1, 2, 3$.

1. Given a fixed mesh, with values $\Delta s_1 = \Delta s_2 = \Delta s_3 = 0.1 \text{ m}$, we vary the time-steps as $\Delta t \in \{2 \times 10^{-4}, 10^{-4}, 5 \times 10^{-5}, 2.5 \times 10^{-5}\}$. We compute the $L^2$-errors in the position $\varphi_d$ at time $t^N = 0.1 \text{ s}$, by comparing $\varphi_d$ with an “exact solution” obtained with the time-step $\Delta t_{\text{ref}} = 6.25 \times 10^{-6} \text{ s}$. That is, for each value of $\Delta t$ we calculate

$$ \|\varphi_d - \varphi_{\text{ref}}\|_{L^2} = \left( \sum_a \sum_b \sum_c \|\varphi_{a,b,c}^N - \varphi_{\text{ref};a,b,c}^N\|^2 \right)^{1/2} . \quad (76) $$

This yields the following convergence with respect to $\Delta t$

<table><thead><tr><th>Δt</th><th>2 × 10<sup>-4</sup></th><th>10<sup>-4</sup></th><th>5 × 10<sup>-5</sup></th><th>2.5 × 10<sup>-5</sup></th></tr></thead><tbody><tr><th>||<i>φ</i><sub>d</sub> − <i>φ</i><sub>ref</sub>||<sub>L<sup>2</sup></sub></th><td>5.23 × 10<sup>-4</sup></td><td>2.55 × 10<sup>-4</sup></td><td>1.18 × 10<sup>-4</sup></td><td>5.035 × 10<sup>-5</sup></td></tr><tr><th>rate</th><td></td><td>1.036</td><td>1.12</td><td>1.29</td></tr></tbody></table>

2. Given a fixed time-step $\Delta t = 2 \times 10^{-5} \text{ s}$, we vary the space-steps as $\Delta s_1 = \Delta s_2 = \Delta s_3 \in \{0.1, 0.05, 0.025, 0.0125\}$. The “exact solution” is chosen with $\Delta s_{1;\text{ref}} = \Delta s_{2;\text{ref}} = \Delta s_{3;\text{ref}} = 0.00625 \text{ m}$. We compute the $L^2$-errors in the position $\varphi_d$ at time $t^N = 0.05 \text{ s}$. We get the following convergence with respect to $\Delta s_1 = \Delta s_2 = \Delta s_3$.

<table><tr><th>Δs<sub>1</sub> = Δs<sub>2</sub> = Δs<sub>3</sub></th><td>0.1</td><td>0.05</td><td>0.025</td><td>0.0125</td></tr><tr><th>||<i>φ</i><sub>d</sub> − <i>φ</i><sub>ref</sub>||<sub>L<sup>2</sup></sub></th><td>0.0375</td><td>0.0284</td><td>0.0177</td><td>0.008</td></tr><tr><th>rate</th><td></td><td>0.40</td><td>0.682</td><td>1.1457</td></tr></table>

<figure><img src="image_12.png" /><figcaption>Figure 11: From left to right: Configuration, relative energy, and momentum maps after 0.05s, when $\Delta t = 2 \times 10^{-5}$ and $\Delta s_1 = \Delta s_2 = \Delta s_3 = 1.25 \times 10^{-2}$.</figcaption></figure>
In Fig. 11 on the right, the components $J_3, J_4, J_5$ of the discrete momentum map are preserved, which corresponds to the symmetry subgroup $SE(2)$.

# 5 Concluding remarks and future directions

In this paper we introduced a structure preserving numerical integrator for hyperelastic solids, which is spacetime multisymplectic, symplectic in time, preserves exactly the momenta associated to symmetries, and nearly preserves total energy. The integrator is derived from a discrete variational formulation of hyperelasticity in the Lagrangian description. Thanks to its variational nature and its adequacy with a previously derived variational discretization for fluids, the scheme can naturally be extended to accommodate fluid-structure interaction problems via the inclusion of impenetrability penalty terms in the discrete action functional. Important steps in our approach were the definition of discrete deformation gradients, discrete Cauchy-Green deformation tensors and discrete Jacobians. The resulting integrator and its conservative properties were illustrated with numerical tests both in 2D and 3D. This paper only presented the very first steps towards the application of variational integrators to fluid-structure interaction problems. The discrete setting developed here opens several directions of research such as the inclusion of dynamic mesh update, inspired by arbitrary Lagrangian-Eulerian methods, and the treatment of h-adaptivity, by appropriate extensions of the variational discretization. One can also explore extensions of our model to take into account of the possible disconnection of fluid cells, which is useful to model the fragmentation of the flow after the impact. Finally, appropriate extensions of the discrete variational formulation need to be developed in order to treat other classes of continuum models, such as viscous-plastic materials.

# A Appendix

## A.1 2D discrete hyperelastic models

To compute the partial derivatives $D_k \mathcal{L}_{a,b}^j$ of the discrete Lagrangian (34) we note

$$D_k \left( \rho_0^e \frac{1}{4} \sum_{\ell=1}^{4} W^e(\mathbf{C}_\ell(\mathcal{D}_{a,b}^j)) \right) \cdot \delta\varphi = \frac{1}{4} \sum_{\ell=1}^{4} \mathrm{Tr} \left( \mathbf{S}_\ell(\mathcal{D}_{a,b}^j) \mathbf{F}_\ell(\mathcal{D}_{a,b}^j)^\mathrm{T} D_k \mathbf{F}_\ell(\mathcal{D}_{a,b}^j) \cdot \delta\varphi \right)$$

for all $\ell = 1, ..., 4$ and all $k = 1, 3, 5, 7$. We observe that the 16 expressions $D_k \mathbf{F}_\ell(\mathcal{D}_{a,b}^j)$ are independent of the value of the field $\varphi_d$ and on the location of the spacetime node (if the spacing $\Delta s_a$ and $\Delta s_b$ are constant), i.e., we can write $D_k \mathbf{F}_\ell(\mathcal{D}_{a,b}^j) = D_k \mathbf{F}_\ell$. For instance, for $\ell = 1$, we have

$$D_1 \mathbf{F}_1 \cdot \delta\varphi = \left[ -\frac{1}{\Delta s_a} \delta\varphi, -\frac{1}{\Delta s_b} \delta\varphi \right]$$

$$D_3 \mathbf{F}_1 \cdot \delta\varphi = \left[ \frac{1}{\Delta s_a} \delta\varphi, 0 \right]$$

$$D_5 \mathbf{F}_1 \cdot \delta\varphi = \left[ 0, \frac{1}{\Delta s_b} \delta\varphi \right]$$

$$D_7 \mathbf{F}_1 \cdot \delta\varphi = [0, 0]$$
and for $\ell = 2$, we have

$$D_1 \mathbf{F}_2 \cdot \delta\varphi = \begin{bmatrix} 0, & \frac{1}{\Delta s_a} \delta\varphi \\ \frac{1}{\Delta s_b} \delta\varphi, & 0 \end{bmatrix}$$ $$D_7 \mathbf{F}_2 \cdot \delta\varphi = \begin{bmatrix} 0, & \frac{1}{\Delta s_b} \delta\varphi \\ \frac{1}{\Delta s_a} \delta\varphi, & 0 \end{bmatrix}$$

$$D_3 \mathbf{F}_2 \cdot \delta\varphi = \begin{bmatrix} -\frac{1}{\Delta s_b} \delta\varphi, & -\frac{1}{\Delta s_a} \delta\varphi \\ \frac{1}{\Delta s_a} \delta\varphi, & 0 \end{bmatrix}$$ $$D_5 \mathbf{F}_2 \cdot \delta\varphi = [0, 0],$$

similarly for $\ell = 3, 4$. We can thus write

$$D_k \left( \rho_0^e \frac{1}{4} \sum_{\ell=1}^{4} W^e (\mathbf{C}_\ell(\mathbb{D}_{a,b}^j)) \right) \cdot \delta\varphi = \frac{1}{4} \sum_{\ell=1}^{4} \left\langle \mathbf{F}_\ell(\mathbb{D}_{a,b}^j) \mathbf{S}_\ell(\mathbb{D}_{a,b}^j), D_k \mathbf{F}_\ell \cdot \delta\varphi \right\rangle \\ = \frac{1}{4} \sum_{\ell=1}^{4} \mathfrak{D}_k^\ell (\mathbf{F}_\ell(\mathbb{D}_{a,b}^j) \mathbf{S}_\ell(\mathbb{D}_{a,b}^j)) \cdot \delta\varphi,$$

where in the last equality we defined the operator $\mathfrak{D}_k^\ell$, $k = 1, 3, 5, 7$, $\ell = 1, 2, 3, 4$, by duality.

From this, the partial derivatives $D_k \mathcal{L}_{a,b}^j$ of the discrete Lagrangian can be computed. For instance for $k = 1$ one has

$$\frac{1}{\mathrm{vol}(\mathbb{D})} D_1 \mathcal{L}_{a,b}^j = -\rho_0^e \frac{1}{4\Delta t} \mathbf{v}_{a,b}^j - \frac{1}{4} \sum_{\ell=1}^4 \mathfrak{D}_1^\ell (\mathbf{F}_\ell(\mathbb{D}_{a,b}^j) \mathbf{S}_\ell(\mathbb{D}_{a,b}^j)) - \frac{1}{4\Delta t} \rho_0^e \nabla \Pi(\varphi_{a,b}^j).$$

Using this result, one then derives (42) from (69).

## A.2 3D discrete Jacobian

$$J_2(\mathbb{D}_{a,b,c}^j) = (\mathbf{F}_{2;a+1,b,c}^j \times \mathbf{F}_{4;a+1,b,c}^j) \cdot \mathbf{F}_{3;a+1,b,c}^j,$$

$$J_3(\mathbb{D}_{a,b,c}^j) = (\mathbf{F}_{5;a+b+1,c}^j \times \mathbf{F}_{1;a+b+1,c}^j) \cdot \mathbf{F}_{3;a+b+1,c}^j,$$

$$J_4(\mathbb{D}_{a,b,c}^j) = (\mathbf{F}_{2;a,b,c+1}^j \times \mathbf{F}_{1;a,b,c+1}^j) \cdot \mathbf{F}_{6;a,b,c+1}^j,$$

$$J_5(\mathbb{D}_{a,b,c}^j) = (\mathbf{F}_{4;a+1,b+1,c}^j \times \mathbf{F}_{5;a+1,b+1,c}^j) \cdot \mathbf{F}_{3;a+1,b+1,c}^j,$$

$$J_6(\mathbb{D}_{a,b,c}^j) = (\mathbf{F}_{1;a+b+1,c+1}^j \times \mathbf{F}_{5;a+b+1,c+1}^j) \cdot \mathbf{F}_{6;a+b+1,c+1}^j,$$

$$J_7(\mathbb{D}_{a,b,c}^j) = (\mathbf{F}_{4;a+1,b,c+1}^j \times \mathbf{F}_{2;a+1,b,c+1}^j) \cdot \mathbf{F}_{6;a+1,b,c+1}^j,$$

$$J_8(\mathbb{D}_{a,b,c}^j) = (\mathbf{F}_{5;a+1,b+1,c+1}^j \times \mathbf{F}_{4;a+1,b+1,c+1}^j) \cdot \mathbf{F}_{6;a+1,b+1,c+1}^j.$$
### A.3 3D discrete Cauchy-Green deformation tensor

$$ C_2(\mathbf{F}_{a,b,c}^j) = \begin{bmatrix} \langle \mathbf{F}_{2;a+1,b,c}^j, \mathbf{F}_{2;a+1,b,c}^j \rangle & \langle \mathbf{F}_{2;a+1,b,c}^j, \mathbf{F}_{4;a+1,b,c}^j \rangle & \langle \mathbf{F}_{2;a+1,b,c}^j, \mathbf{F}_{3;a+1,b,c}^j \rangle \\ \langle \mathbf{F}_{4;a+1,b,c}^j, \mathbf{F}_{2;a+1,b,c}^j \rangle & \langle \mathbf{F}_{4;a+1,b,c}^j, \mathbf{F}_{4;a+1,b,c}^j \rangle & \langle \mathbf{F}_{4;a+1,b,c}^j, \mathbf{F}_{3;a+1,b,c}^j \rangle \\ \langle \mathbf{F}_{3;a+1,b,c}^j, \mathbf{F}_{2;a+1,b,c}^j \rangle & \langle \mathbf{F}_{3;a+1,b,c}^j, \mathbf{F}_{4;a+1,b,c}^j \rangle & \langle \mathbf{F}_{3;a+1,b,c}^j, \mathbf{F}_{3;a+1,b,c}^j \rangle \end{bmatrix}. $$

$$ C_3(\mathbf{F}_{a,b,c}^j) = \begin{bmatrix} \langle \mathbf{F}_{5;a,b+1,c}^j, \mathbf{F}_{5;a,b+1,c}^j \rangle & \langle \mathbf{F}_{5;a,b+1,c}^j, \mathbf{F}_{1;a,b+1,c}^j \rangle & \langle \mathbf{F}_{5;a,b+1,c}^j, \mathbf{F}_{3;a,b+1,c}^j \rangle \\ \langle \mathbf{F}_{1;a,b+1,c}^j, \mathbf{F}_{5;a,b+1,c}^j \rangle & \langle \mathbf{F}_{1;a,b+1,c}^j, \mathbf{F}_{1;a,b+1,c}^j \rangle & \langle \mathbf{F}_{1;a,b+1,c}^j, \mathbf{F}_{3;a,b+1,c}^j \rangle \\ \langle \mathbf{F}_{3;a,b+1,c}^j, \mathbf{F}_{5;a,b+1,c}^j \rangle & \langle \mathbf{F}_{3;a,b+1,c}^j, \mathbf{F}_{1;a,b+1,c}^j \rangle & \langle \mathbf{F}_{3;a,b+1,c}^j, \mathbf{F}_{3;a,b+1,c}^j \rangle \end{bmatrix}. $$

$$ C_4(\mathbf{F}_{a,b,c}^j) = \begin{bmatrix} \langle \mathbf{F}_{2;a,b,c+1}^j, \mathbf{F}_{2;a,b,c+1}^j \rangle & \langle \mathbf{F}_{2;a,b,c+1}^j, \mathbf{F}_{1;a,b,c+1}^j \rangle & \langle \mathbf{F}_{2;a,b,c+1}^j, \mathbf{F}_{6;a,b,c+1}^j \rangle \\ \langle \mathbf{F}_{1;a,b,c+1}^j, \mathbf{F}_{2;a,b,c+1}^j \rangle & \langle \mathbf{F}_{1;a,b,c+1}^j, \mathbf{F}_{1;a,b,c+1}^j \rangle & \langle \mathbf{F}_{1;a,b,c+1}^j, \mathbf{F}_{6;a,b,c+1}^j \rangle \\ \langle \mathbf{F}_{6;a,b,c+1}^j, \mathbf{F}_{2;a,b,c+1}^j \rangle & \langle \mathbf{F}_{6;a,b,c+1}^j, \mathbf{F}_{1;a,b,c+1}^j \rangle & \langle \mathbf{F}_{6;a,b,c+1}^j, \mathbf{F}_{6;a,b,c+1}^j \rangle \end{bmatrix}. $$

$$ C_5(\mathbf{F}_{a,b,c}^j) = \begin{bmatrix} \langle \mathbf{F}_{4;a+1,b+1,c}^j, \mathbf{F}_{4;a+1,b+1,c}^j \rangle & \langle \mathbf{F}_{4;a+1,b+1,c}^j, \mathbf{F}_{5;a+1,b+1,c}^j \rangle & \langle \mathbf{F}_{4;a+1,b+1,c}^j, \mathbf{F}_{3;a+1,b+1,c}^j \rangle \\ \langle \mathbf{F}_{5;a+1,b+1,c}^j, \mathbf{F}_{4;a+1,b+1,c}^j \rangle & \langle \mathbf{F}_{5;a+1,b+1,c}^j, \mathbf{F}_{5;a+1,b+1,c}^j \rangle & \langle \mathbf{F}_{5;a+1,b+1,c}^j, \mathbf{F}_{3;a+1,b+1,c}^j \rangle \\ \langle \mathbf{F}_{3;a+1,b+1,c}^j, \mathbf{F}_{4;a+1,b+1,c}^j \rangle & \langle \mathbf{F}_{3;a+1,b+1,c}^j, \mathbf{F}_{5;a+1,b+1,c}^j \rangle & \langle \mathbf{F}_{3;a+1,b+1,c}^j, \mathbf{F}_{3;a+1,b+1,c}^j \rangle \end{bmatrix}. $$

$$ C_6(\mathbf{F}_{a,b,c}^j) = \begin{bmatrix} \langle \mathbf{F}_{1;a,b+1,c+1}^j, \mathbf{F}_{1;a,b+1,c+1}^j \rangle & \langle \mathbf{F}_{1;a,b+1,c+1}^j, \mathbf{F}_{5;a,b+1,c+1}^j \rangle & \langle \mathbf{F}_{1;a,b+1,c+1}^j, \mathbf{F}_{6;a,b+1,c+1}^j \rangle \\ \langle \mathbf{F}_{5;a,b+1,c+1}^j, \mathbf{F}_{1;a,b+1,c+1}^j \rangle & \langle \mathbf{F}_{5;a,b+1,c+1}^j, \mathbf{F}_{5;a,b+1,c+1}^j \rangle & \langle \mathbf{F}_{5;a,b+1,c+1}^j, \mathbf{F}_{6;a,b+1,c+1}^j \rangle \\ \langle \mathbf{F}_{6;a,b+1,c+1}^j, \mathbf{F}_{1;a,b+1,c+1}^j \rangle & \langle \mathbf{F}_{6;a,b+1,c+1}^j, \mathbf{F}_{5;a,b+1,c+1}^j \rangle & \langle \mathbf{F}_{6;a,b+1,c+1}^j, \mathbf{F}_{6;a,b+1,c+1}^j \rangle \end{bmatrix}. $$

$$ C_7(\mathbf{F}_{a,b,c}^j) = \begin{bmatrix} \langle \mathbf{F}_{4;a+1,b,c+1}^j, \mathbf{F}_{4;a+1,b,c+1}^j \rangle & \langle \mathbf{F}_{4;a+1,b,c+1}^j, \mathbf{F}_{2;a+1,b,c+1}^j \rangle & \langle \mathbf{F}_{4;a+1,b,c+1}^j, \mathbf{F}_{6;a+1,b,c+1}^j \rangle \\ \langle \mathbf{F}_{2;a+1,b,c+1}^j, \mathbf{F}_{4;a+1,b,c+1}^j \rangle & \langle \mathbf{F}_{2;a+1,b,c+1}^j, \mathbf{F}_{2;a+1,b,c+1}^j \rangle & \langle \mathbf{F}_{2;a+1,b,c+1}^j, \mathbf{F}_{6;a+1,b,c+1}^j \rangle \\ \langle \mathbf{F}_{6;a+1,b,c+1}^j, \mathbf{F}_{4;a+1,b,c+1}^j \rangle & \langle \mathbf{F}_{6;a+1,b,c+1}^j, \mathbf{F}_{2;a+1,b,c+1}^j \rangle & \langle \mathbf{F}_{6;a+1,b,c+1}^j, \mathbf{F}_{6;a+1,b,c+1}^j \rangle \end{bmatrix}. $$

$$ C_8(\mathbf{F}_{a,b,c}^j) = \begin{bmatrix} \langle \mathbf{F}_{5;a+1,b+1,c+1}^j, \mathbf{F}_{5;a+1,b+1,c+1}^j \rangle & \langle \mathbf{F}_{5;a+1,b+1,c+1}^j, \mathbf{F}_{4;a+1,b+1,c+1}^j \rangle & \langle \mathbf{F}_{5;a+1,b+1,c+1}^j, \mathbf{F}_{6;a+1,b+1,c+1}^j \rangle \\ \langle \mathbf{F}_{4;a+1,b+1,c+1}^j, \mathbf{F}_{5;a+1,b+1,c+1}^j \rangle & \langle \mathbf{F}_{4;a+1,b+1,c+1}^j, \mathbf{F}_{4;a+1,b+1,c+1}^j \rangle & \langle \mathbf{F}_{4;a+1,b+1,c+1}^j, \mathbf{F}_{6;a+1,b+1,c+1}^j \rangle \\ \langle \mathbf{F}_{6;a+1,b+1,c+1}^j, \mathbf{F}_{5;a+1,b+1,c+1}^j \rangle & \langle \mathbf{F}_{6;a+1,b+1,c+1}^j, \mathbf{F}_{4;a+1,b+1,c+1}^j \rangle & \langle \mathbf{F}_{6;a+1,b+1,c+1}^j, \mathbf{F}_{6;a+1,b+1,c+1}^j \rangle \end{bmatrix}. $$
A.4 3D derivatives of impenetrability constraints (71)

$$
D_1 \Psi_{\text{im1}} = -(\varphi_{a-1,b,c}^j - \varphi_{a,b,c}^j) \times (\varphi_{a,b-1,c}^j - \varphi_{a,b,c}^j) - (\varphi_{a,b-1,c}^j - \varphi_{a,b,c}^j) \times (\varphi_{d,e,f}^j - \varphi_{a,b,c}^j) \\ - (\varphi_{d,e,f}^j - \varphi_{a,b,c}^j) \times (\varphi_{a-1,b,c}^j - \varphi_{a,b,c}^j),
$$

$$
D_2 \Psi_{\text{im1}} = (\varphi_{a,b-1,c}^j - \varphi_{a,b,c}^j) \times (\varphi_{d,e,f}^j - \varphi_{a,b,c}^j),
$$

$$
D_3 \Psi_{\text{im1}} = (\varphi_{d,e,f}^j - \varphi_{a,b,c}^j) \times (\varphi_{a-1,b,c}^j - \varphi_{a,b,c}^j),
$$

$$
D_4 \Psi_{\text{im1}} = (\varphi_{a-1,b,c}^j - \varphi_{a,b,c}^j) \times (\varphi_{a,b-1,c}^j - \varphi_{a,b,c}^j),
$$

$$
D_1 \Psi_{\text{im2}} = -(\varphi_{a,b-1,c}^j - \varphi_{a,b,c}^j) \times (\varphi_{a+1,b,c}^j - \varphi_{a,b,c}^j) - (\varphi_{d+1,e,f}^j - \varphi_{a,b,c}^j) \times (\varphi_{a,b-1,c}^j - \varphi_{a,b,c}^j) \\
- (\varphi_{a+1,b,c}^j - \varphi_{a,b,c}^j) \times (\varphi_{d+1,e,f}^j - \varphi_{a,b,c}^j),
$$

$$
D_2 \Psi_{\text{im2}} = (\varphi_{d+1,e,f}^j - \varphi_{a,b,c}^j) \times (\varphi_{a,b-1,c}^j - \varphi_{a,b,c}^j),
$$

$$
D_3 \Psi_{\text{im2}} = (\varphi_{a+1,b,c}^j - \varphi_{a,b,c}^j) \times (\varphi_{d+1,e,f}^j - \varphi_{a,b,c}^j),
$$

$$
D_4 \Psi_{\text{im2}} = (\varphi_{a,b-1,c}^j - \varphi_{a,b,c}^j) \times (\varphi_{a+1,b,c}^j - \varphi_{a,b,c}^j),
$$

$$
D_1 \Psi_{\text{im3}} = -(\varphi_{a,b+1,c}^j - \varphi_{a,b,c}^j) \times (\varphi_{a-1,b,c}^j - \varphi_{a,b,c}^j) - (\varphi_{d,e+1,f}^j - \varphi_{a,b,c}^j) \times (\varphi_{a,b+1,c}^j - \varphi_{a,b,c}^j), \\
- (\varphi_{a-1,b,c}^j - \varphi_{a,b,c}^j) \times (\varphi_{d,e+1,f}^j - \varphi_{a,b,c}^j),
$$

$$
D_2 \Psi_{\text{im3}} = (\varphi_{d,e+1,f}^j - \varphi_{a,b,c}^j) \times (\varphi_{a,b+1,c}^j - \varphi_{a,b,c}^j),
$$

$$
D_3 \Psi_{\text{im3}} = (\varphi_{a-1,b,c}^j - \varphi_{a,b,c}^j) \times (\varphi_{d,e+1,f}^j - \varphi_{a,b,c}^j),
$$

$$
D_4 \Psi_{\text{im3}} = (\varphi_{a,b+1,c}^j - \varphi_{a,b,c}^j) \times (\varphi_{a-1,b,c}^j - \varphi_{a,b,c}^j),
$$

$$
D_1 \Psi_{\text{im4}} = -(\varphi_{a+1,b,c}^j - \varphi_{a,b,c}^j) \times (\varphi_{a,b+1,c}^j - \varphi_{a,b,c}^j) - (\varphi_{d+1,e+1,f}^j - \varphi_{a,b,c}^j) \times (\varphi_{a+1,b,c}^j - \varphi_{a,b,c}^j) \\
- (\varphi_{a,b+1,c}^j - \varphi_{a,b,c}^j) \times (\varphi_{d+1,e+1,f}^j - \varphi_{a,b,c}^j),
$$

$$
D_2 \Psi_{\text{im4}} = (\varphi_{a,b+1,c}^j - \varphi_{a,b,c}^j) \times (\varphi_{d+1,e+1,f}^j - \varphi_{a,b,c}^j),
$$

$$
D_3 \Psi_{\text{im4}} = (\varphi_{d+1,e+1,f}^j - \varphi_{a,b,c}^j) \times (\varphi_{a+1,b,c}^j - \varphi_{a,b,c}^j),
$$

$$
D_4 \Psi_{\text{im4}} = (\varphi_{a+1,b,c}^j - \varphi_{a,b,c}^j) \times (\varphi_{a,b+1,c}^j - \varphi_{a,b,c}^j).
$$

References

[1] Angoshtari, A. and Yavari, A. [2013], A geometric structure-preserving discretization scheme for incompressible linearized elasticity. *Comput. Methods in Appl. Mech. Eng.*, **259**, 130–153.

[2] Ariza, M.P. and Ortiz, M. [2005], Discrete crystal elasticity and discrete dislocations in crystals, *Arch. Ration. Mech. Anal.*, **178**(2), 149–226.

[3] Cirak, F. and M. West [2005], Decomposition contact response (DCR) for explicit finite element dynamics, *Int. J. Numer. Meth. Engng*, **64**, 1078–1110.

[4] Cottet, G.-H., Maitre, E., and Milcent, T. [2008], Eulerian formulation and level set models for incompressible fluid-structure interaction, *ESAIM-Math. Model. Numer. Anal.* **42**(3),471–492.

[5] Demoures, F., Gay-Balmaz, F., and Ratiu, T. S. [2014], Multisymplectic variational integrator and space/time symplecticity, *Anal. Appl.*, **14**(3), 341–391.

[6] Demoures, F., Gay-Balmaz, F., and Ratiu, T. S. [2016], Multisymplectic variational integrators for nonsmooth Lagrangian continuum mechanics, *Forum Math. Sigma*, **4**, e19, 54p.

[7] Demoures, F., Gay-Balmaz, F., Desbrun, M., Ratiu, T. S., and Alejandro, A. [2017], A multisymplectic integrator for elastodynamic frictionless impact problems, *Comput. Methods in Appl. Mech. Eng.*, **315**, 1025–1052.
[8] Demoures, F., Gay-Balmaz, F., Kobilarov, M. and Ratiu, T. S. [2014], Multisymplectic Lie group variational integrators for a geometrically exact beam in $\mathbb{R}^3$, *Commun. Nonlinear Sci. Numer. Simulat.*, **19**(10), 3492–3512.

[9] Demoures, F. and Gay-Balmaz F. [2021], Multisymplectic variational integrators for barotropic and incompressible fluid models with constraints, submitted arXiv:2102.10907

[10] Gay-Balmaz, F., Marsden, J. E., and Ratiu, T. S. [2012], Reduced variational formulations in free boundary continuum mechanics, *Journal of nonlinear science* **22**(4), 463–497.

[11] Gotay, M. J., Isenberg, J., Marsden, J. E., Montgomery, R., Sniatycki, J., and Yasskin, P. B. [1997], Momentum maps and classical fields. Part I: Covariant field theory. arXiv:physics/389801019v2.

[12] Glowinski, R., Pan, T-W., and Péraux, J. A. [1997], Lagrange multiplier/fictitious domain method for the numerical simulation of incompressible viscous flow around moving rigid bodies: (i) case where the rigid body motions are known a priori. *CR Acad Sci Paris*, **25**(5), 361–9.

[13] Hirt, C. W., Amsden, A. A., and Cook, J. L. [1997], An arbitrary Lagrangian-Eulerian computing method for all flow speeds, *Journal of Computational Physics*, **135**(2), 203–216.

[14] Hou, G., Wang, J., and Layton, A. [2012], Numerical Methods for Fluid-Structure Interaction - A Review, *Commun. Comput. Phys.*, **12**(2), 337–377.

[15] Kamrin, K., C. H. Rycroft, and J. C. Nave [2012], Reference map technique for finite-strain elasticity and fluid-solid interaction, *Journal of the Mechanics and Physics of Solids*, **60**(11), 1952–1969.

[16] Kanso, E., Marsden, J. E., Rowley, C., and Melli-Huber, J.-H. [2005], Locomotion of Articulated Bodies in a Perfect Fluid. *J. Nonlinear Sci.*, **15**, 255–289.

[17] Khayyer, A., Gotoh, H., Falahaty, H., and Shimizu, Y. [2018], An enhanced ISPH-SPH coupled method for simulation of incompressible fluid-elastic structure interactions, *Computer Physics Communications*, **232**, 139–164.

[18] Lefrançois, E. and Boufflet, J.-P. [2010], An introduction to fluid-structure interaction: Application to the piston problem, *SIAM Review*, **52**, 747–767.

[19] Lew A., Marsden, J. E., Ortiz, M., and West, M. [2003], Asynchronous variational integrators, *Arch. Rational Mech. Anal.*, **167**(2), 85–146.

[20] van Loon, R., Anderson, P.D., van de Vosse, F.N., and Sherwin, S.J. [2007], Comparison of various fluid-structure interaction methods for deformable bodies, *Comput. Struct.*, **85**(11-14), 833–843.
MaVe2006

[21] Marsden, J. E. and Hughes, T. J. R. *Mathematical foundations of elasticity*. Prentice-Hall, 1983.

[22] Marsden, J. E., Patrick, G. W., and Shkoller, S. [1998], Multisymplectic geometry, variational integrators and nonlinear PDEs, *Comm. Math. Phys.*, **199**, 351–395.

[23] Marsden, J. E., Pekarsky, S. Shkoller, S., and West, M. [2001], Variational methods, multisymplectic geometry and continuum mechanics, *J. Geom. Phys.*, **38**, 253-284.

[24] Marsden, J. E. and West M. [2001], Discrete mechanics and variational integrators, *Acta Numer.* **10**, 357–514.

[25] Mooney, M. [1940], A theory of large elastic deformation, *J. Appl. Phys.*, **11**(6), 582–592.

[26] Nocedal, J. and Wright, S. J. [2006], *Numerical Optimization*, Springer Series in Operations Research and Financial Engineering, Springer.

[27] Peskin, C. S. and McQueen, D. M. [1989], A three-dimensional computational method for blood flow in the heart i. immersed elastic fibers in a viscous incompressible fluid. *J. Comp. Phys.*, **81**, 372–405.

[28] Peskin, C. S. [2002], The immersed boundary method, *Acta Numerica*, **11**, 479–517.

[29] Rivlin, R.S. [1948], Large elastic deformations of isotropic materials. IV. Further developments of the general theory, *Philos. Trans. Royal Soc.*, **A24I**(6), 379–397.
[30] Rivlin, R.S. [1949], Large elastic deformations of isotropic materials. V. The problem of flexure, *Proc. Royal Soc. Lond.*, **A195**(6), 463–473.

[31] Rivlin, R.S. [1949b], Large clastic deformations of isotropic materials. VI. Further results in the theory of torsion, shear and flexure, *Philos. Trans. Royal Soc.*, **A242**(6), 173–195.

[32] Rockafellar, R. T. [1993], *Lagrange multipliers and optimality*, *SIAM Review* **35**(2), 183–238.

[33] Rockafellar, R. T. and Wets, R. J.-B. [1998], *Variational Analysis*, Grundlehren der Mathematischen Wissenschaften, **317**, Springer-Verlag, Berlin, 1998.

[34] Simo, J. C. and Laursen, T. A. [1992], An augmented Lagrangian treatment of contact problems involving friction, *Computers and Structures*, **42**(1), 97–116.

[35] Yavari, A. [2008], On geometric discretization of elasticity, *J. Math. Phys.*, **49**, 022901.