Corotational Hinge-based Thin Plates/Shells

Computational modeling and simulation of bending behaviors in thin flexible objects is an active research topic in the graphics community. In 1987, the seminal work [TPBF87] introduced tensorial treatments of the second fundamental form, discretized using the finite difference method on a regular quadrilateral grid, to model cloth and deformed surfaces. Subsequently, particle and mass-spring methods [BHW94, VCMT95, CK02] were explored as alternatives to these complex tensorial treatments, aiming to improve efficiency. However, these methods often compromised physical accuracy, resulting in material parameters that were mesh-dependent and not easily transferable across different mesh topologies. Baraff and Witkin [BW98] utilized the hinge angle to model the bending constraint, similar to the approach in [VCMT95], but they focused on the rest-flat configuration for cloth modeling. The hinge angle is measured on an edge-based hinge stencil, which composes a hinge edge and its two adjacent triangles. Building on this stencil, Grinspun et al. [GHDS03] described the bending energy of the rest-curved thin shell on a discrete differential geometry view. In the same published volume, Bridson et al. [BMF03] offered a productive formulation of the hinge bending model. When the curvature is small, these two models can be related by scaling a coefficient [RLR^∗21, FHXW22]. Both models have been adopted by well-known cloth simulators, C-IPC [LKJ21] and Arcsim [NSO12], respectively, due to their simplicity. However, the edge-based hinge bending model is limited in its ability to capture complete local curvature behavior and suffers from mesh dependency [GGRZ06].

To enhance the consistency of the edge-based hinge bending model, a hinge-averaged shape operator [GSH^∗04] was described to model the bending strain, which has been applied to simulate plasticity and fracture [PNdJO14]. Despite this improvement, the convergence to the ground truth solution remains slow [CSvRV18, GGRZ06]. This issue was effectively addressed by introducing an additional degree of freedom (DoF) for midedge normal rotation to correct the hinge-averaged shape operator, resulting in what is known as the midedge operator [Zor05, GGRZ06]. This operator has been integrated into libshell [CSvRV18], which is based on the shear-rigid Koiter shell model, and is also included in the shearable Cosserat shell model [Wei12]. While the midedge operator improves consistency and convergence in thin-shell simulations, the extra midedge DoFs introduce a higher computational burden. In our models, the extra midedge DoFs are not involved.

To improve the efficiency of the edge-based hinge bending model, the hinge angle is linearized under the quasi-isometry (small in-plane strain) condition, resulting in a constant bending energy Hessian and a linear bending force for rest-flat thin shells, known as the Quadratic Shell model [BWH^∗06b, BWH^∗06a]. This model can be generalized from the Cubic Shell model [GGWZ07], which is suitable for rest-curved thin shell configurations. Although both models offer improved efficiency, the inherent drawbacks of the edge-based hinge bending model [GHDS03] persist. The quadratic biharmonic energy [WBH^∗07] has a triangle-centered stencil and tends to perform better in terms of convergence to ground truth than models on the edge-based stencil, though it is limited to isometric, pure bending of plates. Based on the triangle-centered stencil, we introduce an FVM hinge curvature operator and a smoothed hinge curvature operator derived in the corotational frame [Cri97] to address the limitations of the abovementioned edge-based models.

More recently, Le et al. [LDB^∗23] proposed a second-order Discrete Shells [GHDS03] model, demonstrating superior efficiency from the second-order triangle. Similarly benefited from the second-order tessellation, Löschner et al. [LFFJB24] showcased the effectiveness of a second-order three-director finite element shell with microrotation fields. However, computational modeling with second-order elements requires abandoning the piecewise linear triangle structure. Wen and Barbič [WB23] focused on deriving the KL thin-shell mechanical energy for arbitrary 3D volumetric hyperelastic materials, building their computational model from the foundation in the discrete geometry shell [CSvRV18, Wei12].

Another relevant topic in the computational mechanics community is the concept of a "rotation-free shell", where the shell is characterized without nodal rotational DoF. For those interested, a comparison study by Gärdsback et al. [GT07] provides insight, although it focuses on linear shell analysis. More recently, Zhou et al. [ZS12] extended a rotation-free beam model to a rotation-free shell model. While this approach offers certain efficiency advantages, its complex and laborious boundary condition treatment limits its potential applications. The pursuit of high-fidelity and high-performance (accurate, low computational cost, robustness, low sensitivity to poorly shaped meshes) thin-shell simulations increasingly blurs the lines between different research communities. Our models mainly draw inspiration from both discrete geometry shells [BWH^∗06b, GDP^∗06, GGWZ07, GGRZ06, Wei12] and rotation-free shells [OZ00, ZS12].

To maintain the convergence properties of a linear approach while accommodating arbitrarily large rigid body transformations, the corotational approach is commonly employed to measure pure deformation. Müller et al. [MDM^∗02] first introduced the corotational formulation to handle the geometric non-linearity for deformed body simulations in graphics community. Since then, this formulation has been widely adopted in computer graphics applications for stable and efficient solid simulations [ZSTB10, KKB18].

In thin shell formulations with the corotational approach, Etzmuß et al. [EKS03] extract the rotational component from the deformation gradient for each element and apply it to compute the bending stiffness matrices. In a similar vein to extract the rotation field, Thomaszewski et al. [TWS06] adopt subdivision basis functions [COS00] to improve accuracy in cloth simulation, but this basis function comes with significant computational costs. More recently, a smoothed hinge model [Lia24] based on the corotational formulation has been developed to measure pure deformation in the corotational frame, which undergoes only rigid-body motion, for cloth simulation. However, this model is limited to the rest-flat thin shell configuration. In our work, we provide all formulations by providing curvature operators for both rest-flat and -curved configurations, along with detailed boundary treatments. Leveraging the corotational approach, all our models feature constant bending energy Hessians, which enhances the efficiency and stability for implicit simulations.

