published or not. The documents may come

From fluid and structure dynamical behaviors to vascular pathologies

From fluid and structure dynamical behaviors to vascular pathologies.

Deplano Valérie. CNRS, AMU, ECM, IRPHE UMR7342, 13384 Marseille

Introduction: The aim of my work is to develop biomechanical models to broaden the current knowledge of fluid and structure dynamics which are associated with the function and dysfunction of some vascular systems. This research enables to identify the potential correlations between fluid and structure behaviors and dysfunctions to bring insights on some pathologies. Over the last ten years three complementary approaches have been considered at the macroscopic level: the fluid, the structure and the tissue, allowing to describe their couplings. The works carried out at this scale have highlighted the necessity of assessing the microscopic scale to better understand some mechanisms involved in the genesis and development of pathologies. Numerous works on fluid dynamics have not only been achieved in pathological geometrical models such as arterial or valvular stenoses, arterial grafts, abdominal aortic aneurysms, aortic dissections but also in stents, endovascular and mechanical valvular prostheses. Besides, studies at the microscopic scale related to the behavior of red blood cell concentrated suspensions under flow as well as those linked with the influence of mechanical solicitations on the cellular response that acts on tissue remodeling have also been conducted. Others works on the mechanical behavior of the aortic wall have been carried out at different scales. They investigated the arterial wall microstructure at the fibers level to better understand the mechanisms that drive the macroscopic mechanical response and also to design biomimetic analog materials.

Methods: Several original experimental devices have been developed. 3D in vitro models of blood circulation taking into account the fluid structure interaction and the rheological behaviors have thus been achieved for different pathophysiological arterial segments. Set ups able to perform mechanical characterizations of biological material and tissue have also been realized such as biaxial tensile bench. The associated metrologies are stereo particle image velocimetry, micro PIV, stereo digital image correlation, multiphoton microscopic… In addition, numerical modeling as well as theoretical models have also been developed.

Results and discussion:

After a brief description of some works, the results will be focused on the abdominal aortic aneurysm biomechanics [1, 2, 3], the mechanical characterization of the aortic wall [4,5] and on the behavior of concentrated red blood cell under flow [6]. Discussion will be carry out on the links between biomechanical studies and pathologies understanding. references 1. Deplano, V. et al, (2013), Med. Eng. Phys, 35, p800-809. 2. Deplano, V. et al, (2014), J. of Biomech., 47(6), p1262-1269. 3. Deplano, V. et al, (2016), Exp. Fluids, 57:167. 4. Deplano, V. et al, (2016), J. of Biomech., 49(10), p2031-2037. 5. Bailly et al, (2012), J Mech behaviour Biomed Mater., 10, p151-165. 6. Knapp et al, (2016)

  , 22nd Congress of the ESB, Lyon Bio Valérie Deplano is a Director of Research at CNRS. She received her PhD in Fluids Mechanics at the University of Marseille in 1992. She is an expert in biofluid mechanics, cardiovascular biomechanics, development of multi physics modeling and multi modal in vitro experiments for the study of the fluid structure interactions existing in biological systems. She was in charge of a CNRS national research group (~ 80 persons) on biomechanics and on the biological fluid structure interactions for 8 years. She is currently the leader of the Biomechanics team of IRPHE the 'Institut de Recherche sur les Phénomènes hors Equilibre'.


domain: spi.meca math phys

text:
Jack S Hale

Matteo Brunetti

Stéphane Bordas

Corrado Maurini
email: corrado.maurini@upmc.fr

P A Bordas

Simple and extensible plate and shell finite element models through automatic code generation tools

Keywords: thin structures, plates, shells, finite element methods, domain specific language, FEniCS

published or not. The documents may come

Simple and extensible plate and shell finite element models through automatic code generation tools

Introduction

Plates and shells are solids occupying a spatial domain with one dimension, the thickness, much smaller than the others. This implies the possibility of experiencing large changes of shape even with small material deformations and an approximately linear elastic material behaviour. Their study has received renewed attention in the last decade because of their unique nonlinear behaviour, mainly caused by geometric effects [START_REF] Audoly | Elasticity and Geometry. From hair curls to the non-linear response of shells[END_REF], leading to d-cone singularities and crumpling [START_REF] Conti | Confining thin elastic sheets and folding paper[END_REF], fracture [START_REF] Roman | Fracture path in brittle thin sheets: a unifying review on tearing[END_REF], or multistability [START_REF] Seffen | Morphing bistable orthotropic elliptical shallow shells[END_REF][START_REF] Seffen | Growth and shape control of disks by bending and extension[END_REF]. Examples of modern application fields in engineering include shape control through active materials [START_REF] Klein | Shaping of elastic sheets by prescription of non-Euclidean metrics[END_REF][START_REF] Hamouche | Multi-parameter actuation of a neutrally stable shell: a flexible gear-less motor[END_REF], stretchable electronics [START_REF] Rogers | Materials and Mechanics for Stretchable Electronics[END_REF], soft robotics [START_REF] Lu | Flexible and Stretchable Electronics Paving the Way for Soft Robotics[END_REF], and thin nano-structures e.g. graphene sheets and nanotubes [START_REF] Arroyo | Finite crystal elasticity of carbon nanotubes based on the exponential Cauchy-Born rule[END_REF].

Effective models for plates and shells are two-dimensional, the through-the-thickness kinematics being described by including a suitable microstructure in the model. Their mathematical modelling leads to set of non-linear partial differential equations (PDEs) defined on a two-dimensional manifold embedded in three-dimensional space. The resulting PDEs can be discretised using a variety of numerical techniques, e.g. finite element methods (FEM) [START_REF] Bathe | A four-node plate bending element based on Mindlin/Reissner plate theory and a mixed interpolation[END_REF][START_REF] Nguyen-Thanh | A smoothed finite element method for shell analysis[END_REF], isogeometric analysis [START_REF] Kiendl | Isogeometric shell analysis with Kirchhoff-Love elements[END_REF][START_REF] Benson | Isogeometric shell analysis: The Reissner-Mindlin shell[END_REF], and meshfree methods [START_REF] Hale | A locking-free meshfree method for the simulation of shear-deformable plates based on a mixed variational formulation[END_REF]. FEM-based discretisations of the plate and shell models are available in both commercial (see e.g. Abaqus [2], AceFem [START_REF] Korelc | Multi-language and Multi-environment Generation of Nonlinear Finite Element Codes[END_REF] and Adina [START_REF] Bathe | Finite Element Procedures[END_REF] and Ls-Dyna) and open-source (e.g. GetFem++ [START_REF] Pommier | Getfem++, an open source generic C++ library for finite element methods[END_REF], IgaFem [START_REF] Nguyen | Isogeometric analysis: An overview and computer implementation aspects[END_REF], Code Aster [START_REF] De | Finite element code aster , analysis of structures and thermomechanics for studies and research[END_REF], Mat-Fem [START_REF]Mat-fem: Learning the finite element method with matlab and gid[END_REF] and Elmer [START_REF] Råback | Elmer models manual[END_REF]) software packages. However, formulating the mathematical models and implementing finite element solvers for custom plates and shells models remains a complex, highly technical, and time-consuming task, requiring advanced knowledge in differential geometry, numerical analysis, and mechanics. Shell and plate finite element models are regarded as very advanced topics in structural engineering and applied mathematics curricula, and rarely are graduate students trained in their implementation.

Our work aims to overcome some of the above difficulties by using modern automatic code generation tools and suitable mathematical abstractions. We present easily extensible and customisable methods for the implementation of finite element solvers for various models of thin structures, ranging from linear plates to non-linear shells. We leverage the tools offered by the finite element library FEniCS Project [START_REF] Alnaes | The FEniCS Project Version 1.5[END_REF], in particular its Unified Form Language (UFL) [START_REF] Alnaes | Unified Form Language: A Domain-specific Language for Weak Formulations of Partial Differential Equations[END_REF] and the associated FEniCS Form Compiler (FFC) [START_REF] Kirby | A Compiler for Variational Forms[END_REF], to allow users to straightforwardly formulate complex custom nonlinear shell models with less than one hundred lines of Python code and syntax that closely mirrors the mathematical abstractions of the variational formulation. We exploit the symbolic processing capabilities of UFL to specify the potential energy functionals of various thin structural models in a high-level way and select suitable locking-free finite element spaces for their discretisation. Custom constitutive models are straightforward to implement. Consistent linearisation of the potential energy functional is performed automatically using the symbolic differentiation tools in UFL [START_REF] Alnaes | Unified Form Language: A Domain-specific Language for Weak Formulations of Partial Differential Equations[END_REF]. These symbolic expressions are compiled automatically to finite element code using FFC [START_REF] Kirby | A Compiler for Variational Forms[END_REF]. The differential geometry concepts [START_REF] Ciarlet | Theory of Shells[END_REF] are expressed directly in a very simple and expressive format, which we believe is of great pedagogical value. We mimic the geometrically-exact shell model [START_REF] Ciarlet | Theory of Shells[END_REF], where the reference configuration is described by a mapping from R 2 to a surface embedded in Euclidean space. Although we do not discuss here the implementation of specific preconditioners, our finite element solvers can run without modification on high-performance computing architectures using MPI and PETSc [START_REF] Balay | PETSc Users Manual[END_REF] allowing large-scale calculations to be performed. State-of-the-art linear and non-linear solution strategies are immediately available through PETSc [START_REF] Balay | PETSc Users Manual[END_REF]. These aspects sum up to a unique approach that we believe will be of value to researchers, students, and practitioners working on frontier problems in the mechanics of thin structures. The outcome of our work is distributed in the form of an open-source (LGPLv3) Python library, FEniCS-Shells [START_REF] Hale | FEniCS-Shells[END_REF], collecting an implementation of the models and discretisation techniques presented in this paper and including a large set of examples and documented demos.

We consider here several shearable plate and shell models, including linear plates (Reissner-Mindlin model [START_REF] Reissner | On bending of elastic plates[END_REF][START_REF] Mindlin | Influence of rotatory inertia and shear on flexural motions of isotropic, elastic plates[END_REF]), linear and nonlinear shells (Nagdhi model [START_REF] Naghdi | The theory of shells and plates[END_REF][START_REF] Naghdi | Foundations of elastic shell theory[END_REF]), and weakly non-linear shallow shells (Marguerre-von Kármán model [START_REF] Marguerre | Zur theorie der gekrümmten platte grosser formänderung[END_REF]). For these structural models, the discretisation strategy is particularly important, because shear and membrane locking can lead to an unacceptably slow convergence rate of the finite element solution with the mesh size [START_REF] Chapelle | The Finite Element Analysis of Shells -Fundamentals[END_REF]. After a deep analysis of the methods available in the literature and the possible software tools, we selected two main discretisation techniques: the popular Mixed Interpolation of Tensorial Components (MITC) approach [START_REF] Bathe | A four-node plate bending element based on Mindlin/Reissner plate theory and a mixed interpolation[END_REF][START_REF] Bathe | The MITC7 and MITC9 Plate bending elements[END_REF][START_REF] Durán | On mixed finite element methods for the Reissner-Mindlin plate model[END_REF] and a high polynomial-order Partial Selective Reduced Integration (PSRI) approach [START_REF] Arnold | Locking-free finite element methods for shells[END_REF], both of which can be applied in a uniform manner across user-defined models.

