# A mixed variational principle in nonlinear elasticity using Cartan's moving frame and implementation with finite element exterior calculus

This article offers a new perspective for the mechanics of solids using moving Cartan's frame, specifically discussing a mixed variational principle in non-linear elasticity. We treat quantities defined on the co-tangent bundles of reference and deformed configurations as additional unknowns. Such a treatment invites compatibility of the fields with base-space (configurations of the body), so that the configuration can be realised as a subset of the Euclidean space. We first rewrite the metric and connection using differential forms, which are further utilised to write the deformation gradient and Cauchy-Green deformation tensor in terms of frame and co-frame fields. The geometric understanding of stress as a co-vector valued 2-form fits squarely within our overall program. We show that, for a hyperelastic solid, an equation similar to the Doyle-Erciksen formula may be written for the co-vector part of stress. Using these, we write a mixed energy functional in terms of differential forms, whose extremum leads to the compatibility of deformation, constitutive rules and equations of equilibrium. Finite element exterior calculus is then utilised to construct a finite dimensional approximation for the differential forms appearing in the variational principle. These approximations are then used to construct a discrete functional which is then numerically extremised. This discertization leads to a mixed method as it uses independent approximations for differential forms related to stress and deformation gradient. The mixed variational principle is then specialized for 2D case, whose discrete approximation is applied to problems in nonlinear elasticity. An important feature of our FE technique is the lack of additional stabilization. From the numerical study, it is found that the present discretization also does not suffer form locking and related convergence issues.

## Authors

• 3 publications
• 2 publications
• 4 publications
• 5 publications
08/31/2021

### A mixed method for 3D nonlinear elasticity using finite element exterior calculus

This article discusses a mixed FE technique for 3D nonlinear elasticity ...
02/03/2022

### A novel four-field mixed variational approach to Kirchhoff rods implemented with finite element exterior calculus

A four-field mixed variational principle is proposed for large deformati...
05/05/2022

### Geometrical Modelling and Numerical Analysis of Dislocaion Mechanics

This study undertakes the mathematical modelling and numerical analysis ...
11/15/2020

### Quantitative towers in finite difference calculus approximating the continuum

Multivector fields and differential forms at the continuum level have re...
08/20/2020

### Modeling flexoelectricity in soft dielectrics at finite deformation

This paper develops the equilibrium equations describing the flexoelectr...
08/22/2021

### Fast MATLAB evaluation of nonlinear energies using FEM in 2D and 3D: nodal elements

Nonlinear energy functionals appearing in the calculus of variations can...
09/08/2020

### Three-field mixed finite element methods for nonlinear elasticity

In this paper, we extend the tangential-displacement normal-normal-stres...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Mixed or complementary variational principles have a long history in solid mechanics [22]. These variational principles are of significant usefulness in developing approximation techniques where conventional or single field approximations perform poorly [8, 31]. The Hu-Washizu (HW) principle, for instance, is a three field variational approach commonly used to construct FE approximations in nonlinear solid mechanics. The HW variational principle takes deformation gradient and stress as additional inputs along with deformation. From a geometric perspective, deformation gradient and stress are infinitesimal (tangent space based) quantities whose origins can be traced to the geometric hypothesis placed on the configurations of the body. However, scarce little has been achieved in understanding these principles from a geometric standpoint. Electromagnetism is a classical example of a field theory which was reformulated using ideas from geometry. The theory of differential forms [14] played a pivotal role in this reformulation, leading to a deeper understanding of the field equations and novel numerical techniques to approximate their solution [6]. Later, within the broad program of geometrisation of physics, differential forms became an indispensable tool in the mathematical description of physical phenomena. Even though linear elasticity has a lot in common with the theory of electromagnetism [11], the scenario is quite different for nonlinear elasticity. While electric and magnetic fields are the quantities of interest in electromagnetism, the interest in elasticity is on stresses and strains. The distinction however, arises from the nature of the base-spaces (configurations). In the case of electromagnetism, the base-space is fixed and only fields (like electric and magnetic fields) have an evolution. In contrast, for an elastic solid, the base manifold (or configuration) evolves with the deformation process. This not-so-subtle difference prevents one from mimicking the techniques used to geometrically formulate electromagnetism in the context of nonlinear elasticity. Moreover, while electromagnetism can be formulated entirely using differential forms, nonlinear elasticity requires both vector fields and differential forms for a complete description.

The method of moving frames developed by Cartan [7, 13]

is an effective tool that works seamlessly with both vector fields and differential forms. Built on the theory of exterior calculus, Cartan’s moving frames can also encode the connection information on a manifold, which is indispensible for differentiating vector fields on a manifold. To each point of the configuration, the method of moving frames assign a bunch of orthonormal vectors called frame fields. The rate at which these frame fields vary across the configuration defines the connection 1-forms. These connection 1-forms must however satisfy the structure equations so that the parallel transport they encode conforms to the underlying manifold structure, which will presently be assumed Euclidean. Apparently, attaching a set of vectors to a material point is nothing new in continuum mechanics. Many such models have been put forth, starting from Cosserat to micro-morphic theories; these theories are sometimes referred to as micro-continuum theories

[10]

. None of them however encode the connection information of the deforming body whilst evolving the frames or the directors. For these models, directors are just additional degrees of freedom to hold energy. A major difficulty with this point of view is that it does not clarify the geometry within which the model is working. An immediate consequence is that it is impossible to give a co-ordinate independent meaning to the derivatives appearing in the equations of motion. A similar scenario holds, for instance, with shell theories, including the computational schemes