The shell stencil in the corotational edge-based hinge bending model is defined by an edge-based stencil, consisting of one edge and its two adjacent triangles (see Figure 1). On the other hand, the shell stencil for the corotational FVM/smoothed hinge bending model uses a triangle-centered stencil, which includes one central triangle and its three neighboring triangles (see Figure 2 and Figure 3). For a given point $\mathbf{x}\in R^{3}$ within one shell stencil, the position in the current configuration is computed as:

$$\mathbf{x}=\mathbf{X}+\mathbf{u},$$

(3)

where $\mathbf{X}\in R^{3}$ is the position in the initial configuration, and $\mathbf{u}\in R^{3}$ is the displacement.

Terminologies. In this context, $(\cdot)^{{e}}$ denotes a vector that collects the quantities associated with a shell stencil. For instance, ${\mathbf{x}}^{e}=\begin{bmatrix}{\mathbf{x}}_{1}^{T}&{\mathbf{x}}_{2}^{T}&{\mathbf{x}}_{3}^{T}&{\mathbf{x}}_{4}^{T}\end{bmatrix}^{T}$ aggregates the current nodal positions for an edge-based stencil, while ${\mathbf{x}}^{e}=\begin{bmatrix}{\mathbf{x}}_{1}^{T}&{\mathbf{x}}_{2}^{T}&{\mathbf{x}}_{3}^{T}&{\mathbf{x}}_{4}^{T}&{\mathbf{x}}_{5}^{T}&{\mathbf{x}}_{6}^{T}\end{bmatrix}^{T}$ aggregates the current nodal positions for a triangle-centered stencil. The notation $(\cdot)_{0}$ denotes the quantity in the initial configuration. The notation $(\cdot)_{ij}$ denotes a directed line segment from point $(\cdot)_{i}$ to point $(\cdot)_{j}$ in any coordinate frame. The tilde $\tilde{(\cdot)}$ represents the quantity defined in the corotational frame ($\tilde{X}-\tilde{Y}-\tilde{Z}$ is the initial corotational frame and $\tilde{x}-\tilde{y}-\tilde{z}$ is the current corotational frame). Quantities in the corotational frame can be transformed from those in the world frame, as illustrated in Appendix A. $(\cdot)_{{c}}$ denotes projection. More graphical illustration can be found in the accompanying Figure 1, Figure 2 and Figure 3.

The terms "EP", "ES", "FP", "FS", "SP", and "SS" are used to distinguish between different bending formulations, when necessary. Specifically, "E" refers to the edge-based hinge, "F" to the FVM hinge, and "S" to the smoothed hinge. "P" indicates a plate (rest-flat configuration), while "S" denotes a shell (rest-curved configuration).

Remarks. Transitioning from simpler cases to more complex scenarios, we provide more derivation details for the bending models based on edge-centered stencils using corotational approach, enabling a seamless generalization of this derivation process to bending models on triangle-centered stencils. To facilitate understanding of the derivation process, we outline several key points.