Both the MITC and PSRI approches presented here include some original aspects with respect to the formulations previously introduced in the literature. In order to fit within the abstractions of UFL and the associated code generation tools, MITC is reformulated as a mixed hybridisable variational form with element-wise local projection. The PSRI technique proposed in [START_REF] Arnold | Locking-free finite element methods for shells[END_REF] is extended to nonlinear shells with a modified reduced integration rule and optimised weighting factor for the energy splitting. To provide an efficient parallel implementation of the MITC technique, we also extended the FEniCS automatic assembly tools to include the possibility of eliminating local degrees of freedom through static condensation, a feature that could be potentially ported to the FEniCS library itself in the near future.

The effectiveness of our MITC and PSRI discretisations are shown through a large set of numerical verification tests. Beside many classical examples, we propose a novel numerical verification test based on the exact solution of Mansfield [START_REF] Mansfield | Bending, buckling and curling of a heated thin plate[END_REF] for the thermal buckling of a lenticular plate. We believe that this latter example can usefully complement the traditional gallery of verification tests for plates and shells.

The paper is organised as follows. First, we give an overview of the structural models considered in our work, presenting the research of their equilibrium configurations as an energy minimisation problem (Section 2). Hence we discuss the mathematics behind their finite element discretisation and some details about their practical implementation (Section 3). Finally we show some comprehensive verification examples demonstrating what is possible with our approach and the performance of the proposed methods (Section 4). Three detailed documented demos are included as supplementary material and several more are available in the online repository of the python library FEniCS-Shells, companion to this paper.

A permanent DOI [START_REF] Hale | FEniCS-Shells[END_REF] has been created with links to the latest code and documentation. Table 1 summarises the mathematical models and the numerical methods presented in each of the documented demos currently available online. We suggest the reader consult the demos to have an overview of the main capabilities and features of FEniCS-Shells, and in particular the Clamped Reissner-Mindlin plate under uniform load demo for linear plates, the Buckling of a heated von-Kármán plate for weakly nonlinear plate models, and the Clamped semi-cylindrical Naghdi shell under point load demo for fully nonlinear shells. FEniCS-Shells follows the Best Practices for Scientific Computing [START_REF] Wilson | Best Practices for Scientific Computing[END_REF] as closely as possible, including using version control, continuous integration and testing, repeatable computing environments [START_REF] Hale | Containers for Portable, Productive, and Performant Scientific Computing[END_REF] and providing thorough documentation.

Structural models

This section presents the variational formulation of the equilibrium problem for the basic plate and shell models considered in this work. Starting from the most general model, the nonlinear Naghdi shell model [START_REF] Naghdi | Foundations of elastic shell theory[END_REF], the linear Reissner-Mindlin plate model [START_REF] Reissner | On bending of elastic plates[END_REF][START_REF] Mindlin | Influence of rotatory inertia and shear on flexural motions of isotropic, elastic plates[END_REF] is obtained as a special case of the linear Naghdi shell model [START_REF] Naghdi | The theory of shells and plates[END_REF]. In passing, we present also the so-called Marguerre-von Kármán shallow shell model [START_REF] Marguerre | Zur theorie der gekrümmten platte grosser formänderung[END_REF], because (in its unshearable version) it is extensively used in the literature to study nonlinear phenomena in plates and shallow shells [START_REF] Mansfield | The bending and stretching of plates[END_REF][START_REF] Lewicka | Models for elastic shells with incompatible strains[END_REF][START_REF] Audoly | Elasticity and Geometry. From hair curls to the non-linear response of shells[END_REF]. We review below the basic notations used in this work. [START_REF] Wells | A c0 discontinuous galerkin formulation for kirchhoff plates[END_REF][START_REF] Nguyen | Geometrically nonlinear formulation for thin shells without rotation degrees of freedom[END_REF], not discussed in the present paper.

Notation. We equip the Euclidean space with a Cartesian frame {O; a 1 , . . . , a n } with orthogonal unit basis vectors a i ∈ V n , V n being the translation space. Let x = X -O ∈ V n denote the position vector of the point X whose rectangular Cartesian coordinates are (X 1 , . . . , X n ) ∈ R n . We denote by R m×n the space of m × n matrices and by S n = {S ∈ R n×n : S = S T } the space of symmetric n × n matrices, S n + its subset of positive definite matrices. With a notational abuse, for any tensor field

T ∈ V n ⊗ • • • ⊗ V n , let T ∈ R nו••×n
denote also the collection of its Cartesian components. The same applies to vectors and tensors on manifolds. We denote by a • the inner product between two tensors. In particular if T is a tensor of order 2n and v a tensor of order n, we denote by T v the application of the tensor T on v and by T v • v the quadratic form on the space of v having T as metric. For v ∈ R n and x ∈ R m , ∇v ∈ R n×m denotes the linear operator such that dv = ∇v dx + o(|dx| 2 ). We often use also index notation, where Latin indices (e.g. i, j) are assumed to take the values 1, 2, 3, whilst Greek indices (e.g. α, β) take the values 1, 2. For the partial derivatives of a function f we use the abridged notation f, α = ∂f /∂x α . We define the space L 2 (ω) as the usual Sobolev space of square-integrable functions on the domain ω with respect to the Lebesgue measure dx. Then H s (ω) is the usual Sobolev space of L 2 (ω) functions whose weak derivatives of order s ∈ N are also in L 2 (ω). We use the notation f ∈ H s (ω; R n ) to denote a vector-valued function whose components (f 1 , • • • , f n ) are each functions in H s (ω).

Non-linear Naghdi shell model

Kinematics. We describe the reference configuration of shell-like bodies as:

p 0 (x α , x 3 ) = φ 0 (x α ) + x 3 n 0 (x α ) , x α ∈ ω ⊂ R 2 , x 3 ∈ [-t/2, t/2],
where (x α , x 3 ) are curvilinear coordinates, φ 0 : ω → R 3 is a (piece of) regular surface, and

n 0 = (∂ 1 φ 0 × ∂ 2 φ 0 )/ |∂ 1 φ 0 × ∂ 2 φ 0 |
(× refers to the vector cross-product in only this context) is a continuous unit vector field normal to the middle surface φ 0 (ω), Fig. 1. The vectors g 0,α = ∂ α φ 0 define a covariant basis on the tangent plane at each point of the reference configuration. The slenderness condition requires the thickness t to be much smaller that the diameter of the middle surface ω. The geometry of the middle surface is completely described by the metric and curvature tensors, a 0 and b 0 respectively, whose covariant components are computed as: where ∇φ 0 = [g 0,1 , g 0,2 ]. Whenever the reference configuration is chosen to be flat, that is when φ 0 (ω) is a planar surface, n 0 ≡ a 3 (and x 3 ≡ X 3 ) and b 0 = 0. In such a case a body is said plate-like.

a 0 = ∇φ T 0 ∇φ 0 , b 0 = - 1 2 ∇φ T 0 ∇n 0 + ∇n T 0 ∇φ 0 , a 0 ∈ S 2 + , b 0 ∈ S 2 , n O 0 n d x 3 x 3 Reference configuration Current configuration p p 0 u Parametric domain x 1 x 2 (x 1 , x 2 ) o X 1 X 2 X 3 x 3 X 3 d X 1 X 2 x 3 β 1 d β 2
In the nonlinear Naghdi shell model the current configuration of the shell is written as a sum of the placement of the middle surface φ : ω → R 3 (or alternatively the displacement field u = φ -φ 0 ) and of the director field d : ω → R 3 (see Figure 1):

p(x α , x 3 ) = φ(x α ) + x 3 d(x α ).
In general, |d| = 1 and d

• n = 0, where n = (∂ 1 φ × ∂ 2 φ)/ |∂ 1 φ × ∂ 2 φ|
is the unit normal to the current surface and the vectors g α = ∂ α φ define a covariant basis on the tangent plane at each point of the deformed configuration. When the normal fibres are assumed to be unstretchable, |d| = 1 and the director can be parametrised by two angles β = (β 1 , β 2 ) : ω → R 2 , e.g., d(β) = (cos β 1 sin β 2 , -sin β 1 , cos β 1 cos β 2 ), giving rise the so-called five-parameter shell model [START_REF] Betsch | On the parametrization of finite rotations in computational mechanics: A classification of concepts with application to smooth shells[END_REF] (the three-components of the displacement field and the two angles (β 1 , β 2 )), Figure 1.

In a Lagrangian description, we use the following measures of deformation for the shell [START_REF] Naghdi | Foundations of elastic shell theory[END_REF][START_REF] Mardare | On the derivation of nonlinear shell models from three-dimensional elasticity[END_REF]:

e N (u) = 1 2 ∇φ(u) T ∇φ(u) -a 0 , k N (u, β) = - 1 2 ∇φ(u) T ∇d(β) + ∇d(β) T ∇φ(u) -b 0 , γ N (u, β) = ∇φ(u) T d(β), (1)
representing the stretching of the middle plane (a second order tensor), the bending curvature (a second order tensor), and the shear strain vector, respectively. The shear strain vanishes if and only if the director d coincides with the unit normal to the surface.

Strain energy and constitutive equations. We consider elastic shells characterised by an elastic strain energy density in the reference configuration ψ(e, k, γ), a scalar function of the deformations. The stress measures are the membrane stress resultant N = ∂ e ψ, a second order tensor, the couple stress resultant M = ∂ k ψ a second order tensor, and the vector T = ∂ γ ψ representing the transverse shear stress. We will focus in our examples on linearly elastic constitutive behaviour for which the strain energy density is a quadratic function of the strain measures. Thus, using standard notation:

ψ(e, k, γ) = 1 2 A(e -ē) • (e -ē) + B(e -ē) • (k -k) + 1 2 D(k -k) • (k -k) + 1 2 Sγ • γ,
where A, D are fourth-order and positive-definite tensors providing the membrane and bending stiffness, B is a fourth order tensor providing the membrane-to-bending constitutive coupling, S is a second order tensor providing the shear stiffness and ē, k are symmetric second order tensor fields representing inelastic stretching and bending strains. In what follows, unless otherwise stated, ē, k = 0. This leads to the following linear constitutive laws:

N = A(e -ē) + B(k -k), M = B(e -ē) + D(k -k), T = Sγ.
The stored energy density is the sum of membrane, ψm , bending, ψb , and shearing, ψs , contributions:

ψ(e, k, γ) = ψm + ψb + ψs = 1 2 N • (e -ē) + 1 2 M • (k -k) + 1 2 T • γ.
For a shell made of a single homogeneous layer of St. Venant-Kirchhoff isotropic material the contravariant components of the elastic tensors are:

A αβστ t = 12
D αβστ t 3 = 2λµ λ + 2µ a αβ 0 a στ 0 + µ a ασ 0 a βτ 0 + a ατ 0 a βσ 0 , B αβστ = 0, S αβ = µa αβ 0 ,
where a αβ 0 are the contravariant components of the metric tensor a 0 and λ, µ the Lamé constants. Denoting by = t/L the small thinness parameter, |A|, |S| ∝ , while |B| ∝ 2 and |D| ∝ 3 . Pure bending deformations (e = 0, γ = 0) are energetically cheaper and, whenever possible, they are preferred.

Variational formulation of the equilibrium condition. We focus here on the solution of quasi-static problems under the action of conservative loads, although this is not a limitation of the FEniCS-Shells library. The equilibrium configurations are found by solving for the stationary points of the potential energy:

Π N (u, β) = ω ψ(e N (u), k N (u, β), γ N (u, β)) j 0 dx -W ext , (2)
where W ext is the external work of the conservative forces and j 0 = det a 0 . Considering as Dirichlet boundary conditions imposed displacements ū on the part of the boundary γ u and on the rotation β on γ β , the space of admissible displacements and rotations are, respectively:

U ≡ {u ∈ H 1 (ω, R 3 ) | u = ū on γ u }, R ≡ {β ∈ H 1 (ω, R 2 ) | β = β on γ β }. (3)
The variational formulation of the equilibrium equation consist in finding (u, β) ∈ U × R such that

D ũ [Π N (u, β)] = 0, D β [Π N (u, β)] = 0, ∀(ũ, β) ∈ U 0 × R 0 , (4)
where U 0 and R 0 are the vector spaces associated to U and R, obtained by setting ū = 0 and β = 0 in (3). The stability of an equilibrium can by assessed by studying the sign of the second derivative of the energy, a quadratic form of the variations ũ, β.

Remark. Shell models with energy in the form (2) are shearable, since flexure is due to both bending and shearing. However, for very thin shells and plates the shearing strains turn out to be negligible. These structures are quite naturally modelled by means of unshearable or pure bending theories obtained by imposing the inner constraint γ = 0. Such theories can usually be recovered as limit models of shearable theories when the thinness parameter goes to zero, → 0. In what follows the shell models in the library are introduced in their 'shearable version'; their unshearable counterpart can be easily derived. We remark that the distinction between shearable and unshearable theories is of paramount importance, since the numerical solution procedures exhibit quite different challenges.

Remark. By assuming the inner constraint γ N = 0, from Eq. (1) 3 d = n ∈ ker ∇φ T (u). Moreover, from Eq. (1) 2 , the change of curvature can be written as k(u) = b(u) -b 0 , b(u) being the curvature tensor of the current surface. In such a way the nonlinear Koiter shell model [START_REF] Koiter | A consistent first approximation in the general theory of thin elastic shells. The theory of thin elastic shells[END_REF] is recovered as the unshearable counterpart of the Naghdi's one for thin shells. Since b αβ = n i , β (φ 0 + u), α = n i (φ 0 + u), αβ , (i = 1, 2, 3), the displacement u must be sought in a space U ⊆ H 2 (ω, R 3 ). This high regularity requirement is a common feature of the unshearable theories.

Marguerre-von Kármán shallow shell model

Some widely studied and used weakly non-linear versions of the Naghdi shell equations are the Föppl-von Kármán plate model [START_REF] Foppl | Vorlesungen uber technische Mechanik[END_REF][START_REF] Kármán | Festigkeitsprobleme im maschinenbau[END_REF] and the Marguerre shallow shell model [START_REF] Marguerre | Zur theorie der gekrümmten platte grosser formänderung[END_REF]. They can be seen as the simplest nonlinear models able to properly describe the moderately large deflections of thin shallow shells or plates, since they retain a minimal geometrical coupling between the membrane and flexural behaviour. Here, we introduce them in their shearable version [START_REF] Kere | Reissner-mindlin-von kármán type plate model for nonlinear analysis of laminated composite structures[END_REF].

Consider a shell with a shallow initial configuration φ 0 (x α ) = x α a α + ζ(x α )a 3 , with the shallowness assumption |∇ζ| ∼ η 1 and let φ(x α ) = x α a α + v(x α ) + w(x α )a 3 be the current configuration, where e 3 of the middle plane with respect to a flat configuration x = x α e α ∈ ω, and decompose the v ⊥ ω is the in-plane component and w on the a 1 -a 2 plane, with v • a 3 = 0, and the transverse component wa 3 . Similarly, we denote by θ the (small) rotations of the fibres orthogonal to the middle-plane, with θ • a 3 = 0.

The displacement of an arbitrary point of the three-dimensional body can be represented as:

U M (x α , x 3 ) = v(x α ) + w(x α )a 3 -x 3 θ(x α ). (5)
The Marguerre-von Kármán shallow shell model can be derived from the nonlinear Naghdi model by assuming

|v|/L ∼ η 2 , |w|/L ∼ η and |θ| ∼ η, so that d(x α ) -n 0 (x α ) = θ(x α ) + h.o.t., where θ(x α ) = β(x α ) × a 3 .
Then, the Marguerre-von Kármán strain measures are obtained as the leading order terms of the Naghdi strain measures, Eqs. 1,

e M (v, w) = ∇v + ∇v T 2 + ∇w ⊗ ∇w 2 - ∇ζ ⊗ ∇ζ 2 , k M (θ) = ∇θ + ∇θ T 2 -∇∇ζ, γ M (w, θ) = ∇w -θ. (6)
Then, solutions are found as the minimisers over a smooth set of vector fields

v ∈ V ⊆ H 1 (ω, R 2 ), w ∈ W ⊆ H 1 (ω), θ ∈ R ⊆ H 1 (ω, R 2
) satisfying appropriate boundary conditions, of the nonlinear Marguerre-von Kármán energy functional:

Π M (v, w, θ) = 1 2 ω N M (v, w) • e M (v, w) -ē dx + 1 2 ω M M (w, θ) • k M (w, θ) -k dx + 1 2 ω T M (w, θ) • γ M (w, θ) dx -W ext . (7)
Remark. By assuming the inner constraint γ M = 0, from Eq. (6) 3 , θ = ∇w. Moreover, from Eq. (6) 2 , k(w) = ∇∇w. In such a way the classical (unshearable) Marguerre-von Kármán shallow shell model is recovered. As for the Koiter model [START_REF] Koiter | A consistent first approximation in the general theory of thin elastic shells. The theory of thin elastic shells[END_REF], the transverse displacement w must be sought in a space W ⊆ H 2 (ω).

Linear Naghdi shell model

Consider a shell body with reference surface φ 0 (ω). Again, let U Ñ (x α , x 3 ) = u(x α )a 3 -x 3 θ(x α ) be the representation formula for the displacement vector field, with θ ∈ T P φ 0 (x α ), T P φ 0 (x α ) being the tangent plane to the reference surface at the point φ 0 (x α ); u represents an (infinitesimal) displacement from the reference surface, θ denotes the (infinitesimal) rotation of the fibers orthogonal to the middle surface. The linear Naghdi strain measures are then obtained by linearisation from Eqs. ( 1),

e Ñ (u) = 1 2 ∇φ T 0 ∇u + ∇u T ∇φ 0 , k Ñ (u, θ) = - 1 2 ∇φ 0 T ∇θ + ∇θ T ∇φ 0 - 1 2 ∇n 0 T ∇u + ∇u T ∇n 0 , γ Ñ (u, θ) = ∇φ 0 T θ + ∇u T n 0 . (8)
Then, solutions are found as the minimisers over a smooth set of vector fields

u ∈ U ⊆ H 1 (ω, R 3 ), θ ∈ R ⊆ H 1 (ω, R 2
) satisfying appropriate boundary conditions, of the linear Naghdi energy functional:

Π Ñ (u, θ) = 1 2 ω N Ñ (u) • e Ñ (u) dx + 1 2 ω M Ñ (u, θ) • k Ñ (u, θ) dx + 1 2 ω T Ñ (u, θ) • γ Ñ (u, θ) dx -W ext . (9)
Remark. For a plate-like body, Eqs. (8) reduce to:

e RM (v) = 1 2 ∇v + ∇v T , k RM (v, θ) = 1 2 ∇θ + ∇θ T , γ RM (w, θ) = ∇w -θ, (10)
where u(x α ) = v(x α )+w(x α ) a 3 . In such a case the membrane and flexural problem are geometrically decoupled and a linear coupling turns out to be possible only at the constitutive level. Eqs. (10) 2,3 are the Reissner-Mindlin plate strain measures, depending only by the transverse displacement w and the rotation θ. Again, solutions are found as the minimisers over a smooth set of vector fields

w ∈ W ⊆ H 1 (ω), θ ∈ R ⊆ H 1 (ω, R 2 )
, satisfying appropriate boundary conditions, of the Reissner-Mindlin energy functional:

Π RM (w, θ) = 1 2 ω M RM (w, θ) • k RM (w, θ) dx + 1 2 ω T RM (w, θ) • γ RM (w, θ) dx -W ext . ( 11
)
Remark. By assuming the inner constraints γ Ñ = 0 and γ RM = 0, the linear Koiter shell model [START_REF] Warner | On the foundations of linear theory of thin elastic shells. 1[END_REF] and the Kirchhoff-Love [START_REF] Kirchhoff | Uber das gleichgewicht und die bewegung einer elastischen scheibe[END_REF][START_REF] Love | A treatise on the mathematical theory of elasticity[END_REF] plate model are recovered, respectively. Solutions of such unshearable models must be sought in spaces U ⊆ H 2 (ω, R 3 ), W ⊆ H 2 (ω), respectively.

Finite element formulations and implementation details

Notation. We use nomenclature from [START_REF] Logg | Automated Solution of Differential Equations by the Finite Element Method[END_REF] related to finite element methods. We let T h (the mesh) be a partition of the domain ω into a finite set of triangular cells T h = {T } with disjoint interiors:

∪ T ∈T h T = ω.
Each triangular cell T has circumradius h T , and h := max T ∈T h h T . We then denote the set of edges of T h as E = {E}. E int h ⊂ E h is the subset of edges on the interior of the mesh (i.e. not on the boundary). Normals to the edges of each cell are denoted n. Functions u evaluated on opposite sides of an edge E are indicated by the subscripts + and -, giving the jump operator

u = u + • n + + u -• n -. We denote CG p ⊂ H 1 (ω) the standard continuous Lagrangian finite element space of polynomial order p constructed on the mesh T h . The notation f ∈ [CG p ] n ⊂ H s (ω; R n ) denotes a vector-valued function whose components (f 1 , • • • , f n ) are each a member of CG p .

Background

As pointed out in the previous section, shearable models give rise to weak formulations with solutions in H 1 while unshearable models give rise to weak formulations with solutions in H 2 . It is well-known that producing conforming finite element method discretisations of H 2 is not straightforward, although several elegant solutions are possible [START_REF] Cirak | Subdivision surfaces: a new paradigm for thin-shell finite-element analysis[END_REF][START_REF] Nguyen | Geometrically nonlinear formulation for thin shells without rotation degrees of freedom[END_REF][START_REF] Benson | Isogeometric shell analysis: The Reissner-Mindlin shell[END_REF][START_REF] Hale | A locking-free meshfree method for the simulation of shear-deformable plates based on a mixed variational formulation[END_REF]. To bypass this issue, one can choose to use shearable models even for very thin structures, even though in such cases shear strains are negligible. Shearable models only requires H 1 -conforming finite element methods which are significantly easier to implement than either discontinuous Galerkin or H 2 -conforming finite element methods.

Unfortunately it is also the case that applying standard H 1 conforming finite element methods leads to shear-locking in the thin limit, since the finite element spaces are unable to properly represent the zero shear-strain γ = 0 constraint asymptotically enforced within the mathematical model as the thinness parameter goes to zero. A similar problem arises in shell problems whenever pure bending displacements are non-inhibited: H 1 conforming finite element methods for such bending-dominated shells suffer of membrane-locking [START_REF] Chapelle | The Finite Element Analysis of Shells -Fundamentals[END_REF][START_REF] Sanchez-Hubert | Coques élastiques minces: propriétés asymptotiques[END_REF][START_REF] Sanchez-Palencia | Singular Problems in Shell Theory: Computing and Asymptotics[END_REF]. In other words, the discrete spaces are incompatible with the structural models.

