Background

The Reissner-Mindlin shell

Hermit solves the linear-static Reissner-Mindlin (first-order shear deformation) shell [Min51, Rei45]. The kinematics carry a mid-surface displacement \(\mathbf{u}\) and an independent rotation \(\boldsymbol{\theta}\) of the shell director, so transverse shear is retained (unlike Kirchhoff-Love).

Working in the local element frame \((\mathbf{e}_0, \mathbf{e}_1, \mathbf{e}_2)\) built from the cell normal and Jacobian, the membrane strain \(\boldsymbol{\varepsilon}\), bending curvature \(\boldsymbol{\kappa}\) and transverse shear \(\boldsymbol{\gamma}\) are linear in \((\mathbf{u}, \boldsymbol{\theta})\). The stored energy is

\[\begin{split} \Pi = \tfrac12 \int_\Omega \begin{Bmatrix}\boldsymbol{\varepsilon}\\ \boldsymbol{\kappa}\end{Bmatrix}^{\!\top} \begin{bmatrix}\mathbf{A} & \mathbf{B}\\ \mathbf{B} & \mathbf{D}\end{bmatrix} \begin{Bmatrix}\boldsymbol{\varepsilon}\\ \boldsymbol{\kappa}\end{Bmatrix} \, \mathrm{d}\Omega + \tfrac12 \int_\Omega \boldsymbol{\gamma}^{\!\top} \mathbf{A}_s \boldsymbol{\gamma} \, \mathrm{d}\Omega , \end{split}\]

with the ABD constitutive matrices: \(\mathbf{A}\) membrane, \(\mathbf{D}\) bending, \(\mathbf{B}\) membrane-bending coupling, \(\mathbf{A}_s\) transverse shear. A small drilling-stabilization term regularises the in-plane rotation.

Constitutive input

  • Isotropic (hm.isotropic): \(\mathbf{A}, \mathbf{B}, \mathbf{D}, \mathbf{A}_s\) are assembled in plain CSDL from \(E\), \(\nu\) and the thickness, with a \(0.833\) shear-correction factor. Each of \(E\), \(\nu\), thickness, density may be a single value or a Field on its own FE space; a scalar is broadcast to the mesh.

  • Composite (hm.laminate): classical lamination theory — the ABD of a ply stack, differentiable with respect to ply angles and thicknesses.

  • Pre-computed (hm.composite): ABD fields supplied directly.

The residual always takes ABD, so there is no isotropic/composite branching in the UFL form.

Fibre orientation

hm.laminate / hm.composite produce the ABD in the laminate’s own 0-deg frame. On a real shell that frame is not the element local frame (which follows the mesh parametrisation). Pass fiber_direction (a global vector, or a per-cell field) or fiber_angle (per-cell radians from the element e0); Hermit rotates the ABD into each element’s local frame before the solve. A fiber_angle csdl.Variable is a differentiable design field. Orientation may use any supported field space.

Discretization

A mixed space: displacement in CG2 (quadratic Lagrange), rotation in CG1, on quad or triangle meshes (ShellDomain(element="CG2CG1"), the default). Alternatives: "CG1CG1" (linear displacement), and "CG2CR1" (Crouzeix–Raviart rotation, a non-conforming element that resists shear locking on thin shells — triangle meshes only).

Constant-normal bending curvature

The bending curvature is a shallow-shell, flat-facet measure: within each cell the normal is held constant in the bending kinematics. With \(t_{ij}=\mathbf e_i\mathbin{\cdot}\partial_{\mathbf e_j}\boldsymbol\theta\),

\[\begin{split} \boldsymbol\kappa = \begin{bmatrix} -t_{10} & (t_{00}-t_{11})/2\\ (t_{00}-t_{11})/2 & t_{01} \end{bmatrix}, \end{split}\]

so curvature comes from \(\nabla\boldsymbol\theta\) alone. Under an infinitesimal rigid motion \(\boldsymbol\theta\) is constant and \(\boldsymbol\kappa\) vanishes structurally, in warped quadrilateral cells as well as planar triangles. Stress recovery uses the same kinematics, \(\boldsymbol\varepsilon(\xi)=\boldsymbol\varepsilon_\mathrm{mid} -\xi\boldsymbol\kappa\), so displacement, strain energy, von Mises stress and the stress aggregates agree about rigid-body motion.

Constant-normal flat-facet approximations belong to the established shallow-shell element lineage, of which Belytschko–Tsay is a representative [BLT84].

Limitations

The measure drops a term coupling initial curvature to membrane stretch, so it is exact only in the shallow-shell limit. The error grows with the thickness-to-radius ratio and with laminate unsymmetry, and is largest when the shell is also membrane-loaded. Bending-dominated response is insensitive to the dropped term: the hypar cantilever and the twisted beam are unaffected by it, and both meet their references in the verification suite.

In practice: doubly curved geometry does not require triangles to obtain rigid-body objectivity, and warped quadrilaterals are supported by the same theory. For a moderately thick shell with strong initial curvature and significant membrane loading — especially an unsymmetric laminate — run a thickness and mesh study, and do not read the result as a full covariant shell solution. The Reissner–Mindlin formulation of Schöllhammer and Fries is the relevant direction for that regime [SF19].

Loads

Every load builder returns a composable Loads; a + b merges them, and each term reaches the residual on its own FE space with no interpolation between terms.