The kinematics of a shell stencil deformed from its initial configuration to the current configuration can be expressed as $\mathbf{x}^{e}=\mathbf{X}^{e}+\mathbf{u}^{e}$. In the initial configuration, a shell stencil with ${\mathbf{X}}^{e}$ projected onto the tangential plane of the initial corotational frame is referred to as a corotational shell stencil with ${\mathbf{X}}^{e}_{c}$. By projecting the deviation vector $\mathbf{X}^{e}-\mathbf{X}^{e}_{c}$ along the normal direction of the initial corotational frames, and introducing curvature operators constructed from ${\mathbf{X}}^{e}_{c}$, the discretized curvature in the initial configuration can be defined. Under the small (in-plane) strain and curvature assumption, the relative positions of the projected positions ${\mathbf{x}}^{e}_{c}$ in the current corotational frame and ${\mathbf{X}}^{e}_{c}$ in the initial corotational frame are approximately identical. (Another view is that the corotational shell stencil transitions from the initial configuration to the current configuration as ${\mathbf{x}}^{e}_{c}={\mathbf{X}}^{e}_{c}+{\mathbf{u}}^{e}_{c}$, approximately undergoing rigid-body motion). So, the discretized curvature in the current configuration can be defined using the curvature operators of the discretized curvature in the initial configuration. The bending deformation is quantified by the change in curvature between the initial and current configurations. Based on this, the constant bending energy Hessians can be rationally derived. Further details are provided in the subsequent subsections. Our source code (https://github.com/liangqx-hku/libThinPlateShells) is also made available to support practitioners.

For an edge-based stencil based on the triangular mesh (see Figure 1), there is no bending along the hinge edge. In the current configuration, the bend angle $\theta$ is small under the small strain/curvature assumption. The directional curvature $\kappa_{{m}}$ across the hinge edge in the current configuration can be expressed as

$$\kappa_{{m}}=\frac{2sin\theta}{\sqrt{h_{1}^{2}+2h_{1}h_{4}\cos\theta+h_{4}^{2}}}\simeq\frac{2\theta}{h_{1}+h_{4}}$$

(4)

with $\theta\simeq\sin\theta$ and $\cos\theta\simeq 1$. $h_{1}$ and $h_{4}$ are the heights of the triangle $T_{123}$ and $T_{432}$ approximated in the initial configuration, respectively. Under the small curvature assumption, the bend angle in the current configuration can be approximated by

$$\theta=\alpha_{1}+\alpha_{2}\simeq\sin\alpha_{1}+\sin\alpha_{2}\simeq\frac{\mathbf{n}_{\tilde{z}}^{T}({\tilde{\mathbf{x}}}_{1}-{\tilde{\mathbf{x}}}_{P})}{h_{1}}+\frac{\mathbf{n}_{\tilde{z}}^{T}({\tilde{\mathbf{x}}}_{4}-{\tilde{\mathbf{x}}}_{Q})}{h_{4}}.$$

(5)

Here, the perpendicular feet $\tilde{\mathbf{x}}_{P}$ and $\tilde{\mathbf{x}}_{Q}$ are given by

$$\tilde{\mathbf{x}}_{P}=\frac{\lVert\tilde{\mathbf{x}}_{P3}\rVert}{\lVert\tilde{\mathbf{x}}_{23}\rVert}\tilde{\mathbf{x}}_{2}+\frac{\lVert\tilde{\mathbf{x}}_{P2}\rVert}{\lVert\tilde{\mathbf{x}}_{23}\rVert}\tilde{\mathbf{x}}_{3},\ \tilde{\mathbf{x}}_{Q}=\frac{\lVert\tilde{\mathbf{x}}_{Q3}\rVert}{\lVert\tilde{\mathbf{x}}_{23}\rVert}\tilde{\mathbf{x}}_{2}+\frac{\lVert\tilde{\mathbf{x}}_{Q2}\rVert}{\lVert\tilde{\mathbf{x}}_{23}\rVert}\tilde{\mathbf{x}}_{3}.$$

(6)

Under the small strain assumption, these perpendicular feet in Eq. (6) can be approximated as

$$\tilde{\mathbf{x}}_{P}=\frac{\lVert\tilde{\mathbf{X}}_{P3}\rVert}{\lVert\tilde{\mathbf{X}}_{23}\rVert}\tilde{\mathbf{x}}_{2}+\frac{\lVert\tilde{\mathbf{X}}_{P2}\rVert}{\lVert\tilde{\mathbf{X}}_{23}\rVert}\tilde{\mathbf{x}}_{3},\ \tilde{\mathbf{x}}_{Q}=\frac{\lVert\tilde{\mathbf{X}}_{Q3}\rVert}{\lVert\tilde{\mathbf{X}}_{23}\rVert}\tilde{\mathbf{x}}_{2}+\frac{\lVert\tilde{\mathbf{X}}_{Q2}\rVert}{\lVert\tilde{\mathbf{X}}_{23}\rVert}\tilde{\mathbf{x}}_{3},$$

(7)

and $\mathbf{n}_{\tilde{z}}=\mathbf{n}_{c}/\lVert\mathbf{n}_{c}\rVert$ is the direction of the $\tilde{z}$-axis, with $\mathbf{n}_{c}$ expressed as

$$\mathbf{n}_{c}=\frac{\mathbf{x}_{P1}}{\lVert\mathbf{x}_{P1}\rVert}+\frac{\mathbf{x}_{Q4}}{\lVert\mathbf{x}_{Q4}\rVert}.$$

(8)

It should be noted that we use the direction $\mathbf{n}_{\tilde{z}}$, which bisects the bend angle, to approximate the normal direction of the smoothed shell surface in the current configuration. $\mathbf{n}_{\tilde{z}}$ is perpendicular to the hinge edge and points from the shell surface to the mesh. To avoid the numerical issue, when $T_{123}$ and $T_{432}$ are coplanar, $\mathbf{n}_{\tilde{z}}$ aligns with the normal of triangle $T_{123}$ ($T_{432}$).

By substituting Eq. (7) and Eq. (5) into Eq. (4), the directional curvature in the current configuration is obtained as

$$\kappa_{m}={\mathbf{L}}_{m}\mathbf{N}^{T}\tilde{\mathbf{x}}^{e}=\mathbf{L}_{m}\mathbf{d}.$$

(9)

Here ${\mathbf{L}}_{m}=2\mathbf{L}_{\theta}/(h_{1}+h_{4})$ denotes the corotational edge-based hinge curvature operator, where $\mathbf{L}_{\theta}$ is given by

$$\begin{bmatrix}\frac{1}{h_{1}}&-\left(\frac{\lVert\tilde{\mathbf{X}}_{P3}\rVert}{\lVert\tilde{\mathbf{X}}_{23}\rVert h_{1}}+\frac{\lVert\tilde{\mathbf{X}}_{Q3}\rVert}{\lVert\tilde{\mathbf{X}}_{23}\rVert h_{4}}\right)&-\left(\frac{\lVert\tilde{\mathbf{X}}_{P2}\rVert}{\lVert\tilde{\mathbf{X}}_{23}\rVert h_{1}}+\frac{\lVert\tilde{\mathbf{X}}_{Q2}\rVert}{\lVert\tilde{\mathbf{X}}_{23}\rVert h_{4}}\right)&\frac{1}{h_{4}}\end{bmatrix}$$

(10)

and $\mathbf{N}=\mathbf{I}_{4\times 4}\otimes\mathbf{n}_{\tilde{z}}$. $\otimes$ represents the Kronecker product, $\mathbf{I}_{4\times 4}$ is the fourth order identity matrix and the transverse displacement vector of the edge-based stencil $\mathbf{d}=\begin{bmatrix}\tilde{w}_{1}&\tilde{w}_{2}&\tilde{w}_{3}&\tilde{w}_{4}\end{bmatrix}^{T}$ measures the deviations of the nodes away from the corotational shell stencil in the current configuration (see Figure 1).

To conveniently get the derivatives of curvature, we express the transverse displacement vector using the world coordinates, leading to the directional curvature

$$\kappa_{{m}}=\mathbf{L}_{m}\mathbf{N}^{T}(\mathbf{x}^{e}-\mathbf{x}_{c}^{e}).$$

(11)

In the geometry view, $\mathbf{x}^{e}_{c}$ lies in the tangential plane of the current corotational frame. Since $\mathbf{N}$ contains the normal directions of the current corotational frame, it’s obvious that $\mathbf{N}^{T}\mathbf{x}^{e}_{c}=0$. We can obtain that

$$\mathbf{L}_{m}\mathbf{N}^{T}\mathbf{x}^{e}_{c}=\mathbf{0}.$$

(12)

Thus, the curvature of the edge-based shell stencil in the current configuration is

$$\kappa_{{m}}=\mathbf{L}_{m}\mathbf{N}^{T}\mathbf{x}^{e}.$$

(13)

Similarly, the curvature of the edge-based shell stencil in the initial configuration is

$$(\kappa_{{m}})_{0}=\mathbf{L}_{m}\mathbf{d}_{0}=\mathbf{L}_{m}\mathbf{N}_{0}^{T}\mathbf{X}^{e},$$

(14)

where $\mathbf{N}_{0}=\mathbf{I}_{4\times 4}\otimes\mathbf{n}_{\tilde{Z}}$. In fact, $\mathbf{d}_{0}=\mathbf{N}_{0}^{T}\mathbf{X}^{e}$ measures the deviations of the edge-based shell stencil away from the corotational shell stencil in the initial configuration.

The curvature derivation process for the corotational edge-based hinge bending model can also be generalized to the corotational FVM/smoothed hinge bending model. The difference lies in quantifying curvature from the edge-based shell stencil to the triangle-centered shell stencil. For the corotational FVM/smoothed hinge bending model, the direction of the $\tilde{z}$-axis in the current corotational frame is exactly given by the normal of the central triangle $T_{123}$, i.e., $\mathbf{n}_{c}$ in Eq. (8) should be replaced by

$$\mathbf{n}_{c}=\mathbf{x}_{12}\times\mathbf{x}_{13}.$$

(15)

For any point within the shell stencil under the small strain/curvature assumption, the relation

$$\mathbf{n}\mathbf{n}^{T}(\mathbf{x}-\mathbf{x}_{c})=\mathbf{n}\tilde{w}\simeq\mathbf{x}-\mathbf{x}_{c}$$

(19)

holds, and with $\mathbf{L}_{m}\mathbf{N}^{T}\mathbf{x}^{e}_{c}=\mathbf{0}$ in Eq. (12). Therefore, the bending energy can be simplified to

$$\Psi_{b}^{EP}=\frac{A_{\mathcal{E}}}{2}k_{b}(\mathbf{x}^{e})^{T}(\mathbf{L}_{m}^{T}\mathbf{L}_{m}\otimes\mathbf{I})\mathbf{x}^{e}.$$

(20)

This shows that the bending energy is quadratic in terms of the nodal positions, which implies that the Hessian of the bending energy is constant and given by

$$\frac{\partial^{2}\Psi_{b}^{EP}}{\partial\mathbf{x}^{e}{\partial(\mathbf{x}^{e})^{T}}}=k_{b}A_{\mathcal{E}}\mathbf{L}_{m}^{T}\mathbf{L}_{m}\otimes\mathbf{I}.$$

(21)

The gradient of the bending energy, being linear with respect to the nodal positions, is expressed as

$$\frac{\partial\Psi_{b}^{EP}}{\partial\mathbf{x}^{e}}=(\frac{\partial^{2}\Psi_{b}^{EP}}{\partial\mathbf{x}^{e}{\partial(\mathbf{x}^{e})^{T}}})\mathbf{x}^{e}.$$

(22)

It is important to note that quantifying the edge-based curvature operator in the world frame leads to a discrete expression similar to that of the Quadratic Shell model [BWH^∗06b]. However, this formulation overestimates the bending energy by a factor of three, as numerical results will be discussed in detail in Section 6.1. The fundamental difference and accuracy discrepancy has been detailed in Appendix D.

The gradient of the bending energy is then

$$\frac{\partial\Psi_{b}^{ES}}{\partial\mathbf{x}^{e}}={A_{\mathcal{E}}}k_{b}\frac{\partial\epsilon_{b}}{\partial\mathbf{x}^{e}}\epsilon_{b},$$

(25)

where the gradient of the curvature change is

$$\frac{\partial\epsilon_{b}}{\partial\mathbf{x}^{e}}=({\mathbf{x}}^{e})^{T}\frac{\partial\mathbf{N}}{\partial\mathbf{x}^{e}}\mathbf{L}_{m}^{T}+\mathbf{N}\mathbf{L}_{m}^{T}.$$

(26)

Here, ${\partial\mathbf{N}}/{\partial\mathbf{x}^{e}}=\mathbf{I}_{4\times 4}\otimes({\partial\mathbf{n}_{\tilde{z}}}/{\partial\mathbf{x}^{e}})$ with the gradient of the normal ${\partial\mathbf{n}_{\tilde{z}}}/{\partial\mathbf{x}^{e}}$ is detailed in the Appendix B.

Before deriving the Hessian of bending energy, substituting Eq. (24) into Eq. (23) yields

$$\Psi_{b}^{ES}=\Psi_{flat}^{ES}+\Psi_{curved}^{ES},$$

(27)

which decomposes into a flat part

$$\Psi_{flat}^{ES}=\frac{A_{\mathcal{E}}}{2}k_{b}(\mathbf{x}^{e})^{T}(\mathbf{L}_{m}^{T}\mathbf{L}_{m}\otimes\mathbf{I})\mathbf{N}\mathbf{N}^{T}\mathbf{x}^{e}$$

(28)

and a curved part $\Psi_{curved}^{ES}$ is

$$\begin{gathered}\frac{A_{\mathcal{E}}}{2}k_{b}(-2(\mathbf{x}^{e})^{T}(\mathbf{L}_{m}^{T}\mathbf{L}_{m}\otimes\mathbf{I})\mathbf{N}\mathbf{N}_{0}^{T}\mathbf{X}^{e}+(\mathbf{X}^{e})^{T}(\mathbf{L}_{m}^{T}\mathbf{L}_{m}\otimes\mathbf{I})\mathbf{N}_{0}\mathbf{N}_{0}^{T}\mathbf{X}^{e}).\end{gathered}$$

(29)

The flat part in Eq. (27) is identical to the bending energy of the corotational edge-based hinge thin plate in Eq. (18), similar note can be found in the Cubic Shell paper [GGWZ07]. Consequently, the simplified bending energy Hessian from the rest-flat version in Eq. (21) can be used to perform a single bending energy Hessian assembly for the rest-curved version.

The Clamped Boundary Condition. For a boundary edge stencil with a clamped boundary condition, a symmetric virtual transverse displacement is applied, corresponding to the virtual node $L^{\prime}$

$$\tilde{w}_{L^{\prime}}=\tilde{w}_{L},$$

(30)

which ensures the preservation of the normal direction perpendicular to the boundary edge stencil. When this condition is applied to the boundary edge stencil, the corotational edge-based hinge curvature operator $\mathbf{L}_{m}$ is modified to

$$\mathbf{L}_{m}^{{}^{\prime}}=\frac{2}{h_{1}}\begin{bmatrix}\frac{1}{h_{1}}&-\frac{\lVert\tilde{\mathbf{X}}_{P3}\rVert}{\lVert\tilde{\mathbf{X}}_{23}\rVert h_{1}}&-\frac{\lVert\tilde{\mathbf{X}}_{P2}\rVert}{\lVert\tilde{\mathbf{X}}_{23}\rVert h_{1}}&0\end{bmatrix}$$

(31)

This zero-slope condition is also applicable in symmetric structural finite element analysis (FEA).

The Free Boundary Condition. For a free boundary edge, there is no bending energy in the boundary edge stencil, resulting in zero entries for both the gradient and Hessian related to this boundary (zero-curvature condition).

In conclusion, the modification of the curvature operator according to the missing nodes enables the condensation of boundary conditions at the element level. Boundary conditions are crucial in engineering shell simulations to achieve accurate results. However, their impact on the visual effects in computer animation is often negligible.

Smooth Settings. The curvature of the central triangle in the initial configuration can be evaluated using the finite volume method (FVM) within the corotational frame. The generalized curvature vector

$$\bm{\kappa}_{0}=\frac{1}{A}\int_{A}\left(\begin{array}[]{c}\frac{\partial^{2}\tilde{w}_{0}}{\partial\tilde{X}^{2}}\\ \frac{\partial^{2}\tilde{w}_{0}}{\partial\tilde{Y}^{2}}\\ 2\frac{\partial^{2}\tilde{w}_{0}}{\partial\tilde{X}\partial\tilde{Y}}\end{array}\right)dA$$

(32)

is derived through variation the virtual work principle governing the bending deformation. Here $A$ is the area of the central triangle in the initial configuration. By applying the Stokes’s theorem, Eq. (32) can be rewritten as

$$\bm{\kappa}_{0}=\frac{1}{A}\int_{\Gamma}\mathbf{T}\binom{\frac{\partial\tilde{w}_{0}}{\partial\tilde{X}}}{\frac{\partial\tilde{w}_{0}}{\partial\tilde{Y}}}d\Gamma,$$

(33)

where $\Gamma$ is the boundary of the central triangle in the initial configuration, and $\mathbf{T}$ is defined as

$$\mathbf{T}=\left[\begin{array}[]{cc}m_{s}&0\\ 0&m_{t}\\ m_{t}&m_{s}\end{array}\right],$$

(34)

with the normalized normal $\tilde{\mathbf{m}}=\left[\begin{array}[]{cc}m_{s}&m_{t}\end{array}\right]^{T}$ outward to the boundary $\Gamma$ surrounding the central triangle in the $\tilde{X}-\tilde{Y}$ plane of the initial corotational frame. ${s}$-axis is tangential to the boundary $\Gamma$ and ${t}$-axis is perpendicular to the ${s}$-axis. The gradient of the transverse displacement $\tilde{w}$ can be transformed from the directional gradient using

$$\binom{\frac{\partial\tilde{w}_{0}}{\partial\tilde{X}}}{\frac{\partial\tilde{w}_{0}}{\partial\tilde{Y}}}=\left[\begin{array}[]{cc}m_{{{s}}}&-m_{{t}}\\ m_{{t}}&m_{{s}}\end{array}\right]\binom{\frac{\partial\tilde{w}_{0}}{\partial t}}{\frac{\partial\tilde{w}_{0}}{\partial s}},$$

(35)

where ${\partial\tilde{w}_{0}}/{\partial t}$ and ${\partial\tilde{w}_{0}}/{\partial s}$ are respectively the rotation angle about the ${s}$-axis and ${t}$-axis.

Discrete Settings. Given that there is no bending in each hinge edge (see Figure 2), i.e., ${\partial\tilde{w}}/{\partial t}\neq 0$ and ${\partial\tilde{w}}/{\partial s}=0$ in Eq. (35), the generalized curvature of the central triangle in the initial configuration becomes

$$\bm{\kappa}_{0}=\sum_{i=1}^{3}\left(\begin{array}[]{c}(m_{{s}})_{i}^{2}\\ (m_{{t}})_{i}^{2}\\ 2(m_{{s}})_{i}(m_{{t}})_{i}\end{array}\right)(\kappa_{i})_{0},$$

(36)

where $(\kappa_{i})_{0}=({\partial\tilde{w}}/{\partial t})_{i}l_{i}/A$ is the directional curvature, with $l_{i}$ being the length of the hinge edge in the initial configuration. It can be seen that the FVM approach superpositions the directional curvatures of the central triangle onto the generalized curvature in the corotational frame. To discrete the curvature and consider the adjacent triangles’ contribution to the common edge, several discretized schemes can be found in [OZ00, GSH^∗04].

In our approach, we first extend the Eq. (4) to express the directional curvatures, $\bm{\kappa}_{0}$ can then be rewritten as

$$\bm{\kappa}_{0}=\mathbf{R}(\bm{\kappa}_{p})_{0}=\mathbf{R}\left(\begin{array}[]{c}\frac{2(\theta_{1})_{0}}{h_{1}+h_{4}}\\ \frac{2(\theta_{2})_{0}}{h_{2}+h_{5}}\\ \frac{2(\theta_{3})_{0}}{h_{3}+h_{6}}\end{array}\right),$$

(37)

where $\bm{\kappa}_{p}$ denotes the directional curvature vector, $\theta_{i}\ (i=1,2,3)$ is the bend angle across each edge, $h_{i}\ (i=1,2,\cdots,6)$ is the triangle height in the initial configuration (see Figure 9), and the transform matrix $\mathbf{R}$ is

$$\mathbf{R}=\left[\begin{array}[]{ccc}(m_{{s}})^{2}_{1}&(m_{{s}})_{2}^{2}&(m_{{s}})_{3}^{2}\\ (m_{{t}})_{1}^{2}&(m_{{t}})_{2}^{2}&(m_{{t}})_{3}^{2}\\ 2(m_{{s}})_{1}(m_{{t}})_{1}&2(m_{{s}})_{2}(m_{{t}})_{2}&2(m_{{s}})_{3}(m_{{t}})_{3}\end{array}\right].$$

(38)

Finally, under the small strain/curvature assumption, each directional curvature in Eq. (37) is linearized using the corotational edge-based hinge approach described in Section 4.3.

In deriving the derivatives of the bending energy of corotational FVM hinge thin plate, we also obtain the constant bending energy Hessian corresponding to the corotational edge-based hinge thin plate in Section 4.3,

$$\frac{\partial^{2}\Psi_{b}^{FP}}{\partial\mathbf{x^{\mathrm{e}}}{\partial(\mathbf{x}^{e})^{T}}}=A(\mathbf{L}_{p}^{T}\mathbf{D}_{\mathrm{b}}^{p}\mathbf{L}_{p}\otimes\mathbf{I}).$$

(41)

The operation in Eq. (22) can similarly be employed to compute the linear gradient of the corotational FVM hinge thin plate model.

The bending energy Hessian remains the same as in the corotational FVM hinge thin plate model, for reasons similar to those discussed in the corotational edge-based hinge bending model in Section 4.3.

In the smoothed hinge bending model, the generalized curvature vector

$$\bm{\kappa}_{0}=\left(\begin{array}[]{c}\frac{\partial^{2}\tilde{w}_{0}}{\partial\tilde{X}^{2}}\\ \frac{\partial^{2}\tilde{w}_{0}}{\partial\tilde{Y}^{2}}\\ 2\frac{\partial^{2}\tilde{w}_{0}}{\partial\tilde{X}\partial\tilde{Y}}\end{array}\right)$$

(46)

is derived directly from the area integration in Eq. (32), rather than transforming the area integration to line integration using integration by parts in Eq. (33). The generalized curvature vector in Eq. (46) can be discretized as

$$\bm{\kappa}_{0}=\mathbf{L}_{g}\mathbf{d}_{0},$$

(47)

where $\mathbf{d}_{0}=\begin{bmatrix}(\tilde{w}_{1})_{0}&(\tilde{w}_{2})_{0}&(\tilde{w}_{3})_{0}&(\tilde{w}_{4})_{0}&(\tilde{w}_{5})_{0}&(\tilde{w}_{6})_{0}\end{bmatrix}^{T}$ is the vector of transverse displacements in the initial configuration. In the smoothed hinge model, a quadratic fitting surface is used to smooth the triangle-centered stencil (see Figure 3) and $\mathbf{L}_{g}$ is denoted as corotational smoothed hinge curvature operator. So the transverse displacement field on the shell surface can be interpolated as

$$(\tilde{w})_{0}=c_{1}+c_{2}\tilde{X}+c_{3}\tilde{Y}+c_{4}\tilde{X}^{2}/2+c_{5}\tilde{Y}^{2}/2+c_{6}\tilde{X}\tilde{Y}/2,$$

(48)

where $c_{i}\ (i=1,\cdots,6)$ are constants. The first three coefficients correspond to rigid-body motion, while the last three correspond to curvature within the quadratic surface. By introducing Eq. (48) into Eq. (46) with Eq. (47), the generalized curvature vector can be rewritten as $\bm{\kappa}_{0}=\mathbf{L}_{g}\mathbf{d}_{0}=\begin{bmatrix}c_{4}&c_{5}&c_{6}\end{bmatrix}^{T}$. Its three components can be expressed by converting the three directional curvature constituent parts from the directional curvature vector $(\bm{\kappa}_{p})_{0}$. By recalling Eq. (9), Eq. (13) and Eq. (40), $(\bm{\kappa}_{p})_{0}$ can be rewritten as $(\bm{\kappa}_{p})_{0}=\mathbf{L}_{p}\mathbf{d}_{0}$. Consequently,

$$(\bm{\kappa}_{p})_{0}=\mathbf{L}_{p}\mathbf{d}_{0}=\mathbf{L}_{p}\mathbf{C}_{p}\left(\begin{array}[]{c}c_{4}\\ c_{5}\\ c_{6}\end{array}\right)=\mathbf{L}_{p}\mathbf{C}_{p}\mathbf{L}_{g}\mathbf{d}_{0},$$

(49)

where $\mathbf{C}_{p}$ is portrayed as a matrix $\left[\begin{array}[]{cccccc}\tilde{X}_{1}^{2}/2&\tilde{Y}_{1}^{2}/2&\tilde{X}_{1}\tilde{Y}_{1}/2\\ \vdots&\vdots&\vdots\\ \tilde{X}_{6}^{2}/2&\tilde{Y}_{6}^{2}/2&\tilde{X}_{6}\tilde{Y}_{6}/2\end{array}\right]$. We then explicitly get the corotational smoothed hinge curvature operator

$$\mathbf{L}_{g}=(\mathbf{L}_{p}\mathbf{C}_{p})^{-1}\mathbf{L}_{p}$$

(50)

to express the generalized curvature of the smoothed hinge bending model in the initial configuration as

$$\bm{\kappa}_{0}=\mathbf{L}_{g}\mathbf{N}_{0}^{T}\mathbf{X}^{e}.$$

(51)

It should be mentioned that the constant curvature of the quadratic surface can also be computed by the least square method to solve a sixth order linear system for each curvature operator. In practice, we find that our presented method is more stable compared to this approach. Another explicit formulation is presented in a follow-up note by Reisman et al. [RGZ07], which continues the work of Grinspun et al. [GGRZ06], they highlight certain drawbacks of the curvature operator fitted on the quadratic surface. Notably, near-conic degenerate configurations can lead to numerical instabilities and challenges in boundary condition treatment. Thanks to our small strain/curvature assumption within the corotational framework, the bending energy Hessian matrix can be assembled once using the initial geometric data. A good mesh initialization can effectively mitigate the numerical issue. In our numerous numerical exercises, even with no special initialization, we find that our model remains stable. More detailed results can be found in the Section 6. Additionally, we provide a concise treatment of boundary conditions on our corotational smoothed hinge curvature operator, which will be elaborated upon later.

The destination of the corotational smoothed hinge bending energy Hessian

$$\frac{\partial\Psi_{{b}}^{SP}}{\partial\mathbf{x}^{\mathrm{e}}}=A(\mathbf{L}_{g}^{T}\mathbf{D}_{{b}}\mathbf{L}_{g}\otimes\mathbf{I})$$

(53)

is also constant, and the linear gradient can also be computed by the approach in Eq.(22).

The gradient of the bending energy of corotational smoothed hinge thin shell is

$$\frac{\partial\Psi_{\mathrm{b}}^{SS}}{\partial\mathbf{x}^{\mathrm{e}}}=A\frac{\partial(\bm{\varepsilon}_{b}^{SS})^{T}}{\partial\mathbf{x}^{e}}\mathbf{D}_{{b}}\bm{\varepsilon}_{b}^{SS},$$

(56)

where the gradient of the curvature change vector is

$$\frac{\partial(\bm{\varepsilon}_{b}^{SS})^{T}}{\partial\mathbf{x}^{e}}=({\mathbf{x}}^{e})^{T}\frac{\partial\mathbf{N}}{\partial\mathbf{x}^{e}}\mathbf{L}_{g}^{T}+\mathbf{N}\mathbf{L}_{g}^{T}.$$

(57)

The bending energy Hessian is the same as that of the corotational smoothed hinge thin plate model. The rationale behind this similarity is consistent with the explanation provided for the corotational edge-based hinge thin shell model.

The Clamped Edge Boundary Condition. To integrate the zero-slope condition, using Eq. (30), the modified curvature with a clamped edge can be written as

$$\begin{split}\bm{\kappa}_{clamp}=&(\mathbf{L}_{gL}+\mathbf{L}_{g(L+3)})\tilde{w}_{L}+\mathbf{L}_{gM}\tilde{w}_{M}\\ &+\mathbf{L}_{gN}\tilde{w}_{N}+\mathbf{L}_{g(N+3)}\tilde{w}_{N+3}+\mathbf{L}_{g(M+3)}\tilde{w}_{M+3}.\end{split}$$

(59)

The Free Edge Boundary Condition. For the zero-curvature condition, the virtual transverse displacement of the missing node should satisfy

$$\tilde{w}_{L^{\prime}}=2\tilde{w}_{P}-\tilde{w}_{L},$$

(60)

where $\tilde{w}_{P}=(\tilde{w}_{M}+\tilde{w}_{N})/2$ is the transverse displacement of the middle point $P$ of the boundary edge $NM$. We can deduce the modified curvature with a free edge as follows

$$\begin{split}\bm{\kappa}_{free}=&(\mathbf{L}_{gN}+\mathbf{L}_{g(L+3)})\tilde{w}_{N}+(\mathbf{L}_{gM}+\mathbf{L}_{g(L+3)})\tilde{w}_{M}\\ &+(\mathbf{L}_{gL}-\mathbf{L}_{g(L+3)})\tilde{w}_{L}+\mathbf{L}_{g(N+3)}\tilde{w}_{N+3}+\mathbf{L}_{g(M+3)}\tilde{w}_{M+3}.\end{split}$$

(61)

If a boundary triangle includes one more boundary edge, the operations outlined in Eq. (59) and Eq. (61) can be superimposed.

As shown in Figure 5, both the EP and ES models pass the test only on the equilateral triangular mesh. However, the Quadratic Shell and Cubic Shell models fail all tests, even on the equilateral triangular mesh. The primary issue, as mentioned in Appendix D, is that the bending energy formulation in the Quadratic and Cubic Shell models is three times higher than ours. The MidedgeAve model also fails on the equilateral triangular mesh due to numerical issues resulting in "NaN" values. On the regular triangular mesh, our FP, FS, SP, and SS models perform slightly better than the MidedgeSin and MidedgeTan formulations. On the irregular mesh, MidedgeSin and MidedgeTan exhibit better consistency than other methods. Among our models, SP and SS outperform FP and FS, which have comparable performance to the MidedgeAve formulation.

The data from the last column in Table 1 highlights that our FP, FS, SP, and SS models are nearly twice as fast as the MidedgeSin and MidedgeTan formulations in one linear system solve on the equilateral triangular mesh. Additionally, the EP/ES and Quadratic/Cubic Shell models in the edge-based stencil demonstrate exceptional speed. The numerical performance differences across various mesh tessellations primarily arise from the curvature operators used in these formulations.

Cantilever Subjected to End Shear Force. In this test, a flat plate of dimensions $10\times 1$ is subjected to an end shear force, applied as concentrated loads of equal magnitude $F=4/3$ distributed across the nodes on the right side, as illustrated in Figure 6. The concentrated loads are along the $Z$-axis. The geometry is discretized into 51 nodes, with two adjacent rows at the left end clamped to enforce the hard constraints. The material parameters are $E=1.2\times 10^{6}$, $\nu=0.1$, and $h=0.1$. As summarized in Table 1, our FP, FS, SP, and SS models demonstrate superior accuracy compared to others. The MidedgeTan and MidedgeSin formulations, which offer the second-highest accuracy, are nearly five times slower than our FP, FS, SP, and SS models, thanks to their constant bending energy Hessians. Our EP and ES models rank third in accuracy, outperforming the MidedgeAve model. While the Quadratic and Cubic Shell models are very fast, they produce smaller deflections due to the overestimation of bending rigidity. If we scale down the bending rigidity of the Quadratic and Cubic Shell models, they can predict deflections comparable to those of our EP model. Among our models, EP, FP and SP, specifically designed for rest-flat shell configurations, perform slightly better in terms of accuracy and speed compared to ES, FS and SS, which can handle both rest-flat and rest-curved shells.

Hemispherical Shell Subjected to Alternating Radial Forces. To test the performance of the rest-curved shell models, we simulate a hemispherical shell with radius $R=10$ and an $18^{\circ}$ circular cutout at the pole. The shell is subjected to alternating radial point forces of $P=200$ at $90^{\circ}$ intervals (as shown in Figure 6). Two point forces along the $X$-axis induce compression, while two along the $Y$-axis induce tension. To minimize boundary condition effects across different formulations, the entire shell structure is analyzed instead of only a quarter section. Boundary conditions are applied as follows: for nodes lying in the $Y-Z$ plane, the $X$-direction DoFs are fixed; for nodes in the $X-Z$ plane, the $Y$-direction DoFs are fixed. Additionally, for nodes on the top circular cut that lie in the $Y-Z$ plane, the $Z$-direction DoFs are constrained to ensure equivalence with the benchmark case provided in [SLL04]. The shell geometry is discretized using 1088 nodes. The material parameters are $E=6.825\times 10^{7}$, $\nu=0.3$, and $h=0.04$. As shown in Table 1, our ES model outperforms the Cubic Shell models in terms of accuracy. However, the ES model still deviates more from the reference solution. While the MidedgeTan and MidedgeSin models provide more accurate results overall, our FS and SS models deliver competitive accuracy with nearly four times the computational speed of the MidedgeTan and MidedgeSin models.

We introduce two challenging codimensional simulation benchmarks from C-IPC [LKJ21] to demonstrate the robustness and stability of our models. The first test case involves a flat geometry, which is well-suited to all of our models, while the second features a cylindrical geometry, requiring compatible models for rest-curved meshes. Both cases use material properties corresponding to cotton from C-IPC, with $E=0.8\text{MPa}$, $\nu=0.243$, and a cloth density of $472.6\,\text{kg/m}^{3}$. The simulations are time-stepped at $\Delta t=0.04\text{s}$, also consistent with the solver settings from the C-IPC study. Simulation videos can be found in the supplementary materials.

Cloth on Rotating Sphere. This test evaluates the robustness of our formulations (EP, ES, FP, FS, SP and SS) under extreme stress-test conditions, such as tight wrinkling, friction, and contact processing, following the setup used in prior research [BFA02, LKJ21]. A square cloth with 85K nodes (strain limiting up to 1.0608) is dropped onto a sphere and floor, both having a friction coefficient $\mu=0.4$. As the sphere rotates, friction draws the cloth inward, creating a complex structure of wrinkles and folds, effectively capturing fine details of the cloth’s behaviour (see Figure 7).

Twisted Cylinder. In this test, we simulate a cotton cylinder (1m width and 0.25m radius) with 88K nodes. The thickness offset is set to 1.5mm to account for geometric thickness, and the IPC [LFS^∗20] contact force is began at a threshold distance of 1mm. The cylinder is simultaneously twisted at a rate of $72^{\circ}/s$ while the two sides are brought together at 5mm/s. Gravity is excluded from the simulation to prevent sagging. As illustrated frame in Figure 8, global wrinkling and folding effects emerge as the cylinder is deformed, showcasing the ability of our models (ES, FS and SS) to handle rest-curved geometry robustly. It’s worthy to mention that SS generate 20 wrinkles, but ES and FS both give 19 waves.