A wide variety of locking-free numerical methods exist, including the methods of assumed natural strains [START_REF] Bathe | The MITC7 and MITC9 Plate bending elements[END_REF][START_REF] Chapelle | MITC elements for a classical shell model[END_REF][START_REF] Jeon | The MITC3+ shell element in geometric nonlinear analysis[END_REF], enhanced assumed strains [START_REF] Simo | A class of mixed assumed strain methods and the method of incompatible modes[END_REF][START_REF] Bischoff | Shear deformable shell elements for large strains and rotations[END_REF], discrete shear gap [START_REF] Bletzinger | A unified approach for shear-locking-free triangular and rectangular shell finite elements[END_REF], partial selective reduced integration [START_REF] Arnold | Locking-free finite element methods for shells[END_REF] and mixed variational [START_REF] Hale | A locking-free meshfree method for the simulation of shear-deformable plates based on a mixed variational formulation[END_REF] approaches have been developed to overcome this problem. Implementing these locking-free approaches is more complex than implementing the standard H 1 -conforming finite element method as they require mixed formulations and special finite element function spaces. Making these methods easier to code, use and adapt is one of the main goals of the present work and the companion open-source library FEniCS-Shells.

For alleviating shear-locking in the Reissner-Mindlin and Marguerre-von Kármán models we have chosen to implement a Mixed Interpolation Tensorial Components (MITC) reduction operator. The MITC approach was introduced in the paper of Bathe and Dvorkin [START_REF] Bathe | A four-node plate bending element based on Mindlin/Reissner plate theory and a mixed interpolation[END_REF] for curing shear-locking in the finite element discretisation of the Reissner-Mindlin model. Other notable contributions along the same lines include the Durán-Liberman [START_REF] Durán | On mixed finite element methods for the Reissner-Mindlin plate model[END_REF] and the MITC7 and MITC9 [START_REF] Bathe | The MITC7 and MITC9 Plate bending elements[END_REF] elements. The success of the various MITC elements for the Reissner-Mindlin problem rests on solid results in numerical analysis, e.g. [START_REF] Bathe | The MITC7 and MITC9 Plate bending elements[END_REF][START_REF] Ainsworth | The hp-MITC finite element method for the Reissner-Mindlin plate problem[END_REF]. It has also been shown numerically that MITC approach can be successful in solving membrane and shear locking issues in the linear and non-linear Naghdi shell models [START_REF] Jeon | The MITC3+ shell element in geometric nonlinear analysis[END_REF], although without rigorous mathematical justification. A fully general shell element with guaranteed mathematical and numerical behaviour in membrane and bending-dominated regimes is still out of reach and mathematically robust finite element design for shells is an active research topic. We have chosen to implement a high-order partial selective reduced integration (PSRI) method which was analysed for the linear Naghdi model in [START_REF] Arnold | Locking-free finite element methods for shells[END_REF].

The shear and membrane-locking issues

For notational convenience we drop the explicit RM suffix relating to the Reissner-Mindlin model in what follows. Informally, when considering the Reissner-Mindlin plate model, shear-locking is the result of the Kirchhoff constraint (vanishing shear strain γ) [START_REF] Hale | Meshless methods for shear-deformable beams and plates based on mixed weak forms[END_REF]:

γ = ∇w -θ = 0, ( 12
)
being enforced on the discrete finite element spaces w h ∈ W h ⊂ W and θ h ∈ R h ⊂ R for as the small thinness parameter → 0. The approximation power in the discrete finite element basis is inadequate to simultaneously represent the Kirchhoff (bending) dominated mode and also provide a uniformly convergent approximation to the problem with respect to . More formally, denoting the continuous space of pure bending displacements ( = 0) as:

V 0 = {(w, θ) ∈ W × R | ∇w -θ = 0} , (13)
and the discrete solution

(w h , θ h ) ∈ V h := W h × R h , we seek solutions in a subspace V h0 ⊂ (V h ∩ V 0 )
containing discrete pure bending displacements:

V h0 = {(w h , θ h ) ∈ W h × R h | ∇w h -θ h = 0} . ( 14
)
It is not so hard to construct discretisations in which this requirement is quite catastrophic. Consider using standard H 1 (ω)-conforming linear Lagrangian elements for both fields, W h := CG 1 and R h := [CG 1 ] 2 , and suppose fully clamped Dirichlet conditions on all of the boundary, that is

W h ⊂ H 1 0 (ω) and R h ⊂ H 1 0 (ω, R 2 )
. With piecewise linear functions for w h , ∇w h then consists of piecewise constant functions. Furthermore, as θ h must be zero on the boundary, and we are searching for solutions in the space V h0 where ∇w h = θ h , the only possible solution that our discrete space can represent is w h = 0 and θ h = 0. The finite element model is overly stiff, hence the term 'locking'.

The membrane-locking issue arises in the linear and non-linear Naghdi shell models and the Marguerre-von Kármán shallow shell model. The locking issue can be motivated by a similar argument to the shear-locking issue and causes similarly poor results. However, the precise circumstances under which membrane-locking can appear is a more subtle matter, and can depend on a combination of geometry, boundary conditions and loading. We refer the reader to [START_REF] Chapelle | The Finite Element Analysis of Shells -Fundamentals[END_REF] for an excellent discussion of locking in shell finite element models. Here, we limit ourselves to remarks we consider strictly necessary to make our exposé understandable.

For simplicity in what follows, we refer to the linear Naghdi shell model. In this case, the continuous space of pure bending displacements (no membrane or shear strains) reads as:

V 0 = {(u, θ) ∈ U × R | e N (u) = 0 and γ N (u, θ) = 0}. ( 15
)
When the thinness parameter goes to zero two distinct asymptotic behaviours emerge, depending on whether V 0 is empty or not, and thus on the shape of the shell and the Dirichlet boundary conditions applied on it. If V 0 = ∅, pure bending is inhibited and the shell carries the load mainly by membrane actions; this is the so called membrane-dominated regime. Conversely, whenever V 0 = ∅, pure bending is non-inhibited, the membrane and shearing contribution to the elastic energy vanish and the shell carries the load by bending actions; this is the so called bending-dominated regime.

When considering finite element discretisations, this latter case is likely to suffer both membrane and shear locking, unless the discrete space V h0 is chosen so as to properly approximate the space V 0 of pure bending displacements. In any case, since a finite element formulation should perform well regardless of the asymptotic regime, benchmark tests have to be provided for both membraneand bending-dominated shells.

Mixed Interpolation of Tensorial Components (MITC) method

We use the MITC approach to deal with the shear-locking issue. The MITC approach is based on an underlying mixed variational method, where the Reissner-Mindlin shear strains γ are considered as an independent variable. One of the most attractive aspects of the MITC approach is that the final linear system of equations only contains the original primal unknown fields (w h , θ h ) despite the consideration of the shear strain field γ h in the formulation.

The first step in constructing an MITC-type element is the correct identification of the function space S(ω; R 2 ) in which the shear strain γ lives. For any scalar field w ∈ H 1 (ω) it holds that ∇w ∈ H(rot; ω) where H(rot; ω) is the Sobolev space of square-integrable vector-valued functions with square-integrable rot:

H(rot; ω) := γ ∈ L 2 (ω; R 2 ) | rot γ := ∂ 1 γ 2 -∂ 2 γ 1 ∈ L 2 (ω) . (16)
Furthermore, for θ ∈ H 1 (ω; R 2 ) we have rot θ ∈ L 2 (ω) and thus it follows that

γ = ∇w -θ ∈ S ≡ H(rot; ω). (17)
Note that the operator rot is synonymous with the operator curl in R 2 . Finite element approximation spaces that satisfy S h ⊂ S consist of the well-known 'curl/rotconforming' elements, namely the Nédélec [START_REF] Nedelec | Mixed finite elements in R3[END_REF] (or rotated Raviart-Thomas elements [START_REF] Raviart | A mixed finite element method for 2-nd order elliptic problems[END_REF] in R 2 ) of the first NED n 1 and second kind NED n 2 , and the Brezzi-Douglas-Fortin-Marini (BDFM n ) elements [START_REF] Brezzi | Two families of mixed finite elements for second order elliptic problems[END_REF], where n is the element (polynomial) order. All of these finite element spaces are available for use in FEniCS [START_REF] Alnaes | The FEniCS Project Version 1.5[END_REF][START_REF] Logg | DOLFIN: Automated Finite Element Computing[END_REF][START_REF] Alnaes | Unified Form Language: A Domain-specific Language for Weak Formulations of Partial Differential Equations[END_REF]. We refer the reader to [START_REF] Logg | Automated Solution of Differential Equations by the Finite Element Method[END_REF] for a full description of the mathematical construction of these elements, and [START_REF] Rognes | Efficient Assembly of $H(\mathrm{div})$ and $H(\mathrm{curl})$ Conforming Finite Elements[END_REF] as a reference for the issues surrounding their implementation in FEniCS.

For brevity, we shown only the development of the Durán-Liberman element [START_REF] Durán | On mixed finite element methods for the Reissner-Mindlin plate model[END_REF]. This element can be considered the lowest order MITC-type element with rigorous analytical convergence proofs. There is also an implementation of MITC7 [START_REF] Bathe | The MITC7 and MITC9 Plate bending elements[END_REF] in FEniCS-Shells but do not show the development in this paper.

The starting point is the discrete counterpart of eq. ( 11):

Π(w h , θ h ) = Π b (w h , θ h ) + 1 2 ω T (w h , θ h ) • γ(w h , θ h ) dx -W ext , Π b (w h , θ h ) := 1 2 ω M (w h , θ h ) • k(w h , θ h ) dx. ( 18
)
The finite element spaces used in the construction of the Durán-Liberman element are continuous piecewise-linear Lagrangian elements for the transverse displacements w h ∈ W h := CG 1 , continuous piecewise-quadratic vector-valued Lagrangian elements for the rotations θ h ∈ R h := [CG 2 ] 2 and, finally, the lowest-order Nédélec elements of the first kind for the shear strain γ h ∈ S := NED 11 . The product of these spaces then makes up the Durán-Liberman element, see Figure 2.

The element NED 1 is defined [START_REF] Logg | Automated Solution of Differential Equations by the Finite Element Method[END_REF] by the Ciarlet triple (T, V, L) (cell, polynomial space, degrees of freedom):

T = triangle, (19a)
V = [P 0 (T )] 2 + S 1 (T ), (19b)
L = E v • t p ds, for p ∈ P 0 (E), ∀E ⊂ T, ( 19c
)
where T is the reference triangle with edges E with measure ds, P 0 is the set of constant polynomial functions on the triangle T or edge E, and S 1 (T ) is the polynomial space:

S 1 (T ) = s ∈ [P 1 (T )] 2 | s(x) • x = 0 ∀x ∈ T . (20)
The 'trick' of the MITC approach is recognising that by modifying the shear strain field at the discrete level by applying a special operator R h that takes the primal to the conforming space H(rot; ω) for the shear strains

R h : CG 1 × [CG 2 ] 2 → NED 1 (21)
the element can be 'unlocked'. For the Durán-Liberman element, the operator R h can be defined by the following condition for each cell T in the mesh T h :

E [(R h [γ(w h , θ h )] -γ h ) • t] p ds = 0, ∀p ∈ P 0 (E), ∀E. (22)
Note that ( 22) is very similar to the degree of freedom definition L in (19c). In effect, R h 'ties' the degrees of freedom of the transverse displacement and rotation space CG 1 × [CG 2 ] 2 to an underlying NED 1 space. Then, applying the operator R h to [START_REF] Betsch | On the parametrization of finite rotations in computational mechanics: A classification of concepts with application to smooth shells[END_REF] gives:

Π(w h , θ h ) = Π b (w h , θ h ) + 1 2 ω T (R h γ(w h , θ h )) • R h γ(w h , θ h ) dx -W ext . (23)
Remark. For the Reissner-Mindlin strain measure R h [γ(w h , θ h )] = ∇w h -R h θ h because R h w h = w h by definition. This application of the reduction operator to the rotation field is often found in the literature. However, in our work we always apply the reduction operator to the complete shear-strain variable γ.

New implementation of the MITC method in FEniCS-Shells

The formulation and implementation of MITC in FEniCS-Shells is slightly different to the standard approach. In our view it is simpler to understand as it is purely variational and extends the MITC shear-locking cure quite naturally to the non-linear models. Another key advantage is that it can be concisely expressed in UFL [START_REF] Alnaes | Unified Form Language: A Domain-specific Language for Weak Formulations of Partial Differential Equations[END_REF].

Starting with [START_REF] Betsch | On the parametrization of finite rotations in computational mechanics: A classification of concepts with application to smooth shells[END_REF] and this time directly substituting in the conforming strain field γ h ∈ S h gives:

Π(w h , θ h , γ h ) = Π b (w h , θ h ) + ω T (γ h ) • γ h dx -W ext . (24)
We then 'tie' the tangential component of the conforming strain field γ h •t to its counterpart calculated using the primal unknowns γ • t ≡ (∇w -θ) • t using a Lagrange multiplier field p h ∈ NED 1 that lives on all edges E ∈ E h of the mesh T h :

Π(w h , θ h , γ h , p h ) = Π b (w h , θ h )+ ω T (γ h )•γ h dx+ E∈E h E (γ(w h , θ h )•t-γ h •t)•(p h •t) ds-W ext . (25)
It is possible to implement the Durán-Liberman element in FEniCS by directly expressing (25) using UFL. However, this comes at the expense of two extra fields γ h and p h over the standard MITC implementation.

Fortunately it is possible to eliminate both γ h and p h at the cell-local level. To this end, we begin by defining the residual F (q; q) = 0 as the usual Gâteaux derivative of (25) at a point q h := (w h , θ h , γ h , p h ) in the direction of test functions q := ( w, θ, γ, p):

F (q h ; q) := D q [Π(q h )] = 0, ∀q ∈ CG 1 × [CG 2 ] 2 × NED 1 × NED 1 , (26)
which gives:

F (q h ; q) = a b (w h , θ h ; w, θ) + a s (γ h ; γ) + a Π h (w h , θ h , p h ; γ) + a Π h (γ h ; w, θ, p) -W ext (q) = 0, (27a)
where:

a b (w h , θ h ; w, θ) = ω M (w h , θ h ) • k( w, θ) dx, (27b)
a s (γ h ; γ) = ω T (γ h ) • γ dx, ( 27c
)
a Π h (w h , θ h , γ h ; p) = E∈E h E (γ(w h , θ h ) • t -γ h • t) • (p • t) ds. ( 27d
)
Deriving the Jacobian involves taking the Gâteaux derivative again, this time at a point q h in the direction of trial functions q := (w, θ, γ, p): J(q h ; q; q) = D q [F (q h ; q)] .

Because the Reissner-Mindlin equations are linear, this operation is trivial, giving:

J(q h ; q; q) = a b (w, θ; w, θ) + a s (γ; γ) + a Π h (w, θ, p; γ) + a Π h (γ; w, θ, p). ( 29
)
Using a standard Newton method, given an initial guess q 0 we seek for k = 0, . . . , n solutions δq := (δw, δθ, δγ, δp) to the linear systems: J(q k ; q; q) δq = -F (q k ; q), ( 30)

q k+1 = q k + δq. (31)
The left and right-hand side of the discrete Newton linear system corresponding to ( 29) and ( 27) respectively, has the following block structure:

  A 0 C 0 B D C T D 0     δz δγ δp   =   b z b γ b p   , (32)
where δz := (δw, δθ) and A, B correspond to the matrices associated with bilinear forms a b and a s in (29), respectively. The matrix D is a diagonal-form edge mass matrix arising from the bilinear form:

-

E∈E h E (γ • t) • (p • t) ds, (33)
and C arises from the bilinear form:

E∈E h E (γ(w h , θ h ) • t) • (p • t) ds. ( 34
)
The notation b (•) corresponds to a splitting of the residual vector into its sub-blocks on the functions z, γ and p. Solving for δγ with row two of (32):

δp = D -1 (b γ -Bδγ), (35)
and for δγ with row three of ( 32)

δγ = D -1 (b p -C T δz), (36)
and then eliminating δγ from (35) using: [START_REF] Hale | FEniCS-Shells[END_REF] gives

δp = D -1 (b γ -BD -1 b p + BD -1 C T δz), (37)
before finally eliminating δp from the first row of (32) using [START_REF] Hale | Containers for Portable, Productive, and Performant Scientific Computing[END_REF] and rearranging to give: Note again that D is diagonal and can thus be cheaply inverted at the cell-local level allowing [START_REF] Hale | Meshless methods for shear-deformable beams and plates based on mixed weak forms[END_REF] to be formed locally as well. The above procedure is implemented at the element level using a custom assembly function written in C++. FEniCS automatically returns the dense cell-local version of the block structured system (32) corresponding to the Jacobian (29) and residual [START_REF] Chinosi | Numerical analysis of some mixed finite element methods for Reissner-Mindlin plates[END_REF] forms. We then use the Eigen dense linear algebra package [START_REF] Gaël | [END_REF] to split up (32) into its constituent blocks before calculating A s and b s using [START_REF] Hale | Meshless methods for shear-deformable beams and plates based on mixed weak forms[END_REF]. The cell-local contributions A s and b s are then assembled into a sparse matrix and vector associated with global displacement unknowns δz only. Solving for the global unknown δz, we can then reconstruct the eliminated unknowns δγ and δp by post-processing at the cell-local level using ( 36) and [START_REF] Gaël | [END_REF], respectively. A Python code showing how to define the Reissner-Mindlin problem in FEniCS-Shells is given in Figure 3. The new ProjectedFunctionSpace class (line 9 in Figure 3) manages the relationship between full function space space with all unknowns q h := (w h , θ h , γ h , p h ) the reduced function space unknowns with displacement unknowns z = (w h , θ h ). The special custom assembly procedure in FEniCS-Shells is invoked via a call to assemble passing the ProjectedFunctionSpace object as the first argument, resulting in the elimination of δγ and δp at the local level (line 34 in Figure 3).

A s δz = b s , (38a)
A s = (A + CD -1 BD -1 C T ), (38b)
b s = b z + CD -1 BD -1 b p -CD -1 b γ . (38c)
The above procedure applies directly to the case of models with non-linear strain measures, e.g. the non-linear Naghdi shell model (see Section 2.1). The tying term in the Lagrangian is simply re-written using the appropriate strain measure γ N :

(γ N (u h , β h ) • t -γ N h • t) • (p h • t) ds, (39)
before being consistently linearised at the symbolic level using the tools in UFL. The resulting discrete Newton system has precisely the same block form as [START_REF] Durán | On mixed finite element methods for the Reissner-Mindlin plate model[END_REF], and so the same linear algebra assembly procedure can be used to eliminate the auxiliary unknowns γ N h and p h . We refer the reader to the documented demo Non-linear Naghdi roll-up cantilever [START_REF] Hale | FEniCS-Shells[END_REF] for an example applying MITC to a non-linear shell problem. Note, however, that we explicitly cure here only shear locking.

Our high-level approach makes it relatively straightforward to directly apply MITC to user-specified thin structural models.

Partial Selective Reduced Integration (PSRI)

We use the PSRI approach to simultaneously deal with the shear-and membrane-locking issues in the shell-like models. Our method is inspired by, but not identical to the one shown in [START_REF] Arnold | Locking-free finite element methods for shells[END_REF], which contains convergence proofs for the linear Naghdi model. The PSRI approach can be viewed as a reduced integration method [START_REF] Arnold | Locking-free finite element methods for shells[END_REF], a stabilised mixed method [START_REF] Chinosi | Numerical analysis of some mixed finite element methods for Reissner-Mindlin plates[END_REF], or an augmented Lagrangian method [START_REF] Boffi | Analysis of new augmented Lagrangian formulations for mixed finite element schemes[END_REF].

Following [START_REF] Arnold | Locking-free finite element methods for shells[END_REF] we begin with the discrete linear Naghdi energy functional:

Π Ñ (u, θ) = 1 2 ω N Ñ (u h ) • e Ñ (u h ) dx + 1 2 ω M Ñ (u h , θ h ) • k Ñ (u h , θ h ) dx + 1 2 ω T Ñ (u h , θ h ) • γ Ñ (u h , θ h ) dx -W ext , (40)
or in shortened form:

Π Ñ (u, θ) = Π m (u h ) + Π b (u h , θ h ) + Π s (u h , θ h ) -W ext , (41)
where Π m (u h ), Π b (u h , θ h ) and Π s (u h , θ h ) correspond to the membrane, bending and shear energies of the linear Naghdi model, respectively. The finite element spaces used in the construction of the PSRI element are second-order continuous piecewise-quadratic cubic-bubble-enriched Lagrangian elements for the displacements u h ∈ U h := [CG 2 + B 3 ] 3 and second-order continuous piecewise-quadratic Lagrangian elements for the rotations θ ∈ R h := [CG 2 ] 2 . The finite element space of cubic bubbles B 3 is defined by the Ciarlet triple (T, V, L):

T = triangle, V = {v ∈ P 3 (T ) | v| E = 0} , L = v (1/3, 1/3) . ( 42
)
We then introduce a parameter α ∈ R that splits the membrane and shear energy in the energy functional into a weighted sum of two parts:

Π Ñ (u, θ) = Π b (u h , θ h )+αΠ m (u h )+(1-α)Π m (u h )+αΠ s (u h , θ h )+(1-α)Π s (u h , θ h )-W ext , (43)
to which we apply reduced integration to the parts weighted by the factor 1 -α. PSRI is straightforward to implement in FEniCS using UFL [START_REF] Alnaes | The FEniCS Project Version 1.5[END_REF]. We use a fourth-order accurate quadrature scheme on the portions of the energy in eq. 43 weighted by α and a second-order accurate quadrature scheme on the portion of the energy weighted by 1 -α. As discussed in [START_REF] Arnold | Locking-free finite element methods for shells[END_REF] a more refined interpretation of PSRI is a stabilised mixed formulation where the extra stress unknowns are projected out at assembly time. The lowest-order design in [START_REF] Arnold | Locking-free finite element methods for shells[END_REF] suggests using a zero-order accurate (one-point) quadrature scheme, however in practice we have found this to produce spurious (zero-energy) modes.

The paper [START_REF] Arnold | Locking-free finite element methods for shells[END_REF] also suggests setting α = 1. This choice is sufficient to obtain analytical convergence results. However, for the Reissner-Mindlin problem, this leads to suboptimal convergence rates [START_REF] Lovadina | A New Class of Mixed Finite Element Methods for Reissner-Mindlin Plates[END_REF], and a better choice in practice is setting α T = 2 /h 2 T for each T ∈ T h . In numerical experiments we have found that the performance of the method is robust with respect to the parameter α. The combination of a high-order polynomial space with a PSRI approach for alleviating membrane-and shear-locking seems both simple and effective.

Curved Shells