[32] which are built on models that include directors as degrees of freedom in addition to the mid-surface deformation. Clearly, these director degrees of freedom are proxies for the connection induced on the surface. However, since the directors do not satisfy Cartan’s structure equations, it may not be possible to realise the deformed shell as subsets of an Euclidean space, let alone the more general case of shells with intrinsic curvature (i.e. defects).

An attempt at utilizing Cartan’s moving frames to formulate the equations of elasticity was made by Frankel [12]; however his efforts went largely unnoticed and fell short of offering an appropriate computational implementation. In this work, we formulate the kinematics of an elastic body in the language of differential forms. The advantages of having the kinematics formulated in terms of moving frames are twofold. First, important kinematic quantities are described using differential forms that explicate on issues related to compatibility. Second, the geometric hypothesis behind the kinematics is made explicit. The hypothesis that the geometry of non-linear elasticity is Euclidean [18] conforms well with the tensor fields that typically describe the local state of deformation – the right Cauchy-Green deformation tensor to wit, whose roots can be traced to the (Euclidean) metric tensor. Moreover, compatibility equations for both the deformation gradient and Cauchy-Green deformation tensor [35] depend on the geometry of the configuration. Compatibility in terms of the deformation gradient is related to the vanishing of the torsion tensor while compatibility in terms of the Cauchy-Green deformation tensor is related to the vanishing of the curvature tensor [5]

. Both notions of compatibility are however related to the affine connection placed on the configuration. The continuum mechanical definition of stress also has a geometric meaning; Frankel

[12] describes Cauchy stress as a bundle valued differential form. The basic idea in his construction is to decompose the stress tensor into a traction (form) and an area component (form). Such a program was later pursued by Segeve and Rodnay [28] and Kanso et al. [17]. However, both Segeve and Rodnay and Kanso et al. did not pursue a variational principle using this description. From the decomposition of Cauchy stress, it is clear that the constitutive rule needs to be written only for the traction component since the area component is determined kinematically.

Conventionally, numerical techniques for nonlinear elasticity were largely based on single field approximations. It was soon realised that these methods suffered from numerical instabilities, thus affecting the convergence. Displacement based methods with additional stabilization, which barely carries any physical meaning, is a common technique to circumvent these difficulties. Methods like assumed strain and enhanced strain techniques belong to this class; these methods introduce addition terms in the energy function whose origin is purely numerical. Within nonlinear elasticity, techniques based on mixed FE methods are already the preferred choice [29, 30, 1] for large deformation problems since they avoid numerical instabilities like locking and preserve important conserved quantities. However, constructing stable mixed finite elements is difficult even in the case of linear elasticity. A few researchers have turned to ideas from differential geometry to construct stable, well performing mixed-FE techniques for nonlinear elasticity; we cite Yavari [34] as an example of one such attempt where tools from discrete differential forms [16] were utilised to construct stable discretisation schemes. These schemes try to preserve important properties of differential forms even in the discrete setting. And yet, there is barely an attempt at developing a variational principle in solid mechanics to unify the posing and numerical solution of a geometrically conceived model.

Finite element exterior calculus (FEEC) [3] is an FE technique developed by Arnold and co-workers to unify finite elements like Raviart-Thomas, Nédélec [24, 20, 21] and other carefully handcrafted elements under a common umbrella. FEEC relies on the theory of differential forms to achieve this unification [6, 3, 2, 15]. The algebraic and geometric structures brought forth by differential forms are instrumental in achieving this. Recently, Angoshtari et al. [1] and Shojaei and Yavari [29] have brought on ideas from algebraic topology to discretize the equations of nonlinear elasticity. These authors have constructed mixed FE techniques using HW variational principle. Their discretization was based on a differential complex which unfortunately is not the same as the de Rham complex. These methods still require stabilization terms in the three dimensional case [30]. Moreover, from the description of the complex given in Angoshtari et al., it is not clear how the HW variational principle is related to the complex and how the operators defining the complex are affected by the connection placed on the configuration.

The goals of this article are twofold. The first is to develop a mixed variational principle for nonlinear elasticity that takes differential forms as its input argument. Towards this, we first reformulate the kinematics and kinetics of an elastic solid using differential forms. The kinematics of deformation is laid out via Cartan’s method of moving frames. The structure equations associated with the moving frames establish the important relationship between the geometric hypothesis of the configuration and measures of deformation. The geometric understanding of stress as a bundle valued differential form is then exploited to write the kinetics in terms of differential forms. These two ideas are then used to rewrite the conventional HW variational principle, which now has a bunch of differential forms and deformation as its input arguments. The proposed mixed functional is then extremised with respect to these differential forms to arrive at the equations of mechanical equilibrium, constitutive relation and compatibility constraints. The second goal of this work is to use FEEC to discretize the proposed mixed variational principle. Towards this end, the spaces and are used to discretize the differential forms describing the kinetics and kinematics. Using these approximations, a discrete mixed functional is arrived at, which can then be extremized using numerical optimization techniques.

The rest of the article is organized as follows. A brief introduction to Cartan’s moving frames and the associated structure equations are presented in Section 2. The kinematics of an elastic body is then reformulated in Section 3, using the idea of moving frames. In this section, important kinematic quantities like deformation gradient and right Cauchy-Green deformation tensor are described using frame and co-frame fields. This section also contains a discussion on affine-connections using connection 1-forms. In Section 4, we introduce stress as a co-vector valued differential 2-form; this interpretation is originally due to Frankel [12]. We then derive a relationship similar to the Doyle-Ericksen formula relating the stored energy density function to the traction form. In Section 5, we rewrite the mixed variational principle in terms of differential forms describing the kinematics and kinetics of motion. We then show that variation of the mixed functional with respect to different input arguments leads to the compatibility of deformation, constitutive rule and equations of equilibrium. We also remark on the interpretation of stress as a Lagrange multiplier enforcing compatibility of deformation. Section 6 discusses a discrete approximation of differential forms on a simplicial manifold. While these ideas have their roots in the work of Whitney [33], we adopt a description of the FEEC within the framework of Cartan’s moving frames. These are utilized to construct a discrete approximation to the mixed variational principle, which is numerically extremized using Newton’s method. This FE approximation is then applied to standard benchmark problems in nonlinear elasticity in order to assess the performance of the numerical technique against instabilities like volume and bending locking. Finally, in Section 9, we discuss on the usefulness of the moving frames in formulating other theories in nonlinear solid mechanics (like Kirchhoff shells and dislocation mechanics) and the extension of the present numerical techniques to such nonlinear theories.

