{'fraction_non_alphanumeric': 0.04521306915342347, 'fraction_numerical': 0.016165439248500876, 'mean_word_length': 4.199146514935989, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 28, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Unified discrete multisymplectic Lagrangian formulation for hyperelastic solids and barotropic fluids

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.


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.
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 rubberlike 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.

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.

Description of hyperelastic bodies
We assume that the elastic body is represented by a compact manifold B ⊂ R 3 with piecewise smooth boundary, and moves in the ambient space M = R 3 . The motion of the body is described by a time dependent configuration map ϕ : R × B → M which indicates the location m = ϕ(t, X) ∈ M at time t of the material point X ∈ B. The deformation gradient is denoted F(t, X) = ∇ϕ(t, X), given in coordinates by F a i = ϕ a ,i , with X i , i = 1, 2, 3 the Cartesian coordinates on B and m a , a = 1, 2, 3 the Cartesian coordinates on 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 W depending on the points X ∈ B and the gradient deformation F(t, X) such that the first Piola-Kirchhoff stress tensor of the body in the configuration ϕ(t, X) is P = ρ 0 ∂W ∂F , where ρ 0 is the mass density of the body in the Lagrangian description. Given Riemannian metrics G and g on B and M, the stored energy function of an hyperelastic material with configuration map ϕ : B → M (at some fixed time) takes the general form W(X, F(X), G(X), g(ϕ(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 i.e., we have W(X, F(X), G(X), g(ϕ(X))) = W (X, G(X), C(X)), for some function W , where ϕ * denotes the pull-back of a tensor field by ϕ.
The material (or Lagrangian) strain tensor E and the second Piola-Kirchhoff stress tensor S are defined respectively by The elasticity tensor on the reference configuration is defined by For simplicity, we will consider the Euclidean case, where G ij = δ ij and g ab = δ ab in Cartesian coordinates, in which case we can write where ·, · denotes the Euclidean inner product on R n .

Euler-Lagrange equations for an hyperelastic body
The Lagrangian density of an hyperelastic body evaluated on a configuration map ϕ(t, X) is with Π 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 for variations δϕ of the body configuration map with δϕ = 0 at t = 0, T . It yields the Euler-Lagrange equations d dt together with the natural boundary conditions ∂L ∂F a i n i δϕ a = 0 on ∂B for the allowed variations δϕ at the boundary, with F = ∇ϕ. For the Lagrangian (5) we compute hence the Euler-Lagrange equations (6) for the Lagrangian (5) gives the equations of motion for an hyperelastic body as 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.

Isotropic hyperelastic materials
A frame indifferent hyperelastic material is isotropic if and only if the stored energy function can be written as W (X, G(X), C(X)) = Φ(X, λ 1 (X), λ 2 (X), λ 3 (X)), where Φ is a symmetric function of the principal stretches λ 1 , λ 2 , λ 3 1 . This is equivalent to state that W is a function of the invariants I 1 , I 2 , I 3 of C defined by 2 I 1 (C) = Tr(C), The derivative of the invariants with respect to C are given by from which the second Piola-Kirchhoff stress tensor (3) reads We refer to [21, §3.5] for more details.
1 Given the polar decomposition F = RU of the gradient deformation, the principal stretches λ1, λ2, λ3 are the eigenvalues of U. In particular λ 2 1 , λ 2 2 , λ 2 3 are the eigenvalues of C = U 2 . 2 For simplicity we didn't write explicitly the dependence of the invariants on G and assumed Gij = δij.
The Lagrangian density of a frame indifferent isotropic hyperelastic material is therefore while the equations of motion take the form Incompressible materials can be treated by inserting the constraint I 3 = 1 in the Hamilton principle as where µ 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 with positive material constants C 1 , C 2 .

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 π Y,X : Y → X . The configuration fields ϕ of the theory are sections of this fiber bundle, i.e., they are smooth maps ϕ : X → Y such that π Y,X • ϕ = id X , where id X denotes the identity map on X .
The first jet bundle of the configuration bundle π Y,X : Y → X is the field theoretic analogue of the tangent bundle of classical mechanics. It is defined as the fiber bundle π J 1 Y,Y : J 1 Y → Y whose fiber at y ∈ Y consists of linear maps γ : T x X → T y Y satisfying T π Y,X • γ = id TxX , where x = π Y,X (y). The derivative of a field ϕ can be regarded as a section of the fiber bundle π J 1 Y,X = π Y,X • π J 1 Y,Y : J 1 Y → X , by considering its first jet extension x ∈ X → j 1 ϕ( A Lagrangian density is a smooth bundle map L : J 1 Y → Λ n+1 X over X , where Λ n+1 X → X is the vector bundle of (n + 1)-forms on X , with dim X = n + 1. The associated action functional is given by and the Euler-Lagrange equations follow from the stationary condition for δS(ϕ) = 0 for appropriate boundary conditions. For continuum mechanics, the configuration bundle is the trivial fiber bundle where B is the reference configuration of the elastic body and M is the ambient space.
We have the variables x = (t, X) ∈ X and y = (x, m) = (t, X, m) ∈ Y. A section of this bundle is a map ϕ : X → X × M, whose first component is id X . It is canonically identified with a map ϕ : X = R × B → M referred to as the elastic body configuration in §2.1. The first jet extension of ϕ is j 1 ϕ(t, X) = (ϕ(t, X),φ(t, X), ∇ϕ(t, X)) and the Lagrangian density reads Multisymplectic form formula and Noether theorem. A solution ϕ of the Euler-Lagrange equation associated to L : J 1 Y → Λ n+1 X satisfies the multisymplectic form formula given by for all open subset U ⊂ X with piecewise smooth boundary, and for all vector fields V, W solutions of the first variation of the Euler-Lagrange equations at ϕ, see [22]. In this formula, Ω L denotes the Cartan (n + 2)-form on J 1 Y associated to L and j 1 V, j 1 W are the first jet extension of V, W to J 1 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 Y, assume that this action covers a diffeomorphism of X , and assume that the Lagrangian density L is G-equivariant. Then a solution ϕ of the Euler-Lagrange equation satisfies the covariant Noether theorem for all open subset U ⊂ X with piecewise smooth boundary and for all ξ ∈ g, the Lie algebra of G, with J L : J 1 Y → g * ⊗ Λ n J 1 Y the covariant momentum map associated to L and to the Lie group action, see [11] for details. The properties (16) and (17) admit discrete versions that will be considered later.

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.

Multisymplectic discretization and the discrete Cauchy-Green tensor
Discrete Lagrangian setting. Consider the configuration bundle of continuum mechanics Y = X × M → X = R × B. In this section we assume that B is a domain in R 2 and take M = R 2 . To discretize the geometric setting, one first considers a discrete parameter space U d and a discrete base-space configuration, which is a one-to-one map whose image X d is the discrete spacetime. The discrete fields are the sections of the discrete configuration bundle π d : Y d = X d × M → X d , which can be identified with maps ϕ d : X d → M. The discrete field and the discrete base-space configuration can be simultaneously described by introducing the discrete configuration, which is a map 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 and denote by ϕ j a,b := ϕ d (s j a,b ) the value of the discrete field at s j a,b .
Denoting by U d the set of all parallelepipeds in the discrete parameter space U d , we write X d := φ X d U d and recall that the discrete version of the first jet bundle is given by Given a discrete field ϕ d , its first jet extension is the section of (20) defined by which associates to each parallelepiped, the values of the field at its nodes. Finally, we recall that a discrete Lagrangian is a map from which the discrete action functional is constructed as Discrete Cauchy-Green deformation tensor. Given a discrete base-space configuration φ X d and a discrete field ϕ d , the following four vectors F j ;a,b ∈ R 2 , = 1, 2, 3, 4 are defined at each node (j, a, b) ∈ U d , see Fig. 3 on the right: at time t j and the associated vectors F ;a,b , = 1, 2, 3, 4.
Definition 3.1 The discrete deformation gradients of a discrete field ϕ d at the parallelepiped j a,b are the four 2 × 2 matrices F ( j a,b ), = 1, 2, 3, 4, defined by Note that the discrete deformation gradients are defined at each of the four nodes of j a,b that are associated to time t j , see Fig. 3 on the left for the ordering = 1, 2, 3, 4 with respect to the nodes. Also, for each , the discrete deformation gradients F ( ), = 1, 2, 3, 4, depend on the discrete field ϕ d only through its first jet extension j 1 ϕ d ( ) at the given , see (21).
It is assumed that the discrete field ϕ 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 ϕ d at the parallelepiped j a,b are the four 2 × 2 symmetric matrices With the assumption on the discrete field mentioned above, the discrete deformation gradients and the discrete Cauchy-Green deformation tensors are maps defined on X d , with Sym + (2) the set of symmetric positive definite 2 × 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 C ij = ϕ ,i , ϕ ,j . Proposition 3.3 Given a discrete field ϕ d : X d → M, the discrete Cauchy-Green deformation tensors are explicitly given as , in terms of the vector F j ;a,b defined in (22).
Discrete Jacobian. We recall from [9] that the discrete Jacobian of a discrete field ϕ d at the parallelepiped j a,b are the four numbers J ( j a,b ), = 1, 2, 3, 4, defined by 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 for each j a,b ∈ X d and all = 1, 2, 3, 4.

Discrete frame indifferent linear hyperelastic model
In the 2D case, we shall focus on the St. Venant-Kirchhoff model where λ and µ are the first and second Lamé coefficients and E is the material strain tensor defined in (2) with G ij = δ 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 the Poisson ratio and E = 4µ(λ+µ) λ+2µ the Young modulus in 2D.
Discrete stored energy function. The discrete stored energy function W d : is defined in terms of the continuous one as follows Discrete second Piola-Kirchhoff stress. The discrete second Piola-Kirchhoff stress at time t j and spatial position , with the ordering = 1 to = 4, respectively associated to the nodes (j, a, b), for all j a,b ∈ X d .

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 with ρ f 0 the fluid mass density and the specific internal energy barotropic fluid. We recall that Hamilton's principle yields the barotropic fluid equations in the Lagrangian (or material) description as Discrete Lagrangians for the fluid and the elastic body. For the elastic body, we consider discrete Lagrangians of the form where vol( ) is the volume of the parallelepiped ∈ X d . For simplicity, we have assume that ρ e 0 is constant in space. The discrete kinetic energy K d : with with W e the stored energy function of the hyperelastic body. Finally, the discrete potential energy Π d : where Π is the potential energy of the continuous model. We shall focus on the gravitation potential Π(ϕ) = g, ϕ , with 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 where now W f d : J 1 Y d → R is the discrete internal energy of the fluid given by with W f ( f 0 , J) the internal energy of the barotropic fluid and J ( j c,d ) the discrete Jacobians associated to j c,d . We assumed that ρ f 0 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 φ X d (j, a, b) = (j∆t, a∆s 1 , b∆s 2 ), for given ∆t, ∆s 1 , ∆s 2 , similarly for the fluid. In this case, we have vol( ) = ∆t∆s 1 ∆s 2 in the discrete Lagrangian and the mass of each 2D cell in φ We consider the discrete action functional where L e d and L f d are the discrete Lagrangians associated to the elastic body and to the fluid. For notational simplicity we have used the same notation ϕ 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 ϕ d ( ), in the order given in (21). From the discrete Hamilton principle δS d (ϕ d ) = 0, we get the discrete Euler-Lagrange equations in the general form for both L e d and L f d , where we used the abbreviate notation L j a,b := L(j 1 ϕ d ( j a,b )). By using the expressions of the discrete Lagrangians given in (34) and (38), the general form (69) yields n = e, f, see Appendix A.1 for the derivation, as well as for the definition of the operators D k , k = 1, 3, 5, 7. Note that these are the discrete equations of motion corresponding to (8). In (42), S is the discrete second Piola-Kirchhoff tensor of the model considered, see (31). For the St. Venant-Kirchhoff model one has S = λTr(E)l + 2µE while for the barotropic fluid S = −P W JC −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 ϕ 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) (R θ , u) on a discrete field given by ϕ d → R θ ϕ d + u with R θ the rotation of angle θ. The discrete covariant momentum maps J p L d : J 1 Y d → se(2) * , p = 1, ..., 8, associated to the SE(2) action and to a discrete Lagrangian L d : In (43), ϕ (p) is the p-th component of j 1 ϕ d , see (21), D p denotes the partial derivative with respect to the p-th argument of L d , and L j a,b := L d j 1 ϕ d ( j a,b ) . From the discrete covariant momentum maps (43), the discrete classical momentum map is given by with ϕ j the collection of all positions at time t j , see [5]. For the discrete Lagrangians (34) and (38), we get the expressions with similarly for the fluid.
In absence of gravitational potential, the discrete Lagrangians (34) and (38) are SE(2) invariant. Indeed, when ϕ d → R θ ϕ d + u the associated discrete deformation gradients (23) transform as F → R θ F , 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 p L d , given as 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 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 ( ) − 1 = 0, for all ∈ X d and all = 1, 2, 3, 4, which mimics the properties of rubber. The associated quadratic penalty function is given by where r is the penalty parameter. The second constraint is the impenetrability condition between the fluid node ϕ j c,d flowing above the elastic body nodes ϕ j a,b , ϕ j a+1,b , which is defined as follows where R π/2 is a rotation of +π/2. In this expression R π/2 (ϕ j a+1,b − ϕ j a,b ) 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 Ψ im in (47) is evaluated on all the boundary nodes ϕ j c,d of the fluid and on the two boundary nodes of the elastic body that are the closest to ϕ j c,d . The associated quadratic penalty term is given by The derivatives of Φ im with respect to the nodes ϕ j a,b , ϕ j a+1,b , ϕ j c,d are given by which express the directions of the reactions forces.
The resulting discrete action functional is found as 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 ϕ j c,d , the nodes ϕ j a,b , ϕ j a+1,b are chosen as explained above. For the boundary nodes ϕ j c,d that are not in contact, the corresponding term vanishes by definition of Φ 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.

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 We shall use the expression (51) for the treatment of an isentropic perfect fluid, where the value of the constant B = 0 does not affect the dynamics, while it allows to naturally impose the boundary condition P | ∂B = B, with P = A ρ 0 J γ 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 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: 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 J d = (J 1 , J 2 , J 3 ) ∈ 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, ..., D − 1, of the fluid associated to pressures see Appendix A.2 in [9] for details. Similarly, the forces acting on the left at the spacetime nodes (j, 0, d), d = 1, ..., D −1, of the fluid associated to pressures
We consider the explicit integrator described in (42), for the elastic body only, and study the convergence with respect to ∆t and ∆s i , i = 1, 2.
1. Given a fixed mesh, with values ∆s 1 = ∆s 2 = 0.025m, we vary the time-steps as ∆t ∈ {2 × 10 −4 , 10 −4 , 5 × 10 −5 , 2.5 × 10 −5 }. We compute the L 2 -errors in the position ϕ d at time t N = 0.5s, by comparing ϕ d with an "exact solution" obtained with the time-step ∆t ref = 6.25 × 10 −6 s. That is, for each value of ∆t we calculate This yields the following convergence with respect to ∆t The results indicate that the method is second-order accurate with respect to ∆s 1 = ∆s 2 .

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.

3D discrete setting and the discrete Cauchy-Green tensor
Discrete Lagrangian setting. With B a domain in R 3 with piecewise smooth boundary and M = R 3 , the discrete configuration bundle is defined as in the 2D case, see the beginning of §3.1 and Fig. 1 Fig. 8 for the eight pairs of indices in j a,b,c at time t j . We assume that the discrete base-space configuration is of the form and denote by ϕ j a,b,c := ϕ d (s j a,b,c ) the value of the discrete field at s j a,b,c . 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 Discrete Cauchy-Green deformation tensor in 3D. Given a discrete base-space configuration φ X d and a discrete field ϕ d , the following six vectors F j ;a,b,c ∈ R 3 , = 1, ..., 6 are defined at each node (j, a, b, c) ∈ U d , see Fig. 8 on the right:

Definition 4.1
The discrete deformation gradients of a discrete field ϕ d at the parallelepiped j a,b are the eight 3 × 3 matrices F ( j a,b,c ), = 1, ..., 8, defined by F 6 ( j a,b,c ) = F j 1;a,b+1,c+1 F j 5;a,b+1,c+1 F j 6;a,b+1,c+1 , We note that the discrete deformation gradients are defined at each of the eight nodes of j a,b,c that are associated to time t j . We assume that the determinant of these matrices is always positive. The ordering = 1 to = 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 F = RU = VR exactly as in the continuous case.
Discrete Jacobian. The 3D discrete Jacobian of ϕ d associated to j a,b,c at the node (j, a, b, c) is defined by 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.

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.
For each fixed j a,b,c , these invariants depend only of the first jet extension j 1 ϕ d ( j a,b,c ) of the discrete field at the given j a,b,c . Recall that these invariants are related to the coefficients in the characteristic polynomial P (ν) of the 3 × 3 symmetric matrix C j a,b,c as follows . In particular, the eigenvalues ν 1 , ν 2 , ν 3 are defined at each node with respect to each of the eight neighbouring parallelepipeds. For each j a,b,c ∈ X d and = 1, ..., 8, we have and ν 1 = λ 2 1 , ν 2 = λ 2 2 , ν 3 = λ 2 3 , with λ 1 , λ 2 , λ 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. 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 : For each ∈ X 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 at time t j .
3D Discrete second Piola-Kirchhoff stress. It is defined at each node of each ∈ X d as = 1, ..., 8. In particular, for an isotropic material, it takes the form Example 4.6 (Nonlinear Mooney-Rivlin incompressible model) For this model, see (14), the discrete stored energy becomes and the 3D discrete second Piola-Kirchhoff stress is Incompressibility is imposed via a Lagrange multiplier exactly as in the 2D case earlier.

3D Barotropic fluid -elastic body interactions
Discrete Lagrangians for the fluid and the elastic body. The discrete Lagrangian L e d 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 where v j α,β,γ = (ϕ j+1 a,b,c − ϕ j a,b,c )/∆t, we take the discrete stored energy function and the gravitational potential energy density We proceed similarly for the discrete Lagrangian L f d of the fluid, see (38).
Discrete action and discrete Euler-Lagrange equations. As earlier, we assume φ X d (j, a, b, c) = (j∆t, a∆s 1 , b∆s 2 , c∆s 3 ) and the mass of each 3D cell in Φ X d (B d ) is M e = ρ e 0 ∆s 1 ∆s 2 ∆s 3 . The discrete action functional reads Let us denote by D k , k = 1, ..., 16, the partial derivative of the discrete Lagrangians with respect to the k-th component of j 1 ϕ d ( ), in the order given in (55). From the discrete Hamilton principle δS d (ϕ d ) = 0, we get the discrete Euler-Lagrange equations in the general form for both L e d and L f d , where we used the abbreviate notation L j a,b,c := L(j 1 ϕ d ( j a,b,c )). For the discrete Lagrangians based on the discretizations given in (65)-(67), the discrete Euler-Lagrange equations (69) yield 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 ϕ j d,e,f , ϕ j d+1,e,f , ϕ j d,e+1,f , ϕ j d+1,e+1,f and the elastic body nodes are defined as follows, see Fig. 9, The expression (ϕ j a−1,b,c − ϕ j a,b,c ) × (ϕ j a,b−1,c − ϕ j a,b,c ) represents the outward pointing normal vector field to the body, similarly in the other three conditions. The associated penalty terms Φ 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., similarly for Ψ im 2 , Ψ im 3 , Ψ im 4 , see Appendix A.4 for the expressions of these derivatives.
The resulting discrete action functional is 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 ϕ j d,e,f , the nodes of the elastic body in the four terms i = 1, ..., 4 are chosen in accordance with (71) for each functions Φ im i , i = 1, ..., 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 with similarly for the fluid.

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 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 ρ 0 = 997 kg/m 2 , γ = 6, A =Ãρ −γ 0 with A = 3.041 × 10 4 Pa, and B = 3.0397 × 10 4 Pa. The size of the discrete reference configuration at time t 0 is 0.3m×0.3m×0.05m, with space-steps ∆s 1 = ∆s 2 = 0.05, ∆s 3 = 0.025m.
3. The four upper edges of the Mooney-Rivlin incompressible model are fixed. 4. The time-step is ∆t = 2 × 10 −4 and the test is carried out for 0.4s.
The results are reproduced in Fig. 10, where we used the notation J d = (J 1 , J 2 , J 3 , J 4 , J 5 , J 6 ) ∈ se(3) * . Due to the presence of the gravitation potential, the discrete Lagrangian is only invariant under the subgroup 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 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.

Convergence tests for the isotropic nonlinear hyperelastic model
Consider the Mooney-Rivlin incompressible model described above, with properties ρ 0 = 945 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.4m × 0.4m × 0.2m. 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 ∆t and ∆s i , i = 1, 2, 3. 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).

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.1 2D discrete hyperelastic models
To compute the partial derivatives D k L j a,b of the discrete Lagrangian (34)  Tr S ( j a,b )F ( j a,b ) T D k F ( j a,b ) · δϕ for all = 1, ..., 4 and all k = 1, 3, 5, 7. We observe that the 16 expressions D k F ( j a,b ) are independent of the value of the field ϕ d and on the location of the spacetime node (if the spacing ∆s a and ∆s b are constant), i.e., we can write D k F ( j a,b ) = D k F . For instance, for = 1, we have and for = 2, we have D 1 F 2 · δϕ = 0, 1 ∆s a δϕ D 3 F 2 · δϕ = − 1 ∆s b δϕ, − 1 ∆s a δϕ D 7 F 2 · δϕ = 1 ∆s b δϕ, 0 D 5 F 2 · δϕ = [0, 0] , similarly for = 3, 4. We can thus write where in the last equality we defined the operator D k , k = 1, 3, 5, 7, = 1, 2, 3, 4, by duality. From this, the partial derivatives D k L j a,b of the discrete Lagrangian can be computed. For instance for k = 1 one has 1 vol( ) ρ e 0 ∇Π(ϕ j a,b ).
Using this result, one then derives (42)