Distributed over the surface

  • hm.pressure(domain, p) — a scalar load along the shell normal, positive along \(+\mathbf{n}\). The residual uses the live ufl.CellNormal, so it is shape-differentiable. It requires a consistently wound mesh and raises otherwise: on a badly wound import, a uniform pressure would silently become a sign-alternating load.

  • hm.traction(domain, t) — force per unit area in global components. This never touches the normal, so it is the correct choice for a body force such as self weight on a curved roof, where the two are very different loads.

  • hm.moment(domain, m) — moment per unit area, conjugate to the director rotation.

On an edge

hm.edge_traction, hm.edge_moment and hm.edge_pressure apply the same three loads per unit length on the exterior facets a where= predicate selects — the same predicate contract the boundary conditions use:

tip = hm.edge_traction(domain, [0.0, 0.0, -1.0], where=hm.near("x", LENGTH))
lip = hm.edge_moment(domain, [0.0, M, 0.0], where=hm.near("x", LENGTH))

The coefficient stays a CSDL variable, so an edge load is differentiable like any other; the tagged facet measure is structural form data and is fixed at construction.

Concentrated and direct

hm.point_load(domain, at=..., force=..., moment=...) is a consistent point load: the containing cell is located, the state-space basis is evaluated there, and the result is scattered into the right-hand side — the weak form of a Dirac delta, so at need not be a mesh vertex. It is differentiable in force and moment, but not in at or the mesh coordinates. hm.load_vector(domain, vec) supplies a generalized right-hand side directly, in domain.W dof order.

Boundary conditions

Two paths select a region with a coordinate predicate:

  • Penalty (default): a penalty term over the selected facets constrains the masked generalized DOFs. It eliminates no degrees of freedom, and supports several independently masked regions in one form.

  • Strong Dirichlet (hm.clamp(..., method="strong")): standard dirichletbc constraints.

hm.clamp fixes all six DOFs, hm.pin a named subset of ("ux", "uy", "uz", "rx", "ry", "rz"), hm.symmetry applies a plane mask, and hm.gauge adds strong point pins to remove residual null-space modes. hm.near and hm.on_plane build the predicates; any Callable[[ndarray (3, N)], ndarray[bool] (N,)] works.

Prescribed (non-zero) values

Every builder takes value=, which defaults to zero. It accepts a scalar, a (6,) vector over (ux, uy, uz, rx, ry, rz), or a callable that receives (N, 3) coordinates and returns (N, 6):

# a 1 mm enforced settlement, and a linearly varying enforced rotation
settle = hm.pin(domain, where=hm.near("x", L), dofs=("uz",), value=[0, 0, -1e-3, 0, 0, 0])
twist  = hm.pin(domain, where=hm.near("x", L), dofs=("rx",),
                value=lambda x: np.column_stack([np.zeros((len(x), 3)),
                                                 0.01 * x[:, 1], np.zeros((len(x), 2))]))

A prescribed value may not be a csdl.Variable: it is interpolated once into the state space rather than carried as a differentiable solve input, so an enforced displacement cannot be a design variable, while loads and material properties can. Hermit raises a TypeError rather than silently ignoring the dependence. Widening this later is backward-compatible.

Merging

a + b merges two BoundaryConditions on the same domain, and later terms win on overlap: a facet claimed by two regions is governed by the later region’s mask alone, not the union of both masks. That makes hm.clamp(...) + hm.symmetry(...) behave predictably where the two regions touch.

CSDL coupling and adjoints

solve is a CSDL CustomImplicitOperation. Because the shell tangent \(\partial R/\partial w\) is symmetric, a single MUMPS factorization serves both the forward linear solve and the adjoint. The state-based output functions are a CustomExplicitOperation that assembles each scalar form and, in reverse mode, \(\partial(\text{form})/\partial(\text{arg})\).

Reverse-mode sensitivities are available for:

with respect to

how

thickness, density, ABD, \(E\), \(\nu\)

ufl.derivative(form, coefficient)

nodal loads / pressure / moments

linear maps in CSDL

direct load vector

trivial (\(\partial R/\partial b = -I\))

mesh coordinates

ufl.derivative(form, SpatialCoordinate), scattered to node order

See Shape (mesh-coordinate) derivatives for the mesh-coordinate path.

Bibliography

[BLT84]

Ted Belytschko, Jerry I. Lin, and Chen-Shyh Tsay. Explicit algorithms for the nonlinear dynamics of shells. Computer Methods in Applied Mechanics and Engineering, 42(2):225–251, 1984. doi:10.1016/0045-7825(84)90026-4.

[Min51]

R. D. Mindlin. Influence of rotatory inertia and shear on flexural motions of isotropic, elastic plates. Journal of Applied Mechanics, 18(1):31–38, 1951. doi:10.1115/1.4010217.

[Rei45]

Eric Reissner. The effect of transverse shear deformation on the bending of elastic plates. Journal of Applied Mechanics, 12(2):A69–A77, 1945. doi:10.1115/1.4009435.

[SF19]

D. Schöllhammer and T. P. Fries. Reissner–mindlin shell theory based on tangential differential calculus. Computer Methods in Applied Mechanics and Engineering, 352:172–188, 2019. arXiv:1812.05596, doi:10.1016/j.cma.2019.04.018.