Finally, we describe how the curved shell geometry required for the linear and non-linear Naghdi models is implemented in FEniCS. We use an exact geometry approach, describing the mid-surface of the shell as a closed-form map φ 0 : ω → R 3 . For example, for a cylindrical geometry we have (x 0 , x 1 ) ∈ ω := [-π/2, π/2] × [0, L] and the map φ 0 = {ρ sin x 0 , x 1 , ρ cos x 0 }. From this, we can calculate approximations of the metric a 0 and curvature b 0 tensors and the normal n 0 of the middle surface from the map φ 0 using UFL. Further derived quantities required, e.g. the contravariant components of the metric tensor a αβ 0 and the elasticity tensor A σβστ , can be straightforwardly expressed using UFL. An abbreviated Python code giving an indication for how the curved shell model can be implemented is given in Figure 4, extracted from the documented demo Clamped semi-cylindrical Naghdi shell under point load provided in the supplementary material.

Examples

To verify the convergence of the adopted discretisation techniques and illustrate the capabilities of the FEniCS-Shells library, we report some benchmark examples. We select eight test cases, including two linear Reissner-Mindlin plates (Section 4.1), two for linear Naghdi shells (Section 4.2), one for Marguerre-von Kármán shallow shells (Section 4.3), and three for nonlinear Naghdi shells (Section 4.4). The linear plate examples are discussed in more detail for testing the performance of our implementation of the MITC formulation; we also showcase here an adaptive remeshing example included in FEniCS-Shells. For linear shells we consider both membrane-and bending-dominated regimes while the non-linear shell examples serve to check them in large displacements and rotation regimes and are based on the FEniCS-Shells PSRI for curing the membrane and bending locking. The example on the thermal buckling of a Marguerre-von Kármán shallow shell of Section 4.3 is proposed here for the first time as a benchmark test of finite element shell models.

Reissner-Mindlin plate

In this section we demonstrate the correct convergence rate of our implementation of the Durán-Liberman element for the Reissner-Mindlin plate problem towards a manufactured analytical solution. Then, we show an implementation of a residual-based a posteriori error estimator which we use to drive mesh adaptivity, resolving the solution in a re-entrant plate problem with strong boundary layers. ','x [START_REF]Mat-fem: Learning the finite element method with matlab and gid[END_REF]','r*cos(x[0])'), r=rho, degree = 4) V_phi = FunctionSpace(mesh, VectorElement("Lagrange", triangle, degree = 2, dim = 3)) phi0 = project(initial_shape, V_phi) # Normal V_n = FunctionSpace(mesh, VectorElement("Lagrange", triangle, degree = 1, dim = 3)) nv = cross(phi0.dx(0), phi0.dx( 1)) n0 = project(nv/sqrt(inner(nv, nv)), V_n) # surface normal in the reference configuration # Director director = lambda beta: as_vector([sin(beta [START_REF]Mat-fem: Learning the finite element method with matlab and gid[END_REF])*cos(beta [0]), -sin(beta [0]), cos(beta [START_REF]Mat-fem: Learning the finite element method with matlab and gid[END_REF])*cos(beta [0])]) beta0_expression = Expression(["atan2(-n [START_REF]Mat-fem: Learning the finite element method with matlab and gid[END_REF], sqrt(pow(n[0],2) + pow(n[2],2)))", "atan2(n [0],n [2])"], n = n0, degree=4) V_beta = FunctionSpace(mesh, VectorElement("Lagrange", triangle, degree = 2, dim = 2)) beta0 = project(beta0_expression, V_beta) d0 = director(beta0) # director in the reference configuration # Finite element space P2, bubble = FiniteElement("Lagrange", triangle, degree = 2), FiniteElement("B", triangle, degree = 3) enriched = P2 + bubble U, Beta = VectorElement(enriched, dim=3), VectorElement(P2, dim=2) # element for displacements and rotations Q = FunctionSpace(mesh, MixedElement([U, Beta])) # Mixed element # Definitions of functions, test functions, and trial functions q_ = Function(Q); q, q_t = TrialFunction(Q), TestFunction(Q); u_, beta_ = split(q_) # Kinematics F = grad(u_) + grad(phi0) # gradient of the transformation d = director(beta_+ beta0) # direction # Initial metric and curvature a0 = grad(phi0).T*grad(phi0); a0_c = inv(a0); j0 = det(a0) b0 = -0.

Clamped square plate

We solve the problem of a unit-square Reissner-Mindlin plate Ω = (0, 1) 2 clamped with uniform thickness t = = 10 -4 , and with isotropic linear elastic material behaviour with Youngs modulus E = 10920 and Poisson ratio ν = 0.3. The plate is clamped on the boundary (homogeneous Dirichlet boundary conditions on the transverse displacements and rotations) and loaded by a uniform pressure. The closed-form polynomial expressions for the loading and the resulting solution fields are given in Appendix A.

Figure 5 shows the solution to the above problem obtained with FEniCS-Shells using a 32 × 32 cross-pattern mesh (left) and the error with respect to the analytical solution on a sequence of uniformly refined meshes in both the H 1 and L 2 norms (right). The convergence result given in the original paper of Durán and Liberman [START_REF] Durán | On mixed finite element methods for the Reissner-Mindlin plate model[END_REF] is:

||θ -θ h || H 1 (ω,R 2 ) + ||w -w h || H 1 (ω) ≤ Ch||f || L 2 (ω) , (44)
where C is some unknown positive constant independent of the cell size h, corresponding to a linear convergence in the H 1 -norm for both the transverse displacements and rotation variables. Although no proof is given in [START_REF] Durán | On mixed finite element methods for the Reissner-Mindlin plate model[END_REF] we achieve quadratic convergence for both variables in the L 2 norm, the canonical 'one order higher convergence in L 2 ' behaviour. Figure 6 shows the convergence in the H 1 -norm of the transverse displacement w and the rotations θ for varying thickness t. We see a very slight deterioration of the error for a thick plate, but the rate of convergence is nearly identical for thick through to very thin plates.

In conclusion, this benchmark indicates that the implementation of the Durán-Liberman element in FEniCS-Shells is locking-free and optimally convergent. We remark that a finite element method designed without consideration for locking would not converge at all for this problem. A commented code for this test case is given in the documented demo Clamped Reissner-Mindlin plate under uniform load (see supplementary material).

Partly-clamped L-shaped plate

This example shows how sophisticated and efficient adaptive remeshing strategies can be easily implemented using FEniCS-Shells. We consider the a posteriori error estimator proposed in [START_REF] Beirão Da Veiga | A-priori and a-posteriori error analysis for a family of Reissner-Mindlin plate elements[END_REF], Left: Error for the transverse displacement w (H 1 -norm) as a function of the mesh size h. Right: Error for the rotation field θ (H 1 -norm). In the latter case all lines lie on top of each other (convergence behaviour independent of t). We achieve optimal convergence in all cases.

10 2 10 1 h 10 1 H 1 ( ) error 1 1 t = 1E+00 t = 1E-02 t = 1E-04 10 2 10 1 h 10 1 H 1 ( , 2 ) error 1 1 t = 1E+00 t = 1E-02 t = 1E-04
which is similar to the one developed in [START_REF] Carstensen | A Posteriori Error Analysis for Conforming Mitc Elements for Reissner-Mindlin Plates[END_REF], but it includes extra terms that control the error at boundaries induced by the strong boundary layer present [START_REF] Arnold | Asymptotic analysis of the boundary layer for the Reissner-Mindlin plate model[END_REF] in the Reissner-Mindlin problem as → 0. The full specification of the error indicator is given in Appendix B. It is possible to implement these complex expressions in UFL with just a few lines of code and a demo in FEniCS-Shells demonstrates the procedure.

We consider the problem of a non-convex L-shaped plate with a re-entrant corner clamped along the two edges that form the re-entrant corner, and the other four sides left free, Figure 7. The plate is loaded with a uniform transverse load f scaled with the cube of the plate thickness t 3 . The material parameters are the same as before. The difficulties of resolving the low regularity solution in the straight corners and in the re-entrant corner, in addition to the strong boundary layers on the free edges, means that a uniform refinement strategy will be very inefficient. Using an a posteriori error estimator to guide the refinement process can lead to orders of magnitude efficiency improvements.

Because no exact analytical solution is known for the studied problem we cannot calculate the exact error. Therefore we resort to comparing the behaviour of the error estimator under both uniform and adaptive refinement strategies.

Figure 8a shows the convergence of the error estimator for a moderately thick plate t = 10 -1 using adaptive (blue solid) and uniform (green dashed) refinement. We can see that using uniform refinement we achieve an algebraic convergence rate of only O(d -1/4 ) where d is the number of degrees of freedom for the discrete system. This corresponds to a classical convergence rate of O(h 1/2 ). This sub-optimal behaviour is expected given the low regularity of the solution particularly in the re-entrant corner. The convergence of the adaptive scheme initially matches that for the uniform refinement before increasing to O(d -1/2 ) due to the strong local refinements in the re-entrant corner region, and also in the other corners during the latter refinement steps. In Figure 7a we show the meshes at the 2nd and 6th adaptive refinements, with 2,499 and 10,494 degrees of freedom, respectively. The strong refinement near the re-entrant corner is clearly evident.

Critically, we can achieve the same level of convergence with significantly fewer degrees of freedom using the adaptive scheme. Taking the most refined uniformly meshed problem, with 2,365,443 degrees of freedom, we can reach the roughly the same level of error using 123,495 degrees of freedom with the adaptive strategy. Figure 8: Plot of error estimator η against the number of cells in the L-shaped plate problem with using adaptive and uniform refinement. The rate of convergence is significantly improved using the adaptive refinement strategy compared with a uniform refinement strategy.

The trends for a very thin plate with t = 10 -4 shown in Figure 8b and Figure 7b are broadly similar to those in Figure 7a. However, in addition, because of the presence of the strong boundary layer on the free edges as t → 0, we can see that the error indicator preferentially refines these regions as well as the re-entrant corner.

Linear Naghdi shell

In this section we demonstrate the performance of the finite element discretisation of the Naghdi linear shell model implemented in the FEniCS-Shells library. The numerical results presented here are based the PSRI formulation. To verify the discretisation in both membrane-and bending dominated regimes two well known benchmark examples have been chosen, namely, a partly clamped hyperbolic paraboloid (bending-dominated test) and an axisymmetric hyperboloid with clamped ends (membrane-dominated test). In both cases, the computational domain is a rectangular domain ω which is mapped to the surface in R 3 through the action of the map φ 0 .

Partly clamped hyperbolic paraboloid

We consider a thin shell whose middle surface is an hyperbolic paraboloid; the shell is clamped on one of its edges and subject to a uniform load in the vertical direction, see Figure 9 for details. This is an excellent bending dominated benchmark problem to test a FE formulation with respect to membrane-locking issues [START_REF] Bathe | An evaluation of the mitc shell elements[END_REF]. The initial shape is the analytical surface

φ 0 = {x 0 , x 1 , x 2 0 -x 2 1 } with (x 0 , x 1 ) ∈ [-1/2, 1/2] × [-1/2, 1/2]
. For the simulation we set Poisson ratio to be ν = 0.3. The convergence analysis is performed for thicknesses t = 10 -n for n = 1, 2, 3, 4 and obtained by setting the uniform pressure load to p to be proportional to t. Convergence rates in strain energy are reported for different values of the dimensionless thinness parameter t = . The deformed shape of the shell shows two zones corresponding to radically different shell behaviours: a zone where geometry and constraints prevent pure bending displacements (inhibited zone) and a zone where they are allowed (non-inhibited zone). For very thin shells this implies a sharp change of the shell stiffness across the boundary separating the two zones, born out by the bending energy localisation [START_REF] Lee | On the asymptotic behavior of shell structures and the evaluation in finite element solutions[END_REF].

Axisymmetric hyperboloid with clamped ends