### 1.1 A remark on notations

We do not use bold face letters to distinguish between a scalar, a vector or a tensor. Indices are used to index objects of the same kind and not the components; for example, if we have three forms, we may denote the th form by . A symbol with one index does not mean that it is a component of a vector or a -form. The objects featured in the theory are defined wherever they appear first. Often in nonlinear elasticity, lower and upper case indices are used to distinguish objects in the reference and deformed configurations. We do not follow this convention since we adopt separate notations for the same object defined in the reference and deformed configurations. We follow the convention of Einstein summation over repeated indices. In places where the summation convention is not followed, we make it explicit.

## 2 Cartan’s moving frame

In this section, we present a brief introduction to the method of moving frames; our motive being to review some basic results so that the kinematics of an elastic solid can be written in terms of moving frames. For a detailed exposition on moving frames, the reader may consult [7] and [13]. Cartan introduced the method of moving frames as a tool to study the geometry of surfaces. These techniques were later extended to study the geometry of abstract manifolds. Common examples of moving frames include the Frenet frame for a curve and Durboux frame for a surface.

We denote the reference and deformed configurations of a body by and ; the respective tangent bundles are denoted by and . Both and are smooth manifolds with boundaries; their boundaries are denoted by and . These configurations are endowed with a chart from which they inherit their smoothness. Positions (placements) of a material point in the reference and deformed configurations are denoted by and respectively.

At each tangent space of a configuration, we choose a collection of orthogonal vectors, which we call the frame. The orthogonality of the frame field is with respect to the Euclidean inner product of the respective configuration. The frame fields of the reference and deformed configurations are denoted by and respectively. A collection of frame fields that span the tangent spaces constitutes a moving frame or simply a frame (allowing for a slight abuse of the terminology). We denote frames for the reference and deformed configurations by and . Given a frame for a tangent bundle, the natural (algebraic) duality between tangent and co-tangent spaces induces a co-frame for the co-tangent bundle as well. These co-frames (at a point) constitute a basis for the co-tangent spaces of the respective configurations. We denote the co-frames of the reference and deformed configurations by and respectively, where and are sections from the cotangent bundles of the reference and deformed configurations. The natural duality between frame and co-frame fields of the reference and deformed configurations may be written as,

 Ei(Ej)=δij;ei(ej)=δij;Ei∈T∗B,ei∈T∗S (1)

For a material point in the reference configuration, the differential of position is denoted by . In terms of the frame and co-frame fields, it can be written as,

 dX=Ei⊗Ei (2)

Similarly, in terms of the frame and co-frame fields in the deformed configuration, the differential of position is given by,

 dx=ei⊗ei (3)

From the definition of (2), it follows that a tangent vectors from the reference configuration is mapped to itself under . To see this, choose with . Substituting the latter and using the definition of , we arrive at . Using the duality between the frame and co-frame fields, we conclude that . Similarly, maps a tangent vector from the deformed configuration to itself. The differential of a position vector in the reference or deformed configuration is thus an identity map on the corresponding tangent space.

Similar to the differential of position, one may also define the differential of a frame, which in the reference configuration is given by,

 dEi=γji⊗Ej (4)

where,

is called the connection matrix; it contains 1-forms as its entries. Because of the orthogonality between the frame fields, the connection matrix is skew symmetric, i.e.

. Similarly, the differential of a frame for the deformed configuration is given by,

 dei=¯ωji⊗ej (5)

is the connection matrix of the deformed frame fields. It is also skew, i.e. .

For a given choice of connection 1-forms and co-frame fields, there are certain compatibility conditions (Poincaré relations) which guarantee the existence of the placements and . These equations are called Cartan’s structure equations. In the present context (of all manifolds being Euclidean), the first compatibility condition establishes the torsion free nature of the configuration. For the reference and deformed configurations, this condition may be written as,

 d2X=0;d2x=0 (6)

Plugging (2) and (3) into the above equation and making use of (4) and (5) leads to,

 dEi=γij∧Ej;dei=¯ωij∧ej (7)

The second compatibility condition presently establishes that the reference and deformed configurations are curvature-free. This leads to the following conditions on the reference and deformed frame fields,

 d2Ei=0;d2ei=0 (8)

Using the differential of the frame fields in the above equations yields,

 dγij=γik∧γkj;d¯ωij=¯ωik∧¯ωkj (9)

For a simply connected body, the structure equations (7) and (9) (for the reference and deformed configurations) provide the necessary kinematic closure to ensure that the configurations can be embedded with an Euclidean space. Indeed, without this closure effected by the structure equations, a model cannot in general produce a deformed configuration which is a subset of an Euclidean space.

### 2.1 Differentials of position and frame

We now present the geometric meaning of the infinitesimal quantities (differentials of position and frame) introduced in the last subsection. Consider the differential of the position for the reference configuration given in (2). For a given coordinate system, the position is a smooth function of its coordinates . Let be a curve parametrized by its arc length, , , where denotes the metric tensor of the reference configuration. The frame fields can be constructed by an orthonormalisation (Gram-Schmidt procedure) of the tangent vectors to the coordinate curves. Fig. 1 shows the coordinate curves and frame at and . The tangent vector to the curve is written as , where are real numbers. Now, the differential of position, which is a vector valued 1-form, can be integrated along the curve to produce a vector; this vector translates the position to . This translation may be formally written as,

 X(b)−X(a) =∫badX(ciEi)ds =∫baciEjEj(Ei)ds, =∫baciEjδjids, =∫baciEids. (10)

In the last equation, and can vary along the curve . The above interpretation of is very similar to that of 1-forms as real numbers defined on curves.

We now consider the differential of frame. Integrating (4) along , we have,

Evaluating

on the tangent vector produces a skew symmetric matrix with real co-efficients. In other words, the above integration can be written as a solution to the ordinary differential equation,

 ˙Ei=γjiEj, (11)

may be understood as the derivative with respect to the parameter . Alternatively, one can understand 11, as a restriction of (4) to a curve. Given the initial , the solution to (11) is a rotation matrix which relates the frames between two points on the curve . This is formally written as,

 Ei(s)=R(s)Ei(a). (12)

Figure 1 (b) depicts this idea by overlaying the frames at and . From this discussion, it is clear that the vector translating the point to is dependent on the frame, the co-frame and the curve chosen for integration. However, the vector is path and frame independent since the body is a subset of an Euclidean space. This path independence is exactly what the structure equations enforce.

### 2.2 Affine connection via frame fields

We now discuss the affine connection and covariant differentiation encoded by the connection 1-forms discussed in the previous subsection. An affine connection on a smooth manifold is a device used to differentiate sections of vector and tensor bundles in a co-ordinate independent manner. Let be an arbitrary section from . The covariant derivative of in the direction of is given by,

 ∇eiw=dwj(ei)ej+wj¯ωkj(ei)ek (13)

It is easy to check that the above definition is a differential satisfying the properties,

 ∇feiw =f∇eiw;f∈Λ0 ∇ei(w+v) =∇eiw+∇eiv;w,v∈TS (14) ∇ei(fw) =ei[f]w+f∇eiw

Using these properties, it is possible to extend the above definition of covariant differentiation to arbitrary tensor fields [18].

## 3 Kinematics

The deformation map sends the placement of material points in the reference configuration to their corresponding placements in the deformed configuration. We denote the deformation map by so that . The differential of the deformation or the deformation gradient, denoted by , maps the tangent space of the reference configuration to the corresponding tangent space in the deformed configuration. For an assumed frame field (for both reference and deformed configurations), the differential of deformation can be obtained by pulling the co-vector part of the deformed configuration’s differential of position back to the reference configuration.

 dφ =ei⊗φ∗(ei) =ei⊗θi;θi∈T∗B,ei∈FS (15)

In writing (15), we have introduced the following definition: . In our construction, the 1-forms contains local information about the deformation map . The 1-forms are a primitive variable in our theory; we call these differential forms, the deformation 1-forms. From (15), we see that the vector leg of the deformation gradient is from the deformed configuration, while the co-vector leg is from the reference configuration, making it a two-point tensor [18]. The action of the deformation gradient on a vector is given by,

 dφ(V)=eiθi(V) (16)

Since are real numbers, the above equation is a linear combination of tangent vectors from the deformed configuration. Conventionally, deformation gradient is introduced as the differential of deformation; we have not taken this perspective since our interest is in constructing mixed variational principles for nonlinear elasticity. Another important aspect of the above construction is that, we have written the deformation gradient using a frame and a co-frame. Contrast this with a conventional mixed method, where the deformation gradient is identified with its components in a particular coordinate system.

### 3.1 Strain and deformation measures

The notion of length is central to continuum mechanics; important kinematic quantities like strain and rate of deformation are derived from it. Indeed, it may not be possible to assess the state of deformation without a metric structure (notions of length and angle) for both reference and deformed configurations. The notion of length is encoded by a symmetric and positive definite tensor, defined on the tangent space of the respective configuration. The metric tensor of the reference and deformed configurations are denoted by and respectively; and . In this work, we assume the metric structures of both reference and deformed configurations to be Euclidean. In terms of the co-frame field, the metric tensor of the reference configuration is given as,

 G =dX.dX =(Ej⊗Ej).(Ei⊗Ei) (17) =Ei⊗Ei. (18)

The dot product introduced in the above equation is the inner product between the vector legs of , which is computed using the Euclidean inner product. Similarly, the metric tensor in the deformed configurations may be written as,

 g =dx.dx =(ei⊗ej).(ej⊗ei) (19) =ei⊗ei. (20)

In terms of the frame fields, the inverses of the metric tensors for the reference and deformed configurations are written as,

 G−1=Ei⊗Ei;g−1=ei⊗ei. (21)

Now, the right Cauchy-Green deformation tensor may be obtained as the pull-back of the deformed configuration’s metric tensor. In terms of the deformation 1-forms, this relationship may be written as,

 C =φ∗(g) =φ∗(ei⊗ei) =θi⊗θi. (22)

An alternative way to compute the is to use the usual definition in continuum mechanics, . Here, is understood to be the adjoint map induced by the metric structure. Using the orthonormality of the frame field we arrive at,

 C =(θi⊗ei)(ej⊗θj) =θi⊗θi. (23)

The calculations leading to (22) and (23) are exactly the same; only the sequence in which pull-back and inner product are computed differs. The Green-Lagrangian strain tensor may now be written as,

 E =12(C−G) =12[(θi⊗θi)−(Ei⊗Ei)]. (24)