We consider a thin shell whose middle surface is an hyperboloid of revolution; the shell is clamped on both its edges and subjected to a surface axisymmetric normal pressure, so that symmetry arguments enable us to solve the problem over a smaller computational domain (1/8th of the whole domain), see Figure 10 for details. This is a well-posed membrane dominated benchmark problem [START_REF] Hiller | Measuring convergence of mixed finite element discretizations: an application to shell structures[END_REF]. The initial shape is the analytical surface φ 0 = {cos x 0 cosh x 1 , sin x 0 cosh x 1 , sinh x 1 } with (x 0 , x 1 ) ∈ [0, π/2] × [-1, 0]. For the simulation we set Poisson ratio ν = 1/3 and the pressure load f = f 0 cos(2x 0 ). Convergence rates in strain energy are reported for different values of the thinness parameter t. The deformed shape of the shell shows a boundary layer of width ≈ 6 √ t near the clamped edge can be recognised, where the mesh must adequately refined for proper convergence. These results are in agreement to what is reported in the literature [START_REF] Hiller | Measuring convergence of mixed finite element discretizations: an application to shell structures[END_REF].

Marguerre-von Kármán shallow shell

We discuss the performance finite element discretisations of the nonlinear Marguerre-von Kármán shallow shell model (see Section 2.2) implemented in the FEniCS-Shells library. We compare here the results obtained with the MITC and PSRI formulations.

Lenticular plate with inelastic curvature