The first invariant of the right Cauchy-Green tensor is given by,

 I1=⟨θi,θi⟩G. (25)

Here denotes the inner product induced by . The area forms induced by the co-frame of the reference configuration are given by,

 A1=E2∧E3;A2=E3∧E1;A3=E1∧E2, (26)

Similarly, the area-forms induced by the co-frame of the deformed configuration are given by,

 a1=e2∧e3;a2=e3∧e1;a3=e1∧e2, (27)

These area forms and serve as a basis for the space of 2-forms defined on their respective configurations. The area forms in the deformed configuration may be pulled back to the reference configuration under the deformation map. These pulled-back area forms are denoted by . In terms of the deformation 1-forms, the pulled back area forms can be written as,

 A1=θ2∧θ3;A2=θ3∧θ1;A3=θ1∧θ2. (28)

In terms of the pulled-back area forms, the second invariant of may now be written as,

 I2=⟨Ai,Ai⟩G (29)

In terms of the co-frame fields, the volume forms of reference and deformed configurations may be written as,

 V=E1∧E2∧E3;v=e1∧e2∧e3. (30)

The pull back of the volume form in the deformed configuration to the reference configuration is denoted by . In terms of the deformation 1-forms, may be written as,

 V=θ1∧θ2∧θ3 (31)

Finally, in terms of the pulled-back volume-form, the third invariant of is given by,

 I3=(⋆V)2 (32)

In the above equation denotes the Hodge star operator, which establishes an isomorphism between the space of forms and forms. We also define or simply .

## 4 Stress as co-vector valued two-form

In the previous section, we reformulated the deformation gradient and right Cauchy-Green deformation tensor in terms of differential forms. We now present a geometric approach to stress, originally due to Frankel [12] and subsequently developed by Segeve and Guy [28] and Kanso et al. [17]. Even though this approach is intuitive and geometric, it was never used to construct a variational principle. This geometric approach to stress has its origins in classical dynamics [19], where force is understood as a co-vector. Identifying force with a co-vector permits us to write power as a pairing between force and velocity (action of a 1-form on a vector) without the use of a metric tensor. Extending this concept to stress in continuum mechanics is non-trivial and requires the machinery of bundle valued differential forms [17]. As in classical dynamics, this interpretation of stress as a bundle valued differential form permits us to write the power expended by stress upon deformation without using a metric (i.e. without explicitly bringing into play the geometric structure of the configuration).

We denote the Cauchy stress tensor by . The traction acting on an infinitesimal area element with unit normal is denoted by , which is given by the well known formula . The traction is a force which depends on the material point and the area sustaining it. The relation between normal and traction is postulated to be linear. In the language of differential forms, an infinitesimal area is regarded as a form, while from classical dynamics we also know that force is a co-vector or a form. Putting these two ideas together, we are led to a geometric definition of Cauchy stress given by,

 σ=ti⊗ai; (33)

Recall that is the area form of the deformed configuration, which sustains the traction vector . The tensor product in the above equation is due to the linearity between traction and area forms. From this equation, it is easy to see that the area form changes orientation if the order of the co-vectors in the area form is reversed. Geometrically, if there are linearly independent area forms on the manifold, the stress tensor assigns to each area form a 1-form called traction. With this understanding, Cauchy stress may now be identified with a section from the tensor bundle which has the deformed configuration as its base space.

In contemporary continuum mechanics, Nanson’s formula describes the transformation of an infinitesimal area under the deformation map [23]. Geometrically, Nanson’s formula is nothing but the pull-back of an area form in the deformed configuration under the deformation map. These pulled-back area forms are given in (28).

The first Piola stress may now be obtained by pulling back the area leg of the Cauchy stress under the deformation map. This partial pull-back of the Cauchy stress is termed the Piola transform. This relationship may be formally written as,

 P =ti⊗φ∗(ai) =ti⊗Ai (34)

Note that in the definition of Piola transform, traction 1-form is left untouched. Thus, contrary to convention, Cauchy and first Piola stresses are now identified as third order tensors, not the usual second order. This ambiguity can be removed if one applies the Hodge star on the area leg of these two stresses. In three dimensions, the Hodge star in question establishes an isomorphism between differential forms of degree 2 and 1. The usual definition of stress may thus be recovered as,

 σ =ti⊗⋆(ai) (35) P =ti⊗⋆(Ai). (36)

Kanso et al. made a distinction between the stress tensors given in (33), (34) and (35), (36). However we do not see a need for it, since both the usual and geometric definitions of stress contain exactly the same information; only the ranks of these tensors are different.

### 4.1 Traction 1-form via stored energy function

The Doyle-Ericksen formula is an important result in continuum mechanics [9], which relates the Cauchy stress and the metric tensor of the deformed configuration. For a stored energy density function , Doyle-Ericksen formula gives us the following relationship,

 σ=2∂W∂g. (37)

In writing (37), we have assumed that is frame-invariant. From the discussion presented so far, it is seen that the area leg of the Cauchy stress tensor is determined by the choice of coordinate system (frame and co-frame fields) for the tangent bundle of the deformed configuration. On the other hand, the area leg of the first Piola stress is determined by both the deformation map and the co-ordinate system for the tangent bundle of the deformed configuration. Clearly, the area leg of a stress tensor does not require a constitutive rule; it is only the traction component that demands a constitutive rule.

We now claim that for a stored energy function , the traction form has the following constitutive rule,

 ti=1J∂W∂ei (38)

The last equation is in the same spirit as the Doyle-Ericksen formula. To establish the result in (38), we first compute the directional derivative of along ,

 ∂W∂ei =∂W∂(dφ)∂(dφ)∂ei (39a) =(tj⊗⋆Ak)(ei⊗ej⊗θk) (39b) =⟨⋆Aj,θj⟩ti (39c) =(⋆V)ti (39d) =Jti (39e)

We used a chain rule to arrive at the right hand side of (

39a). In (39b), the expression for Piola stress (as a two tensor) in terms of and the directional derivative of along are used to get the right hand side and performing the required contractions lead to (39c). The claim is finally established by using the definitions of pull-back and Hodge star for volume forms. In three dimensions, constitutive relations have to be supplied to the three traction forms. From these calculation, it is found that the traction 1-forms are conjugate to the frame fields.

If the Cauchy stress (the usual definition) generated by a stored energy function is known, then the expression for the traction 1-form can be computed using the simple relation , where the vector fields are chosen to be elements from the frame of the deformed configuration.

## 5 Mixed variational principle

We first present the conventional HW variational principle for a finitely deforming elastic body. As mentioned, the HW variational principle takes the deformation gradient, first Piola stress and deformation map as input arguments. In the reference configuration, the HW functional for a non-linear elastic solid can be written as,

 IHW=∫B[W(C)−P:(F−dφ)]dV−∫∂B⟨t,φ⟩dA. (40)

In the above equation, is the traction defined on the surface , whose unit normal is . The integration in the above equation is with respect to the volume and area forms of the reference configuration. In (40), the deformation gradient is assumed to be independent; this tensor field is denoted by . On the other hand, the deformation gradient computed as the differential of deformation is denoted by . It is worthwhile to note that the second term in the above equation is bilinear in the Piola stress and the deformation gradient. The variation of the HW functional with respect to deformation, deformation gradient and first Piola stress leads to the equilibrium equation, constitutive rule and compatibility of deformation gradient. This form of HW variational principle has been previously exploited to formulate numerical solution procedures for non-linear problems in elasticity; see[1, 30, 29].

### 5.1 Mixed variational principle with geometric definitions of stress and deformation

We now use the definitions of Cauchy and Piola stresses given in (33) and (34) respectively to rewrite the HW variational principle such that it takes deformation forms, traction forms and deformation map as inputs. We also assume that compatible frames for the reference and deformed configurations are given. This assumption permits us to eliminate the frame fields from the list of unknowns. The mixed functional may be now written as,

 I(θi,ti,φ)=∫BW(θi)dV−(ti⊗Ai)˙∧(ei⊗(θi−dφi))−∫∂B⟨t♯,φ⟩dA. (41)

In (41), denotes a bilinear map. For , and , , the action of this map is given by . Note that the definition of given here is a little different from the one in Kanso et al. [17]. Specifically, we do not use the metric tensor. From our definition of , it is seen that the work done by stress on deformation is metric independent. This property of our current variational formulation brings the continuum mechanical definition of stress a step closer to the definition of force (as a form) in classical mechanics [19]. We did not write the volume form in the second term on the RHS of (41), since the outcome of is a 3-form which can be integrated over the reference configuration to produce work done by traction 1-forms on deformation. Also note that the second term on the RHS in (41) is equivalent to the second term in (40); however now the relationship between the different arguments is multi-linear. The functional given in (41) can also be discussed within the framework given by Oden and Reddy [22] for the construction of complementary variational principles.

Remark 1: In writing (41), we have postulated that the geometry of the body is Euclidean and it is frozen during the deformation process. Indeed, within the present set-up, this assumption can be relaxed by permitting non-integrability in the connection and deformation 1-forms (i.e. by incorporating source terms in the structure equations).

Remark 2: For a frame to represent Euclidean geometry, it is not required that the connection 1-forms be identically zero. It is only required that the structure equations have a zero source term.

We now proceed to obtain the Euler-Lagrange equations or the condition that determines the critical points of the functional . We use the Gateaux derivative for this purpose. Let denote a small parameter and the direction in which the change in the functional is computed; this change is often referred to as the variation of .

We first calculate the variation of with respect to traction 1-forms; , where are assumed to be from the tangent space of . The perturbed functional in the direction of can be written as,

 I(ϵ)=∫BW(θi)dV−((ti+ϵ^ti)⊗Ai)˙∧(ej⊗(θj−dφj))−∫∂B⟨t♯,φ⟩dA. (42)

Using the definition of and Gateaux derivative, we get a vector valued 3-form for each . These three 3-forms have to be equated to zero to get the condition for critical points in the direction of traction 1-forms. These conditions may be formally written as,

 (43)

Since are orthonormal with respect to a positive definite metric, the above equation can be true only when the coefficient matrix on the LHS is zero, which leads to the following conditions,

 (A1∧(θ1−dφ1))=0;(A1∧dφ2)=0;(A1∧% dφ2)=0 (44a) (A2∧dφ1)=0;(A2∧(θ2−dφ2))=0;(A2∧% dφ2)=0 (44b) (A3∧dφ1)=0;(A3∧dφ2)=0;(A3∧(θ3−d% φ3))=0 (44c)

Using the definition of , the above equations may be recast as,

 ⎡⎢⎣θ1−dφ1dφ2d% φ3θ1θ2−dφ2dφ3θ1dθ2dθ3−φ3⎤⎥⎦∧⎡⎢⎣θ2∧θ3θ3∧θ1θ1∧θ2⎤⎥⎦=⎡⎢⎣000⎤⎥⎦. (45)

For these equations to hold, the following conditions must be met,

 θ1−dφ1=0;θ2−dφ2=0;θ3−dφ3=0. (46)

The above condition simply states that there exist three forms whose exterior derivatives are the deformation 1-forms; or in other words, the deformation 1-forms are exact and are the potentials for the corresponding deformation 1-forms.

We now compute the variation of with respect to the deformation forms. Incremental changes in the deformation 1-forms may be written as, , where is assumed to be an element from the tangent space . The perturbed functional in the direction of deformation forms may be written as,

 I(ϵ)=∫BW(θi+ϵ^θi)dV−(ti⊗A(ϵ)i)˙∧(ej⊗θj(ϵ))−∫∂B⟨t♯,φ⟩dA. (47)

Using the definition of the Gateaux derivative, for each we have,

 ∂W∂θ1 =[t1(e1)⋆(θ2∧θ3)−t2(e1)⋆(dφ1∧θ3)+t2(e2)⋆((θ2−dφ2)∧θ3)−t2(e3)⋆(dφ3∧θ3) −t3(e1)⋆(θ2∧dφ1)+t3(e2)⋆(θ2∧dφ3)+t3(e3)⋆(θ2∧(θ3−dφ3))]♯ (48a) ∂W∂θ2 =[t1(e1)⋆(θ3∧(θ1−dφ1))−t1(e2)⋆(θ3∧dφ2)−t1(e3)⋆(θ3∧dφ3)−t2(e2)⋆((θ3∧dθ1) −t3(e1)⋆(dφ1∧θ1)−t3(e2)⋆(θ2∧θ1)+t3(e3)⋆((θ3−dφ3)∧θ1)]♯ (48b) ∂W∂θ3 =[t1(e1)⋆((θ1−dφ1)∧θ2)−t1(e2)⋆(dφ2∧θ2)−t1(e3)⋆(dφ3∧θ2)−t2(e1)⋆(θ1∧dφ1) −t2(e2)⋆(θ1∧(θ2−dφ2))+t2(e3)⋆(θ1∧dφ3)+t3(e3)⋆(θ1∧θ2)]♯. (48c)

If we now take into account the compatibility equations previously established in (46), the last equations reduce to,

 ∂W∂θ1 =[t1(e1)⋆(θ2∧θ3)+t2(e1)⋆(θ3∧θ1)+t3(e1)⋆(θ2∧θ1)]♯ (49a) ∂W∂θ2 =[t1(e2)⋆(θ2∧θ3)+t2(e2)⋆((θ3∧θ1)+t3(e2)⋆(θ1∧θ2)]♯ (49b) ∂W∂θ3 =[t1(e3)⋆(θ2∧θ3)+t2(e3)⋆(θ3∧θ1)+t3(e3)⋆(θ1∧θ2)]♯. (49c)

From these equations, we see that a form accompanies the components of traction forms; this is indeed true since we use a Piola transform to write the constitutive rule in the reference configuration. From a comparison of (48) and (49), we note that the expressions for traction in the former has additional terms. These additional terms may be related to incompatibilities created by the emergence of defects (such as dislocations) as the deformation evolves. In other words, without an explicit imposition of the kinematic closure conditions (compatibility conditions), the deformed body may never be realized as a subset of an Euclidean space. To our understanding, this is perhaps the most significant finding of this article.

Finally, we compute the variation of with respect to deformation; , where belongs to . Using the definition of the superimposed incremental deformation in and upon computing the Gateaux derivative, we have the following equation,

 ∫B(ti⊗Ai)˙∧(ej⊗% d^φi)=0 (50)

To complete the variation, we need to shift the differential from . We first calculate the following,

 d(φktj(ek)Aj)=dφk∧tj(ek)Aj+φkd(tj(ek))∧Aj+φktj(ek)dAj (51)

This equation invites a few comments. The first is that we are calculating the exterior derivative of a form, with being a scalar. Using the product rule of differentiation, we have expanded the right hand side of (51). The second term in (51) should be evaluated using the connection 1-forms since it involves the exterior derivative of a vector. This terms is relevant when one works with a frame field whose connection 1-forms are different from zero. If we invoke the compatibility of deformation, we have , which leaves (51) in the following form,

 d(φktj(ek)Aj)=dφk∧tj(ek)Aj+φkd(tj(ek))∧Aj (52)

An expression similar to (52) was utilized by Kanso et al. [17] to define the mechanical equilibrium. The expression for the exterior derivative defined in (52) involves the connection 1-forms of the manifold, which is similar to the covariant exterior derivatives used in gauge theories of physics [14]. Using (51) in (50) leads to,

 ∫Bd(^φktj(ek)Aj)−^φkd(tj(ek))∧Aj=0 (53)

The first term in the above equation may be converted to a boundary term via Stokes’ theorem leading to,

 ∫∂B^φktj(ek)Aj−∫B^φkd(tj(ek))∧Aj=0 (54)

Using the arbitrariness of , we conclude that,

 d(tj(ek))∧Aj=0 (55)

This is the condition for the critical point of the mixed functional in the direction of deformation, which represents the balance of forces. Note that the connection forms of the deformed configuration appear through the exterior derivatives of the frame fields of the deformed configuration.

### 5.2 Stress a Lagrange multiplier

For a hyper-elastic solid, stress is derived from the stored energy which may be written as a function of the deformation gradient. This assumption permits us to write the equations of equilibrium as the Euler-Lagrange equation of the stored energy functional. In a certain sense, the stress generated in a hyper-elastic solid should satisfy certain integrability condition (i.e. the existence of the stored energy function). Moreover, if we assume the stored energy function to be translation and rotation invariant, it implies equilibrium of forces and moments. Thus for the hyper-elastic solid, balances of forces and moments are consequences of translation and rotation invariance; stress is only a secondary variable introduced for writing the equations of equilibrium in a convenient way.

When formulated as a mixed problem, the stress tensor or more specifically the traction 1-form has a completely different role. Our mixed functional has deformation, deformation forms and stress forms as inputs. For the stored energy, viewed as a function of deformation 1-forms, translation and rotation invariance cannot be discussed directly, since nothing about the geometry of the co-tangent bundle from which the deformation 1-forms were pulled back is known. In other words, there is nothing in the stored energy function that requires the base space for the deformed configuration to be Euclidean. The second term in (41) is introduced to impose this constraint. Observe that, in (41), the second term is multilinear in the input arguments, i.e., stress 2-form, differential of deformation and deformation 1-form. The traction form may now be thought of as a Lagrange multiplier introduced to impose the equality between the differential of deformation and deformation 1-forms. Alternatively, the equality between the differential of deformation and deformation 1-forms implies that the deformed configuration is Euclidean.

## 6 Discretization of differential forms

In this section, we consider the local finite element spaces which are suitable for approximating differential forms over a simplicial complex. A simplicial complex in , denoted by , is a set with simplices as its elements. By an -simplex, we mean the convex hull of points in . We denote a general -simplex by and a specific -simplex by , where ; these are referred to as the vertices of the simplex. We may also call a simplex , a face of . We expect every face of to be an element in and if two faces intersect, then their intersection is also a face in . The dimension of is defined as the largest dimension of the simplex contained in . Finite element mesh created by triangulating a two dimensional domain is a good example of a two dimensional simplicial complex. Such a finite element mesh has nodes, edges and faces, which are simplices of dimensions , and respectively.

From now on we restrict ourselves to two spatial dimensions; techniques discussed here can be extended to any spatial dimensions. The just stated assumption also permits us to work with a simplicial complex of dimension 2. We denote the barycentric coordinates of a 2-simplex (triangle) by , . These coordinates satisfy the relation . In terms of the Cartesian coordinates, the barycentric coordinates of the reference triangle are given by,

 λ0=1−x1−x2;λ1=x1λ2=x2 (56)

where, and are the Cartesian coordinates of a point within the triangle. The reference triangle, showing the vertices and orientations of edges, is presented in Figure 3; we denote this reference triangle by or simply by .

### 6.1 Spaces PrΛn and P−rΛn

We denote the space of variable polynomials of degree by . The space of polynomial differential forms with form degree and polynomial degree is denoted by , . Often we suppress from our notation and simply denote these spaces by and . For vector fields, ,

 PrΛn={ω∈Λn(Rm)|ω(v1,...,vn)∈Pr} (57)

In other words, for a polynomial differential form , the coefficient functions are polynomials of degree . From the definition of the spaces and , their dimensions can be computed as and respectively. For a differential form of degree , the interior product of with a vector , , is given as,

 κvω=ω(v,v1,...,vn−1) (58)

for any vectors, . From the above definition it is easy to see that is a differential form of degree .

For any point , the vector field translates the origin to . Using this vector field , a Koszul type operator on the space of polynomial differential forms with form degree can be defined. This operator is given as,

 κXω=ω(X,v1,...,vn−1) (59)

where, are arbitrary vector fields form . From the definition, it is easy to see that decreases the degree of a differential form by 1 and increases the polynomial degree by 1. An important property of is . The operator also commutes with affine pull-back. If is an affine linear map, , then . The polynomial spaces , can be defined as,

 P−rΛk={ω∈PrΛk|κX(ω)∈PrΛk−1} (60)

The dimension of the space can be computed as which is larger than that of but smaller than that of . In the case of polynomial differential forms with form degree 0, we have ; these spaces can be identified with the Lagrange family of finite element spaces. The spaces and constitute a large family of finite elements. Well known members of this family include the Raviart-Thomas [24] and Nédélec [21] type vector finite elements. These subspaces of polynomial differential forms were proposed by Arnold and co-workers [3] to unify the vector finite elements used in the construction of mixed finite element techniques. These finite dimensional polynomial spaces form the cornerstone for the finite element techniques developed under the umbrella of finite element exterior calculus.

### 6.2 Degrees of freedom and finite element bases

In the previous subsection, we introduced the polynomial spaces and . We now restrict these polynomial spaces to the reference triangle and construct basis functions suitable for computation. The description of the computational basis functions for the spaces and closely follows the work of Arnold et al. [4] where a geometric decomposition was utilized to construct the computational basis function on a simplex. The idea of the geometric decomposition is to index the DoF’s of a finite element (FE) space using the sub-simplices of the simplex on which the FE space is constructed. For Lagrange finite elements over simplices, this amounts to assigning DoF’s to the vertices of the simplex. For finite elements in the family and , DoF’s may be indexed with edges, faces and volume. The following table gives the basis functions for the spaces , and . From Table 6.2, it can be seen that has DoF only on the vertex, while has two DoF’s on each edge and has one DoF per edge.

FE space Node Edge Face
- -
- -
- -

### 6.3 Whitney forms and P−1Λn

Differential forms of degree are functions that take vectors and produce real numbers. Alternatively, one can also view an form as a real number defined on a plane spanned by tangent vectors. The latter definition is more suitable for the construction of discrete differential forms on a simplicial complex. On a simplex, a discrete differential form of form degree can be thought as a real number defined on the subsimplex. Using this idea, Hirani [16] constructs a discrete analog of the exterior calculus. For such a discrete approximation, Whitney [33]

gave a formula for interpolating these differential forms on a

simplex. By producing basis functions whose degrees of freedom are the real numbers defined on subsimplices, Whitney was able to interpolate these differential forms within the simplex. The basis functions, constructed by Whitney to approximate differential forms over a simplex are now-a-days referred to as Whitney forms. In terms of the barycentric coordinates, these polynomial differential forms are given by the formula,

 ϕni=r!∑ϵi∈C(k,r)(−1)iλϵi(dλϵ0∧...∧ˆdλϵi∧...∧dλϵr). (61)

In the above expression