We consider an initially flat isotropic linear elastic disk of radius a with lenticular thickness subject to a through-the-thickness temperature gradient. Thus, the initial shape of the plate is the disk x 2 0 + x 2 1 ≤ a 2 and the thickness varies as (1 -x 2 0 /a 2 -x 2 1 /a 2 ) t, where t is the maximal thickness at the center. This special shape has been proposed by Mansfield[59,[START_REF] Mansfield | The bending and stretching of plates[END_REF] to avoid boundary layer effects and obtain a prestressed plate allowing for an exact analytical solution in the nonlinear regime. The temperature gradient is modelled as a uniform isotropic inelastic curvature k with kxx = kyy = k T and kxy = 0 (with x 0 = x, x 1 = y). Mansfield has shown that in the framework of the von Kármán model the exact solution of the problem is with uniform curvature in space. The curvature tensor at the equilibrium is with k xy = 0 and

         for k T ≤ k * T , k x = k y : k x + 4 (1 + ν) 2 k 3 x k * 2 T = k T for k T > k * T , k x , k y = k T ± k 2 T -k * 2 T , k * T = t r 2 2 1 + ν 2(7 + ν) (1 + ν) . ( 45
)
For k T k * T , the behaviour is approximately linear. The equilibrium curvature k is spherical and close to k. However, conforming the plate to a spherical surface implies membrane extension. For k T k * T , the extensional deformation becomes too expensive from the energetic point of view and the minimal energy configurations are approximately cylindrical and bending-dominated. The critical value k * T corresponds to a bifurcation from the spherical shape. For k k * T bending and extensional energies are of the same order. In the perfectly isotropic case, for k > k * T the plate is neutrally stable: all the shapes obtained by rotating the curvature axis are still solutions of the problem. This structure was a source of inspiration for several recent works on morphing plates and shells [START_REF] Seffen | Morphing bistable orthotropic elliptical shallow shells[END_REF][START_REF] Seffen | Growth and shape control of disks by bending and extension[END_REF][START_REF] Hamouche | Multi-parameter actuation of a neutrally stable shell: a flexible gear-less motor[END_REF].

We perform simulations increasing the inelastic curvature k T , for a Poisson ration ν = 0.3 and maximum thickness t = 0.01. The Young modulus has no influence here. To break the rotational symmetry, we introduce a slight imperfection in the inelastic curvature imposing kx = k T /ρ, ky = k T ρ with ρ = 0.998. Figure 11 shows the bifurcation diagrams obtained with the MITC and the PSRI discretisations for the shearable version of the model exposed in Section 2.2. We use uniform unstructured triangular meshes and parametrise the mesh size by the number of division n along the radius. To avoid rigid body motions, we applied suitable pointwise Dirichlet boundary conditions. Namely we set to zero all the components of the displacement at the center, the displacement in x direction at (0, a) and the displacement in y direction at (a, 0). The analytical solution [START_REF] Kiendl | Isogeometric shell analysis with Kirchhoff-Love elements[END_REF] is reported in the plot as solid blue line. The results obtained with both the discretisations converge to the analytical solution. However, PSRI converges more quickly than MITC, giving very good results even with n = 5. A possible explanation is that membrane locking is not cured in our MITC implementation. We conclude that the PSRI technique introduced in FEniCS-Shells is a very efficient discretisation to cure both shear and membrane locking in the weakly nonlinear regime. This will be further tested and assessed in the next Section. To introduce fully nonlinear models, we compare in Figure 12 the results obtained with the fully nonlinear Naghdi shell model and the Marguerre-von Kármán shallow shell model. In both cases we use the PSRI discretisation. As expected [START_REF] Lewicka | Models for elastic shells with incompatible strains[END_REF], the results confirm that Marguerre-von Kármán model can safely be used to approximate the Naghdi model for sufficiently small loadings, namely for curvatures of the order of a 2 /t. More surprisingly, for the present test the discrepancies with respect the fully nonlinear model remain tolerable also for k 50 a 2 /t. The documented demo Buckling of a heated von Kármán plate in the supplementary material includes the full commented code for the Marguerre-von Kármán shallow shell with PRSI discretisation. An extension of this problem to the case of an orthotropic material is included in FEniCS-Shells as a further documented demo in order to illustrate some basic tools for modelling composite laminates that we do not explicitly describe in this paper.

Nonlinear Naghdi shell

We demonstrate here the performance of the finite element discretisation of the nonlinear Naghdi shell model implemented in the FEniCS-Shells library. To test the predictive capability of the formulation in the large displacements and rotations regime we selected three well known benchmark problems, see [START_REF] Jeon | The MITC3+ shell element in geometric nonlinear analysis[END_REF] and references therein.

Cantilever plate under end moment

We consider a strip clamped on one of its edges and subjected to a moment M on the other edge, see Figure 13 for details. In its natural configuration the shell occupies a flat rectangular domain of width b, length L. In the numerical simulations we set L = 12b, thickness t = b/10, and Poisson ratio ν = 0. For load M = 2πEI/L, where I = t 3 /12 the plate completely rolls up on itself. The normalised vertical and horizontal displacements at the free end have the analytical expressions w/L = (1 -cos 2πµ)/(2πµ) and v/L = sin (2πµ)/(2πµ) -1, respectively, where µ = M/M max . A comparison between the finite element and analytical solutions in terms of load-deflection curves is reported in Figure 13. The numerical solution is computed by discretising the domain with a 48 × 4 element mesh and shows excellent agreement with the analytical one. The insets show the deformed configurations of the plate for two values of the end moment.

Slit annular plate under end shear force

We consider a slit annular plate of internal radius R 1 = 6 and external radius R 2 = 10 clamped on one of its edges and subjected to a line shear load p on the other edge, see Figure 14 reported in Figure 14, where we set the thickness t = 0.03, Poisson ratio ν = 0 and maximum load p = 0.8E, with E = 21 × 10 6 . The numerical solution is computed by discretising the domain with a 28-division unstructured mesh and the MITC technique. The reference solution [START_REF] Sze | Popular benchmark problems for geometric nonlinear analysis of shells[END_REF] is computed on a 10 × 80 Abaqus S4R-element mesh. We see excellent agreement between the two results. We consider a semi-cylindrical shell clamped on one of its curved edges and subjected to a point load P at the centre of the free curved edge. Along its straight longitudinal edges the transverse displacement and the normal rotation are restrained, see Figure 15 for details.

The computational domain is a rectangle (x 0 , x 1 ) ∈ [-π/2, π/2] × [0, 3.048] which is mapped to the natural configuration semi-cylindrical shell of radius ρ = 1.016 through the initial transformation φ 0 (x 0 , x 1 ) = {ρ sin x 0 , x 1 , ρ cos x 0 }. A comparison between the FE and reference solution in terms of load-deflection curves is reported in Figure 15 for the thickness t = 0.03, Young modulus E = 2.0685 × 10 7 , Poisson ratio ν = 0.3, and a maximum load P max = 2000. The numerical solution is computed by discretising the domain with a 21-division unstructured mesh and shows a very good agreement with the reference solution in [START_REF] Sze | Popular benchmark problems for geometric nonlinear analysis of shells[END_REF] computed on a 40 × 40 abaqus S4R-element mesh. For a complete commented code for this test case we refer to the documented demo Clamped semi-cylindrical Naghdi shell under point load, available in the supplementary material.

Conclusions

In this paper we shown how a high-level approach to specifying thin structural models, in conjunction with automatic code generation tools, can lead to simple, extensible and effective finite element plate and shell formulations. The power of this approach is demonstrated via the open-source FEniCS-Shells library. We have shown the main mathematical aspects of the structural models available and the numerical procedures used to discretise them, with a particular emphasis on shear-and membrane-locking issues. We selected two main discretisation strategies: a special implementation of the MITC technique that provides an effective cure for shear-locking in linear plate models, and an extension to nonlinear shells of the PRSI strategy proposed in [START_REF] Arnold | Locking-free finite element methods for shells[END_REF] for nonlinear shell problems shear-and membrane-locking. The examples section, accompanied by a set of documented demos, while by no means exhaustive, shows the good performance of the proposed methods and the type of analysis that can be performed.

We believe that FEniCS-Shells will be a helpful contribution for researchers and graduate students in the community simulating thin structures. FEniCS-Shells could provide users with efficient and simple tools to implement specific solvers for multilayered or other specific plate and shell models [START_REF] Karam | Homogenization of Heterogeneous Thin and Thick Plates[END_REF], shell models including multiphysical couplings or nonlinear material behaviour [START_REF] Roman | Fracture path in brittle thin sheets: a unifying review on tearing[END_REF][START_REF] Klein | Shaping of elastic sheets by prescription of non-Euclidean metrics[END_REF][START_REF] Hamouche | Multi-parameter actuation of a neutrally stable shell: a flexible gear-less motor[END_REF], as well as one atom thick shell structures such as nanotubes and graphene [START_REF] Arroyo | Finite crystal elasticity of carbon nanotubes based on the exponential Cauchy-Born rule[END_REF].

The present work constitutes an instrumental step in a larger project on the mechanics and physics of shell structures, including the analysis of the nonlinear material behaviour and fracture of isotropic and anisotropic shells through phase-field variational models [START_REF] Bourdin | The variational approach to fracture[END_REF][START_REF] Ibarra | The tearing path in a thin anisotropic sheet from two pulling points: Wulff's view[END_REF][START_REF] Marigo | An overview of the modelling of fracture by gradient damage models[END_REF][START_REF] Amiri | Phase-field modeling of fracture in linear thin shells[END_REF], applications to multistable structures with embedded active materials [START_REF] Hamouche | Multi-parameter actuation of a neutrally stable shell: a flexible gear-less motor[END_REF][START_REF] Seffen | Prestressed morphing bistable and neutrally stable shells[END_REF], and for the understanding of singularities and energy scalings in plates and shells [START_REF] Conti | Confining thin elastic sheets and folding paper[END_REF]. Current and future technical developments of FEniCS-Shells will be aimed at providing anisotropic adaptive remeshing tools [START_REF] Sanchez-Palencia | Singular Problems in Shell Theory: Computing and Asymptotics[END_REF][START_REF] Micheletti | Reliability and efficiency of an anisotropic Zienkiewicz-Zhu error estimator[END_REF] and coupling to advanced nonlinear solvers, such as asymptotic numerical continuation methods [START_REF] Vannucci | An asymptotic-numerical method to compute bifurcating branches[END_REF] and deflation techniques [START_REF] Farrell | Deflation techniques for finding distinct solutions of nonlinear partial differential equations[END_REF] to automatically detect multiple solutions of nonlinear systems. A further interesting extension would be to adapt to shell modelling the special approach proposed by [START_REF] Rognes | Automating the solution of PDEs on the sphere and other manifolds in FEniCS 1.2[END_REF] for solving PDEs on manifolds with FEniCS.

where the mesh-dependant norm |||(θ, w)||| h is defined by: |||(θ, w)|||

Figure 1 :

 1 Figure 1: Kinematics of the nonlinear Naghdi shell model. Left: Parametric domain ω and placement of a generic point of the shell in the reference and current configuration. Right: Parameterisation of the director d in terms of the β1, β2 angles.

Figure 2 :

 2 Figure2: Finite element spaces making up the Durán-Liberman element. Circles denote point evaluations for degrees of freedom. Arrows denote point evaluations of tangential component for degrees of freedom. The Durán-Liberman has a total of 21 degrees of freedom: 3 for the displacement w h , 12 for the rotations θ h , 3 for the shear strain γ h , and 3 for the Lagrange multiplier p h . Thanks to the static condensation at the cell level, the latter 6 degrees of freedom for γ h and p h can be eliminated during the assembling process, reducing the number of effective degree of freedom per element to 15.

Figure 3 :

 3 Figure 3: Complete Python code for solving clamped Reissner-Mindlin problem with uniform transverse loading.

  from dolfin import *; from mshr import *; import numpy as np rho = 1.016; L = 3.048 # radius and length mesh = generate_mesh(Rectangle(Point(-np.pi/2., 0.), Point(np.pi/2., L)), 21) # Material properties and thickness E, nu, t = 2.0685E7, 0.3, 0.03; mu = E/(2.0*(1.0 + nu)); lmbda = 2.0*mu*nu/(1.0 -2.0*nu) # Cylindrical map initial_shape = Expression(('r*sin(x[0])

  5*(grad(phi0).T*grad(d0) + grad(d0).T*grad(phi0)) # Deformations e = lambda F: 0.5*(F.T*F -a0) # membrane deformation k = lambda F, d: -0.5*(F.T*grad(d) + grad(d).T*F) -b0 # bending deformation gamma = lambda F, d: F.T*d-grad(phi0).T*d0 # shear deformation # Contravariant Hooke's tensor. i, j, l, m = Index(), Index(), Index(), Index() A = as_tensor((((2.0*lmbda*mu)/(lmbda + 2.0*mu))*a0_c[i,j]*a0_c[l,m] + mu*(a0_c[i,l]*a0_c[j,m] + a0_c[i,m]*a0_c[j,l])), [i,j,l,m]) # Stress N = as_tensor(t*A[i,j,l,m]*e(F)[l,m],[i, j]) # membrane stress M = as_tensor((t**3/12.0)*A[i,j,l,m]*k(F,d)[l,m],[i, j]) # bending moment T = as_tensor(t*mu*a0_c[i,j]*gamma(F,d)[j], [i]) # shear stress # Energy densities psi_m = 0.5*inner(N, e(F)) # membrane energy psi_b = 0.5*inner(M, k(F,d)) # bending energy psi_s = 0.5*inner(T, gamma(F,d)) # shear energy # Total energy functional including PSRI h = CellDiameter(mesh); alpha = project(t**2/h**2, FunctionSpace(mesh,'DG',0)) dx_h = dx(metadata={'quadrature_degree': 2}) # measure with reduced integration Pi_PSRI = psi_b*sqrt(j0)*dx + alpha*psi_m*sqrt(j0)*dx + alpha*psi_s*sqrt(j0)*dx + (1.0 -alpha)*psi_s*sqrt(j0)*dx_h + (1.0 -alpha)*psi_m*sqrt(j0)*dx_h # First and second directional derivatives to calculate residual and Hessian F = derivative(Pi_PSRI, q_, q_t) J = derivative(F, q_, q)

Figure 4 :

 4 Figure 4: Abbreviated Python code showing construction of energies for half-cylinder non-linear Naghdi model (see Section 4.4). For further details we refer the reader to the documented demo Clamped semi-cylindrical shell under point load provided as the supplementary material.

Figure 5 :

 5 Figure 5: Chinosi clamped plate problem solved using Durán-Liberman element: Transverse displacements (scalar field, density plot) and rotations (vector field, quiver plot, 50% of points discarded) (left). Convergence to Chinosi analytical solution on sequence of uniformly refined meshes in H 1 and L 2 norms for t = 10 -4 (right).

Figure 6 :

 6 Figure 6: Convergence of Reissner-Mindlin problem discretised with Durán-Liberman element to Chinosi analytical solution for varying thickness parameter t.Left: Error for the transverse displacement w (H 1 -norm) as a function of the mesh size h. Right: Error for the rotation field θ (H 1 -norm). In the latter case all lines lie on top of each other (convergence behaviour independent of t). We achieve optimal convergence in all cases.

Figure 7 :

 7 Figure 7: Meshes after adaptive remeshing at different refinement levels (3rd and 6th level) for thick and thin plates.The increased refinement on the free boundaries for thin plates capture the strong boundary layer effect as t → 0.

  Thin plate (t = 10 -4 )

Figure 9 :

 9 Figure 9: Partly clamped hyperbolic paraboloid with linear Naghdi shell model: initial shape (left), convergence plot (centre), and deformed shape (right). We obtain convergence of O(d -1 ) in the energy for this problem.

Figure 10 :

 10 Figure 10: Clamped hyperboloid: initial shape (left), convergence plot (centre), and deformed shape (right).We obtain convergence of O(d -1 ) in the energy for this problem. This convergence requires strong refinement in the region of the boundary layer, shown on the deformed shape (right).

Figure 11 :

 11 Figure 11: Lenticular orthotropic plate. Bottom row: Bifurcation diagram for the curvatures kx and ky at the equilibrium as a function of the temperature loading. The curvature values reported in the plots as dots are the spatial average extracted from the finite element results with the MITC (left) and PRSI (right) models with uniform unstructured meshes with n ∈ [5, 10] divisions along the radius. The top row: Sketch of the lenticular plate (left) and snapshots of the equilibrium shapes before (center) and after (right) the bifurcation. The results are for a 2 /t = 100 and ν = 0.3.

Figure 12 :

 12 Figure 12: Lenticular orthotropic plate. Comparison of the results obtained with the von Kármán and the Naghdi models using the PSRI discretisation and an unstructured uniform mesh with n = 5 division along the radius. The results are for a 2 /t = 100 and ν = 0.3.

Figure 13 :

 13 Figure 13: Cantilever plate under end moment. Left: initial geometry and loading. Right: Comparison between analytical and numerical solution (MITC) for the moment-displacement diagram, where the insets show the deformed configurations for M = 0.5Mmax and M = Mmax. We report here the results for displacement at the free end in the X1 and X3 directions, normalised with the initial length L.

Figure 14 :

 14 Figure 14: Slit annular plate subjected to a line vertical load. The load-deflection diagram compare the numerical solution obtained with FEniCS-Shells MITC discretisation with a reference solution (Sze [87], Abaqus S4R). Right: geometry and loading; Left: final deformed configuration. Excellent agreement between the result of Sze [87] and FEniCS-Shells is achieved.

4. 4 . 3 .Figure 15 :

 4315 Figure 15: Clamped semi-cylindrical shell under point load. Left: Initial configuration and boundary conditions. Right: Comparison of the results of the PRSI implementation with a reference solution (Sze [87], Abaqus S4R). The insets show the deformed configurations for P = Pmax and P = Pmax/2. Notice the stiffening of the shell beyond a load of around 750N.

Table 1 :

 1 Structural models and discretisation techniques presented in main documented demos provided in FEniCS-Shells[START_REF] Hale | FEniCS-Shells[END_REF]. All the demos adopt shareable shell theory, except the Clamped Kirchhoff-Love platewhich is based on the discontinuous Galerkin discretisation technique presented in

	Title

  2 h = ||θ|| 2 H 1 (Ω) + T + t 2 )||f + div γ h || 2 L 2 (T ) + h 2 T || div M h + γ h || 2 L 2 (T ) , (B.3b) η 2 int = h E (h 2 E + t 2 )|| γ h • n || L 2 (E) + h E || M h n || 2 L 2 (E) , (B.3c) η 2 free = h E ||M h n|| 2 L 2 (E) + h E (h 2 E + t 2 )||γ h • n|| 2 L 2 (E) , (B.3d) ρ 2 = || rot(I -R h )θ h || 2 L 2 (T ) + ||(I -R h )θ h || 2 L 2 (T ) , (B.3e)where E free h ⊂ E int h are the subset of edges on the free boundary (natural boundary conditions).Stéphane P. A. Bordas thanks partial funding for his time provided by the European Research Council Starting Independent Research Grant (ERC StG) Grant No. 279578 entitled "Towards real time multiscale simulation of cutting in non-linear materials with applications to surgical simulation and computer guided surgery".Matteo Brunetti and Corrado Maurini acknowledge the financial support of Project ANR-13-JS09-0009 (Agence Nationale de la Recherche, 2014).

			T ∈T h	1 t 2 + h 2 T	||∇w -θ|| 2 L 2 (T ) ,	(B.2)
	can be controlled by an error estimator of the following form:
	η 2 h =	(η 2 + ρ 2 ) +	η 2 int +	η 2 free ,	(B.3a)
	T ∈T h	E∈E int h		E∈E free h
	η 2 = h 2 T (h 2

The original element construction of Durán and Liberman uses a slightly smaller quadratic space constructed as a sum of linear polynomials and quadratic bubbles. As this space is not available in FEniCS we use the full quadratic space instead.

Acknowledgements

We would like to thank Douglas N. Arnold for sharing his implementation of the unreduced Durán-Liberman element in the Unified Form Language with us. Jack S. Hale is supported by the National Research Fund, Luxembourg, and cofunded under the Marie Curie Actions of the European Commission (FP7-COFUND) Grant No. 6693582.

Supporting information

Code (LGPLv3), documentation and demos are available from the links at: https://doi.org/ 10.6084/m9.figshare.4291160.

Three documented demos are formally enclosed as supplementary materials:

(S1) Clamped-Reissner-Mindlin-plate-under-uniform-load.pdf (S2) Buckling-heated-von-Karman-plate.pdf (S3) Clamped-semi-cylindrical-w-shell-under-point-load.pdf

Appendix A. Closed-form solution for a Reissner-Mindlin clamped square plate with polynomial loading

This analytical solution is taken from the paper [START_REF] Chinosi | Numerical analysis of some mixed finite element methods for Reissner-Mindlin plates[END_REF]. With transverse loading:

f (x, y)/D = 12y(y -1)(5x 2 -5x + 1) 2y 2 (y -1) 2 + x(x -1)(5y 2 -5y + 1)

, it is possible to obtain the following closed-form solutions for the rotations and transverse displacements:

Residual-based a posteriori error estimator

The key result found in the paper [START_REF] Beirão Da Veiga | A-priori and a-posteriori error analysis for a family of Reissner-Mindlin plate elements[END_REF] is that the error between the finite element solution and the exact solution defined in the following sense: