research-article

Open access

Analytic rotation-invariant modelling of anisotropic finite elements

Authors:

Huancheng Lin,

Floyd M. Chitalu,

Taku KomuraAuthors Info & Claims

ACM Transactions on Graphics, Volume 43, Issue 5

Article No.: 157, Pages 1 - 20

https://doi.org/10.1145/3666086

Published: 09 August 2024 Publication History

PDF eReader

Abstract

Anisotropic hyperelastic distortion energies are used to solve many problems in fields like computer graphics and engineering with applications in shape analysis, deformation, design, mesh parameterization, biomechanics, and more. However, formulating a robust anisotropic energy that is low order and yet sufficiently non-linear remains a challenging problem for achieving the convergence promised by Newton-type methods in numerical optimization. In this article, we propose a novel analytic formulation of an anisotropic energy that is smooth everywhere, low order, rotationally invariant, and at least twice differentiable. At its core, our approach utilizes implicit rotation factorizations with invariants of the Cauchy-Green tensor that arises from the deformation gradient. The versatility and generality of our analysis is demonstrated through a variety of examples, where we also show that the constitutive law suggested by the anisotropic version of the well-known As-Rigid-As-Possible energy is the foundational parametric description of both passive and active elastic materials. The generality of our approach means that we can systematically derive the force and force-Jacobian expressions for use in implicit and quasistatic numerical optimization schemes, and we can also use our analysis to rewrite, simplify, and speed up several existing anisotropic and isotropic distortion energies with guaranteed inversion safety.

1 Introduction

Simulating elastic materials that exhibit directional effects requires the use of anisotropic (i.e., transversely isotropic) constitutive models that describe this rotationally invariant but direction-dependent behavior. However, existing expressions of these anisotropic energies are highly non-linear, which precludes inversion safety when simulating with finite elements (FEs) due to the occurrence of spurious stable states under large deformations.

A notable exception is the recent \(I_4\)-invariant-based expression of anisotropic elastic energy by Kim et al. [2019], which is appealing for its seemingly simple and low-order invariant-expressible form that can be optimized efficiently. Notably the formulation subsumes a rotation-variant polar factorization, and the gradient (force) can become singular. Polar factorization hinders proper analytic expression of the energy gradient due to the occurrence of a fourth-order rotation-gradient tensor, where the rotation matrix from which this tensor is defined will be computed numerically, for example, via singular value decomposition (SVD).¹ Force singularity will halt a simulation, which happens when an FE enters a physically degenerate configuration such as being collapsed to a point, line, or undergoing inversion with large deformations. Thus, while appealing for its simplicity, the \(I_4\)-invariant-based formulation of anisotropic elastic energy lacks sufficient analytic description due to a fundamental reliance on iterative/numerical procedures to compute polar factorization and can suffer from force singularities, which are the two core problems that we aim to address in this article.

We propose to analyze and extend the \(I_4\)-invariant-based formulation of anisotropic elasticity by using analytic rotations [Lin et al. 2022], where we show that the complete Hessian expression of this invariant can be obtained in proper closed form. Standard Cauchy-Green invariants are adopted to evaluate rotation factors as in the work of Lin et al. [2022], from which we also identify new invariants that are required to arrive at a novel (and fully differentiable) expression of \(I_4\). We then derive a closed-form expression of its Hessian for optimization with Newton-type solvers that require second derivatives. We also propose two solutions to eliminate singularities via mollification or regularization, which are designed to ensure that forces are gracefully bounded, and where the energy function is at least \(C^1\) (and piecewise \(C^2\)) continuous for all geometric configurations of an FE. We thus extend and supplement the definition of any energy written in terms of the \(I_4\)-invariant by using these solutions for addressing the incompleteness of the analytic description and singularities.

Our analysis presents numerous applications and benefits by providing a general approach with which the various forms of anisotropy (cf. Figure 1) can be expressed and simulated. A simple, fast, and robust anisotropic As-Rigid-As-Possible (ARAP) energy is made possible with our \(I_4\) Hessian expression by addressing singularities. We in turn use these expressions to show that this anisotropic ARAP energy is a ‘baseline’ parametric description of both passive and active elastic constitutive models in invariant-expressible forms, much like corotational elasticity is largely established as the prototypical rotation-invariant model of isotropic elasticity. An even deeper implication of our analysis is that our \(I_4\)-invariant formulation can be used to rewrite several common anisotropic and isotropic energies including anisotropic ARAP [Alexa et al. 2000; Kim et al. 2019; Sorkine and Alexa 2007], the Baraff-Witkin cloth model [Kim 2020], and shape targeting for simulating active elasticity models [Ichim et al. 2017; Klár et al. 2020a], in lower-order form, with the added benefit of faster performance in terms of solver convergence rates. Orthotropic- and generally anisotropic materials (see Figure 1), which exhibit unique response along prescribed directions, can likewise be expressed and thus understood in terms of our novel formulation of transverse isotropy.

Fig. 1.

A summary of our primary contributions is as follows:

—

A closed-form expression of the anisotropic invariant \(I_{4}\) to enable proper analytic description of its gradients;

—

an approach to identify and resolve singularities associated with \(I_{4}\) for ensuring numerical stability of energies based on this invariant; and

—

novel expressions for rewriting several common hyperelastic energies in simpler and faster invariant-expressible forms.

2 Background

In this section, we review the formulations of invariants that encompass all possible isotropic and anisotropic hyperelastic² FE energies before presenting our \(I_4\)-invariant formulation. For notation, we will represent scalars as italic letters (\(J\)), vectors as bold lowercase letters (\(\mathbf {x}\)), and matrices as bold uppercase letters (\(\mathbf {M}\)).

Our starting point is the recently proposed fundamental set of strain measures from Smith et al. [2019]:

\begin{equation} I_1 = \text{tr}(\mathbf {S}), \quad I_2 = \Vert \mathbf {S}\Vert _F^2, \quad I_3 = \det (\mathbf {S}), \end{equation}

(1)

which are isotropic invariants of the deformation gradient \(\mathbf {F}=\mathbf {R}\mathbf {S}\) admitting polar factorization to obtain the scaling \(\mathbf {S}= \mathbf {R}^T\mathbf {F}\) and a rotation factor \(\mathbf {R}\). This \(\mathbf {F}\) is a matrix encoding the change-of-shape within an FE and is used to describe local deformation (cf. Figure 2). The strain measures are designed to expand the space of invariant-expressible isotropic hyperelastic FE energies into the inversion-safe regime, leading to a systematic method for obtaining analytic eigensystems for minimization with projected Newton solvers. Building on Equation (1), Kim et al. [2019] extend the fundamental set by advocating

\begin{equation} I_4 = \mathbf {a}^T\mathbf {S}\mathbf {a} \end{equation}

(2)

as the strain measure expanding the aforementioned energy space with inversion-safe anisotropic energies. The vector \(\mathbf {a}\) prescribes the direction of the anisotropic force, which also represents the locally preferred direction of material response to deformation.

Fig. 2.

Our focus is expanding the broader applicability of \(I_4\) by recasting the formulation of Kim et al. [2019] into the framework of analytic rotations [Lin et al. 2022] to address several practical challenges associated with this invariant for arriving at a general rotation-invariant model of anisotropy. For brevity, we refer the readers to the work of Bonet and Wood [2008], Kim and Eberle [2022], and Sifakis and Barbic [2012], as well as that of Marsden and Hughes [1994], for an extensive treatment of background literature and theory about the pertinent topic of continuum mechanics in computer graphics and engineering for simulating hyperelastic materials with FEs.

2.1 The Derivatives of Arbitrary I₄ Energies

Analyzing and simulating the deformation behavior of a given elastic material will require the analytic expressions for the stress and Hessian of the potential energy function associated with this material.³ These stress and Hessian expressions are essential for determining the spatial derivatives of the energy. Thus, for a given set of invariants \(\mathcal {I}\), an energy \(\Psi (\mathcal {I})\) that is expressed in terms of these invariants will have first- and second-order spatial derivatives of the form

\begin{align} \frac{\partial \Psi }{\partial \mathbf {x}} &= \frac{\partial \Psi }{\partial \mathbf {F}} :\frac{\partial \mathbf {F}}{\partial \mathbf {x}}, \end{align}

(3)

\begin{align} &= \sum _{i \in \mathcal {I}}\left(\frac{\partial \Psi }{\partial i} \frac{\partial i}{\partial \mathbf {F}}\right):\frac{\partial \mathbf {F}}{\partial \mathbf {x}}, \end{align}

(4)

and

\begin{align} \frac{\partial ^2 \Psi }{\partial \mathbf {x}^2} &= \frac{\partial \mathbf {F}}{\partial \mathbf {x}}: \frac{\partial ^2 \Psi }{\partial \mathbf {F}^2} : \frac{\partial \mathbf {F}}{\partial \mathbf {x}}, \end{align}

(5)

\begin{align} &=\frac{\partial \mathbf {F}}{\partial \mathbf {x}}: \frac{\partial }{\partial \mathbf {F}}\left(\sum _{i \in \mathcal {I}}\left(\frac{\partial \Psi }{\partial i} \frac{\partial i}{\partial \mathbf {F}}\right) \right) : \frac{\partial \mathbf {F}}{\partial \mathbf {x}}, \end{align}

(6)

where \({\partial \Psi }/{\partial \mathbf {F}}\) and \({\partial ^2 \Psi }/{\partial \mathbf {F}^2}\) are the first Piola-Kirchhoff stress (PK1) and Hessian, respectively. It is from this PK1 and Hessian that one obtains the energy gradient (force) \(\partial \Psi /\partial \mathbf {x}\) and force-Jacobian \(\partial ^2 \Psi /\partial \mathbf {x}^2\) (via tensor contraction ‘\(:\)’ with the change-of-basis tensor \({\partial \mathbf {F}}/{\partial \mathbf {x}}\)) as partial derivatives of the scalar potential \(\Psi\) w.r.t. positions \(\mathbf {x}\). These derivatives are most commonly used for time integration in FE simulations.

However, an energy expression that is written purely in terms of \(I_4\) (i.e., \(\Psi (\mathcal {I} = \lbrace {I_4\rbrace }) \rightarrow \Psi _{4}\)) will have a PK1 of the form (cf. Equation (28) in the work of Kim et al. [2019]):

\begin{align} \frac{\partial \Psi _4}{\partial \mathbf {F}} &= \frac{\partial \Psi _4}{\partial {I}_4} \frac{\partial I_4}{\partial \mathbf {F}}, \end{align}

(7)

\begin{align} &= \frac{\partial \Psi _4}{\partial {I}_4} \left(\frac{\partial \mathbf {R}}{\partial \mathbf {F}} : \mathbf {F}\mathbf {A}+ \mathbf {R}\mathbf {A}\right), \end{align}

(8)

where \(\mathbf {A}= \mathbf {a}\mathbf {a}^T\). Equation (7) presents two impeding difficulties. The first is that it contains a fourth-order rotation gradient tensor \(\tfrac{\partial \mathbf {R}}{\partial \mathbf {F}}\), where the rotation \(\mathbf {R}=\mathbf {U}\mathbf {V}^T\) is assumed to come from a rotation-variant polar decomposition that is computed numerically (e.g., via SVD \(\mathbf {F}=\mathbf {U}\Sigma \mathbf {V}^T)\) [Twigg and Kačić-Alesić 2010; Sorkine-Hornung and Rabinovich 2016; Higham 2008; Irving et al. 2004; McAdams et al. 2011b]. Thus, proper closed-form description of this rotation gradient is at best unwieldy, which also hinders exact derivation of the gradient and Hessian (see also the work of Gao et al. [2009] and Papadopoulo and Lourakis [2000] for examples of numerical approximations). The second challenge relates to the behavior of the energy, which exhibits singularities (infinitely large forces) that can make minimization intractable. This issue occurs, for example, during FE inversion when two singular values \(\sigma _i\) of \(\mathbf {F}\) are equal in magnitude but with opposite signs, where \(\textrm {diag}(\Sigma) = (\sigma _1, \ldots , \sigma _d)\).

2.2 Inversion-Safe Anisotropic ARAP Energy

Faced with these issues, Kim et al. [2019] analyze several existing anisotropic energies to construct an inversion-safe energy of the form

\begin{equation} \Psi _{AA} = \frac{\mu }{2}\left(\sqrt {I_5} - \mathcal {S}(I_4)\right)^2. \end{equation}

(9)

To our knowledge, Equation (9) is the only anisotropic ARAP model that is inversion safe, which Kim et al. [2019] define as the property of having a single stable state under large (e.g., flip-inducing) deformations of FEs. The step function \(\mathcal {S}(x)\) detects element inversion, which is based on the sign information of the singular values of \(\mathbf {F}\) via \(I_4\). Its effect is to mask the high-order rotation-gradient tensor from the PK1

\begin{equation} \frac{\partial \Psi _{AA}}{\partial \mathbf {F}} = \mu \left(\sqrt {I_5} - \mathcal {S}(I_4)\right)\left[\frac{1}{2\sqrt {I_5}}\mathbf {F}\mathbf {A}- \delta (I_4)\frac{\partial I_4}{\partial \mathbf {F}}\right], \end{equation}

(10)

and the respective Hessian, where \(\delta (I_4) = \partial \mathcal {S}(I_4)/\partial \mathbf {F}\) is the Dirac delta function that ensures that forces no longer explode to infinity. Anisotropy of material response is determined by the higher-order invariant

\begin{equation} I_5 = \Vert \mathbf {F}\mathbf {a}\Vert ^2_2 = \mathbf {a}^T\mathbf {C}\mathbf {a}= \mathrm{tr}(\mathbf {C}\mathbf {A}), \end{equation}

(11)

where this response is along the vector \(\mathbf {a}\). It is noteworthy that this invariant \(I_5\) is not inversion safe (information about the sign of the smallest singular value \(\sigma _d\) is lost) due to the squaring implied by the right Cauchy-Green tensor \(\mathbf {C}= \mathbf {F}^T\mathbf {F}= \mathbf {S}^2 = \mathbf {V}\Sigma ^2\mathbf {V}^T\). Thus, inversion-safe anisotropic ARAP relies on the step function for capturing sign information that indicates element inversion from \(I_4\) to arrive at a perhaps more palatable expression of the PK1 (and thereby the Hessian), and with the effect of masking singularities that arise during anisotropic FE inversion or collapse.

2.3 Challenges with \(\Psi_AA\)-Energy Minimization

The constitutive law in Equation (9) is a highly versatile and practical model, but it unfortunately suffers from three key drawbacks that limit broader applicability.

The first drawback is that this energy is discontinuous at either \(I_4 = 0\) or \(I_5 = 0\), which are two new sources of singularity. This follows the fact that \(I_4 = 0\) does not generally imply \(I_5 = 0\) and vice versa since, for example, \(I_5\) is a weighted norm (i.e., involving a linear combination of the columns) of \(\mathbf {S}^2\) rather than \(\mathbf {S}\) as in the case of \(I_4\). Singularity thus shows itself once more due to a scaling by \(\delta (I_4 = 0) = \infty\) or a division-by-zero when \(I_5 = 0\) as shown in Equation (7). Patching the discontinuity with the missing value is proposed for overcoming this drawback by using a PK1 expression without the \(\delta\)-terms (see Equation (36) in the work of Kim et al. [2019]). However, solutions computed at \(I_4=0\) for minimization are then not congruent with the local energy landscape around this configuration due to the aforementioned discontinuity, leading to solver divergence.

The second drawback is that minimization of Equation (9) is numerically unstable due to a potential for division-by-zero as \(\sqrt {I_5} \rightarrow 0\), which will also severely impact an iterative solver’s rate of convergence since forces become intractably large.

The third drawback of Equation (9) is that it does not possess what we call the characteristic symmetry, which we detail in Section 5. Broadly, this property refers to the ability to reproduce isotropy from anisotropy by a summation of energies/terms that are parameterized by orthonormal vectors forming a basis of the same dimension as the FE.

Thus, Equation (9) is potentially unstable, as it still suffers from singularities that are fatal or will severely impact performance during minimization. The model also does not satisfy a key condition of symmetry, which we show to be vital for expanding the broader applicability of this energy.

3 Expressing I₄ in Proper Closed form

Following the preceding discussion, we design a new \(I_4\)-invariant expression with which any energy that is written purely in its terms will have a well-defined PK1 and Hessian in proper closed form for use in optimization with Newton-type methods. This section will specifically examine the role of rotation factor \(\mathbf {R}= \mathbf {S}^{-1}\mathbf {F}\) of the deformation gradient \(\mathbf {F}\). We defer discussions about our treatment of singularities to Section 4.

3.1 Volumetric FEs in 3D

Here we examine the rotation factor under the assumption of working with a volumetric FE that has a \({3\times 3}\) deformation gradient \(\mathbf {F}\) (e.g., a tetrahedron) to address the analytic drawbacks of \(I_4\) in this particular case.

The root cause of the issues posed by \(I_4\) is the rotation \(\mathbf {R}\) from polar factorization of \(\mathbf {F}\), where this presumably numerical factorization hinders a proper closed-form expression of the PK1 and Hessian of \(I_4\). We observe that such problems may be avoided by following Lin et al. [2022] to apply

\begin{equation} \mathbf {R}= \frac{\partial f}{\partial I_{C}}\frac{\partial I_{C}}{\partial \mathbf {F}} + \frac{\partial f}{\partial II_{C}}\frac{\partial II_{C}}{\partial \mathbf {F}} + \frac{\partial f}{\partial J}\frac{\partial J}{\partial \mathbf {F}} \end{equation}

(12)

as the general closed-form expression for this rotation factor, where

\begin{equation} I_{C} = {{\mathrm{tr}}(\mathbf {C})}=I_2 \qquad II_{C} = ||\mathbf {C}||^2_F \end{equation}

(13)

are the first two Cauchy-Green invariants [Bonet and Wood 2008]. The variable \(f\) refers to a root function of the quartic polynomial which is given by (cf. Equation (21) in the work of Lin et al. [2022])

\begin{equation} t^4-2I_Ct^2-8Jt{-I_C^2}+2II_C=0, \end{equation}

(14)

where the precise expression of this \(f\) that is of our particular interest is the one which evaluates to the trace

\begin{equation*} t := \mathrm{tr}(\mathbf {R}^T\mathbf {F}) = \sum _{i=1}^d \sigma _i= I_1. \end{equation*}

We also adopt the lower-order invariant \(J = \text{det}(\mathbf {F}) \equiv I_3\) here instead of the square root of the third Cauchy-Green invariant \(III_{C} = \mathrm{det}(\mathbf {C}) \equiv I_3^2\) because we are concerned with inversion-safe energies (i.e., our elements can flip). Otherwise, we can use \(\sqrt {III_C}\) in Equation (14) and Equation (12) instead of \(J\), like we do in Section 3.3 when discussing co-dimensional FEs.

Thus, revisiting Equation (2), we have

\begin{align} I_4 &= \mathbf {a}^T\mathbf {R}^T\mathbf {F}\mathbf {a}\nonumber \nonumber\\ &= \mathbf {a}^T\left(\underbrace{\frac{\partial f}{\partial I_C}\frac{\partial I_C}{\partial \mathbf {F}}^T\mathbf {F}+\frac{\partial f}{\partial II_C}\frac{\partial II_C}{\partial \mathbf {F}}^T\mathbf {F}+\frac{\partial f}{\partial J}\frac{\partial J}{\partial \mathbf {F}}^T\mathbf {F}}_{\mathbf {S}=\mathbf {R}^T\mathbf {F}}\right)\mathbf {a}, \end{align}

(15)

\begin{align} &= \mathbf {a}^T\left(\frac{\partial f}{\partial I_C}2\mathbf {F}^T\mathbf {F}+\frac{\partial f}{\partial II_C}4\mathbf {F}^T\mathbf {F}\mathbf {F}^T\mathbf {F}+\frac{\partial f}{\partial J}\text{cof}(\mathbf {F})^T\mathbf {F}\right)\mathbf {a}, \end{align}

(16)

\begin{align} &=2\frac{\partial f}{\partial I_C}IV_C+4\frac{\partial f}{\partial II_C}V_C+\frac{\partial f}{\partial J}{{I_8}} \end{align}

(17)

as our revised definition for \(I_4\) using the proper closed-form expression for the symmetric factor \(\mathbf {S}\) for a square \({3\times 3}\) deformation gradient \(\mathbf {F}\), where \(IV_C := I_5\) is the fourth Cauchy-Green invariant and \(V_C = \mathbf {a}^T\mathbf {C}^T\mathbf {C}\mathbf {a}\) is another even higher-order anisotropic invariant found in biomechanics literature [Weiss et al. 1996; Blemker et al. 2005; Chagnon et al. 2015].

The new low-order invariant \({{I_8}} = \mathbf {a}^T{J}\mathbf {I}\mathbf {a}= {J}\text{tr}(\mathbf {A})\) is a measure of scale and reflection of the vector \(\mathbf {a}\) due to the volume-ratio \(J\). The numbering has been shifted to accommodate pre-existing high- and low-order anisotropic invariants that follow below, which measure the change of \(\mathbf {a}\) w.r.t. stretch and shear via \(\mathbf {S}\). This new invariant is obtained using the relation \(\text{cof}(\mathbf {F})^T\mathbf {F}=\text{adj}(\mathbf {F})\mathbf {F}=\det (\mathbf {F})\mathbf {I}\), where \(\text{cof}(\mathbf {F})\) is the cofactor and \(\text{adj}(\mathbf {F}) = J\mathbf {F}^{-T}\) (i.e., the area map from Nanson’s rule [Bonet et al. 2015;, 2016; Poya et al. 2023]) is the adjugate of the deformation gradient.

PK1 and Hessian. We are now able to write out the proper closed-form expressions necessary to compute the exact forces \(\partial \Psi _4/\partial \mathbf {x}\) and their Jacobian \(\partial ^2 \Psi _4/\partial \mathbf {x}^2\) following Equations (3) and (6). The PK1 resolves to

\begin{align} \frac{\partial \Psi _4}{\partial \mathbf {F}}&=\frac{\partial \Psi _4}{\partial I_C} \frac{\partial I_C}{\partial \mathbf {F}} +\frac{\partial \Psi _4}{\partial II_C}\frac{\partial II_C}{\partial \mathbf {F}} +\frac{\partial \Psi _4}{\partial J}\frac{\partial J}{\partial \mathbf {F}}\nonumber \nonumber\\ &+{\frac{\partial \Psi _4}{\partial IV_C} \frac{\partial IV_C}{\partial \mathbf {F}} +\frac{\partial \Psi _4}{\partial V_C}\frac{\partial V_C}{\partial \mathbf {F}} +\frac{\partial \Psi _4}{\partial {I_8}}\frac{\partial {I_8}}{\partial \mathbf {F}}} \end{align}

(18)

with which we evaluate Equation (8) using

\begin{equation*} \mathcal {I} = \lbrace I_C, II_C, J, IV_C, V_C, {I_8}\rbrace \end{equation*}

as the set of invariants that are required to determine \(I_4\) when \(\mathbf {F}\) is a \({3\times 3}\) matrix. For brevity, we defer descriptions about the expressions for the partial derivatives \(\partial I_4/\partial i\) and \(\partial i/\partial \mathbf {F}\), \(i \in \mathcal {I}\) to Appendix A. Following Equation (6), we also can determine the Hessian by

\begin{equation} \frac{\partial ^2 \Psi _4}{\partial \mathbf {F}^2}=\sum _{i }\frac{\partial \Psi _4}{\partial i} \frac{\partial ^2 i}{\partial \mathbf {F}^2}+\sum _{i} \sum _{j} \frac{\partial ^2 \Psi _4}{\partial i\partial j}\frac{\partial j}{\partial \mathbf {F}} \otimes \frac{\partial i}{\partial \mathbf {F}}, \end{equation}

(19)

where our expressions for the partial derivatives are also described in Appendix A with \(i, j \in \mathcal {I}\).

Equations (18) and (19) thus crucially provide the missing PK1 and Hessian expressions (in proper closed form) that are necessary to analyze and simulate any energy written purely in terms of \(I_4\), when \(\mathbf {F}\) is a \({3\times 3}\) matrix.

3.2 Planar FEs in 2D

The proper closed-form expression of \(I_4\) for planar elements in 2D space is also Equation (15), but for a few simplifications given the reduced number of spatial variables that determine \(\mathbf {F}\) which is now a \(2 \times 2\) matrix. These simplifications build on the method of Lin et al. [2022]: we first rewrite the trace \(t\) in terms of the invariants as a root of the depressed quadratic polynomial \({t^2-I_C-2J=0}\). Then using implicit differentiation, we have

\begin{equation*} \frac{\partial f}{\partial I_C}=\frac{1}{2f}, \quad \text{and} \quad \frac{\partial f}{\partial J}=\frac{1}{f}, \end{equation*}

as the non-zero first-order partial derivatives of the root function \(f\) to give⁴

\begin{equation} \mathbf {R}= \frac{\partial f}{\partial \mathbf {F}} =\frac{\partial f}{\partial I_C}\frac{\partial I_C}{\partial \mathbf {F}} + \frac{\partial f}{\partial J}\frac{\partial J}{\partial \mathbf {F}}\in \rm I\!R^{2\times 2}, \end{equation}

(20)

as the corresponding rotation expression which now has two terms (rather than three as in Equation (12)) because the depressed quadratic polynomial is independent of \(II_C\). Equation (15) is thus reduced to

\begin{align} I_4 &= 2\frac{\partial f}{\partial I_C}IV_C+\frac{\partial f}{\partial J}{I_8} = \frac{1}{f}\left(IV_C+{I_8}\right), \end{align}

(21)

as the definition of \(I_4\) for planar elements in 2D space, where the PK1 \(\partial \Psi _4/\partial \mathbf {F}\) is evaluated as in Equation (7) and the Hessian \(\partial ^2 \Psi _4/\partial \mathbf {F}^2\) following Equation (19). These expressions are now dependent on the subset of invariants

\begin{equation*} \mathcal {I} = \lbrace {I_C, J,IV_C, {I_8}\rbrace }, \end{equation*}

and the corresponding the first- and second-order partial derivatives w.r.t. \(\mathbf {F}\) are given in Appendix A.

3.3 Co-dimensional Elements

For co-dimensional elements (planar elements in embedded in 3D space) with a \({3\times 2}\) deformation gradient \(\mathbf {F}\), the depressed quadratic resolves to \(t^2-I_C{-2\sqrt {III_C}}=0\) when rewriting the trace. The root function \(f\) then gives a rotation of the form

\begin{equation} \mathbf {R}= \frac{\partial f}{\partial \mathbf {F}} =\frac{\partial f}{\partial I_C}\frac{\partial I_C}{\partial \mathbf {F}} + \frac{\partial f}{\partial III_C}\frac{\partial III_C}{\partial \mathbf {F}}\in \rm I\!R^{3\times 2}, \end{equation}

(22)

which differs only slightly from Equation (20) because we now use \(\sqrt {III_C}\) instead of \(J\). Moreover, the stretch factor \(\mathbf {S}\) is guaranteed to be positive semi-definite, where the singular values \(\sigma _i \ge 0\) are always greater than or equal to zero because co-dimensional elements (like triangles in 3D space) cannot ever flip/invert. The corresponding definition of \(I_4\) will thus be

\begin{align} I_4&=\mathbf {a}^T\mathbf {S}\mathbf {a}=\mathbf {a}^T(\mathbf {R}^T\mathbf {F})\mathbf {a}, \end{align}

(23)

\begin{align} &= \mathbf {a}^T\left(\frac{\partial f}{\partial I_C}\frac{\partial I_C}{\partial \mathbf {F}}^T\mathbf {F}+ \frac{\partial f}{\partial III_C}\frac{\partial III_C}{\partial \mathbf {F}}^T\mathbf {F}\right) \mathbf {a}, \end{align}

(24)

\begin{align} &= \mathbf {a}^T\left(\frac{\partial f}{\partial I_C}2\mathbf {F}^T\mathbf {F}+ \frac{\partial f}{\partial III_C}2\det (\mathbf {F}^T\mathbf {F})\mathbf {I}\right) \mathbf {a}, \end{align}

(25)

\begin{align} &= \frac{\partial f}{\partial I_C}2IV_C+ \frac{\partial f}{\partial III_C}2{VI_{C}}, \end{align}

(26)

which is dependent on the subset of invariants

\begin{equation*} \mathcal {I} = \lbrace {I_C, III_C,IV_C, {VI_{C}}\rbrace }, \end{equation*}

where \({VI_{C}}={III_C\cdot }\text{tr}(\mathbf {A})\) is another new invariant measuring the scaling of the vector \(\mathbf {a}\) due to the squared area-ratio \(J^2 \equiv III_C\). Notation with Roman numerals is adopted here since this invariant is a function of the Cauchy-Green tensor \(\mathbf {C}\), which is common convention and more clearly distinguishes the invariant from our previous (and lower-order) definition \({I_8}={J}\text{tr}(\mathbf {A})\) in Section 3.1. We provide the PK1 and Hessian expressions of Equation (26) in Appendix A.

4 Fixing the Singularity Problem

As noted in Section 2.1, an energy \(\Psi _4\) that is written purely in terms of \(I_4\) will have singularity as an element is deformed into certain configurations, which is generally when \(\sigma _{d-1} = -\sigma _d\). Kim et al. [2019] address this problem by using a step function filter, but this introduces yet more challenges that affect stability and robustness of energies that are based on this filter. We describe in this section an approach to eliminating the singularities associated with \(I_4\), applying a mollifier (local higher-order smoothing function) or a regularizer on \(I_4\) in regions of the principal-stretch space that are near problematic configurations of a deformed element.

4.1 The Sources of Singularity

The sources of singularity lay hidden within the expression of \(I_4\) and are specifically found in the partial derivative terms \(\partial f/\partial i\) of the root function \(f\). These partial derivative terms share a common factor \(1/\alpha\) in their expressions (see also Equation (24) in the work of Lin et al. [2022]), where

\begin{align} \alpha &= 4f^3 - 4I_Cf - 8J, \end{align}

(27)

\begin{align} &= 8 (\sigma _1 + \sigma _2)(\sigma _1 + \sigma _3)(\sigma _2 + \sigma _3) \end{align}

(28)

is the precise sub-expression whose definition contains the sources of singularity when \(\mathbf {F}\) is a \(3\times 3\) matrix. Thus, we can see (once more, as in the work of Kim et al. [2019]) that singularity occurs as an element approaches \(\sigma _1 = -\sigma _2\), \(\sigma _1 = -\sigma _3\) and/or \(\sigma _2 = -\sigma _3\). Similarly, we have

\begin{align} \alpha &= f = \sigma _1 + \sigma _2 \end{align}

(29)

when \(\mathbf {F}\) is a \(2 \times 2\) matrix (cf. Section 3.2), where singularity occurs as we approach \(\sigma _1 = -\sigma _2\) (element collapsed to a point or fully inverted along one axis as a mirror image of its original shape).

For the particular case of co-dimensional elements, we require partial derivatives of \(f\) w.r.t. \(I_C\) and \({III_C}\) (not \(J\) for reasons discussed in Section 3.3) to give

\begin{equation} \alpha = f\sqrt {III_C} = (\sigma _1 + \sigma _2)\sqrt {\sigma _1\sigma _2}, \end{equation}

(30)

which indicates that \(I_4\) singularity occurs within such elements when they are collapsed to a point or a line (zero area). Equation (30) also illustrates why co-dimensional elements cannot ever flip/invert as having a negative singular value would imply a deformation that cannot be described without evoking the complex plane.⁵

A Note on the Derivatives. A further notable distinction between our analysis and that of Kim et al. [2019] is that we have a direct and analytic measure \(\alpha\) with which to determine when an element is approaching singularity, meaning that we effectively avoid the explicit burden of having to compute SVD. This implication is especially pertinent because optimizing an energy that is expressed in terms of \(\alpha\) (as we show in Section 4.2) will require the respective PK1 \(\partial \alpha /\partial \mathbf {F}\) and Hessian \(\partial ^2\alpha /\partial \mathbf {F}^2\) expressions, which would be rather cumbersome if quantities (e.g., the singular values \(\sigma _i\)) are expressed based on a numerically determined SVD. Conversely, our approach allows for expressions of this PK1 and Hessian in proper closed form, bypassing any need to compute SVD. Readers are referred to Appendix B for these expressions.

4.2 Two Solutions for Fixing Singularity

Leveraging the analytic expressions of the sources of singularity, we can now address this singularity without resorting to numerical differentiation of SVD. In this section, we will outline two robust methods designed specifically to tackle this issue.

Smooth Mollifier. Our first solution resolves singularity by multiplying the \(I_4\) energy by a mollifier, smoothly patching the region at \(\alpha {\le } \epsilon\) using an interpolating polynomial, where \(\epsilon = 10^{-3}\) is our threshold. A natural choice is the smoother step function [Perlin 2002]

\begin{equation} \mathcal {S}_2(x) = {\left\lbrace \begin{array}{ll} 0 & x \le 0\\ 6x^5 - 15x^4 + 10x^3 & 0\le x\le 1\\ 1& 1 \le x \end{array}\right.}, \end{equation}

(31)

which is at least twice differentiable for \(0\le x\le 1\), and from which we introduce the form

\begin{equation} (1 - t(\alpha))g(\mathbf {F}) + t(\alpha)I_4 \end{equation}

(32)

for filtering singularities from \(I_4\) where the blend parameter function \(t(\alpha) = \mathcal {S}_2 (\tfrac{\alpha - \epsilon }{2\epsilon })\) maps \(\alpha\) to the range \((0, 1)\) to give

\begin{equation} \hat{I}_4 = {\left\lbrace \begin{array}{ll} g(\mathbf {F}) & \alpha \lt \epsilon ,\\ (1 - t(\alpha))g(\mathbf {F}) + t(\alpha)I_4 & \epsilon \le \alpha \le 3\epsilon ,\\ I_4 & \alpha \gt 3\epsilon \end{array}\right.} \end{equation}

(33)

as our our final mollified expression of \(I_4\), and the function \(g(\mathbf {F})\) is a proxy measure temporarily standing in for \(I_4\) and with the requirement that its gradients do not tend toward infinity near or at \(\alpha = 0\). There are several choices for this proxy such as a constant \(c\), \(I_2\), \(I_3,\) or \(-\sqrt {I_5}\) from which we pick \(g(\mathbf {F}) = c\) with \(c=0\) as our constant.⁶ The choice is motivated by the simplicity of analyzing and computing the ensuing PK1 and Hessian expressions, and we found that using \(c\) (or \(-\sqrt {I_5}\)) is considerably faster than \(I_2\) and \(I_3\) (over 4\(\times\) less Newton solver iterations). The PK1 and Hessian of Equation (33) follow naturally, which are also then free of singularity.

Aberration Barrier. Our second solution augments an energy that may suffer from singularity with a correction term that serves to nudge collapsed FEs toward a more stable region with numerically tractable forces. Following Li et al. [2020] (see also Section 3.3 in the work of Smith et al. [2018] and Schüller et al. [2013] for related methods), we can add a regularized aberration barrier

\begin{align} b(\alpha) ={\left\lbrace \begin{array}{ll} -(\alpha -\epsilon)^2\ln \left(\frac{\alpha }{\epsilon }\right) & \alpha \le \epsilon ,\\ 0 & \alpha \gt \epsilon , \end{array}\right.} \end{align}

(34)

which inserts a peak of stiffness about the singularity region but uses \(\epsilon\) to smooth away the logarithmic singularity. Equation (34) is added to the local energy of an FE as an additional term, where we use \(\epsilon = 1\mathrm{e}{-2}\). An example is the form

\begin{equation*} \Psi _b(\mathbf {a},\alpha) = (I_4(\mathbf {a}) - 1)^2 + b(\alpha), \end{equation*}

which penalizes deformation along the prescribed direction of anisotropy \(\mathbf {a}\) in local material space, and toward the prescribed region of singularity (\(\alpha \le \epsilon\)).

While singularities can still occur when \(\alpha = 0\), the line search component of a Newton-type optimization scheme actively prevents \(\alpha\) from directly stepping into zero. It accomplishes this by finding a (local) minimum where \(\alpha\) is deterred from zero. To fortify our approach further, we implement an additional safeguard: if the algorithm is about to take a step that would cause \(\alpha\) to become less than another pre-defined threshold (e.g., \(1\mathrm{e}{-6}\)), we halve the step size. This strategy provides an additional layer of protection against the risk of \(\alpha\) nearing zero. In our testing, this method has proven to be robust, successfully passing all of the tests we conduct in Section 6.

5 Applications and Generalizations

To demonstrate the potential of the \(I_4\) expressions from the previous sections, we use them to rewrite three popular energies which have broad application in computer graphics and beyond: anisotropic ARAP [Kim et al. 2019; Alexa et al. 2000], the FE formulation of Baraff-Witkin cloth [Baraff and Witkin 1998; Kim 2020], and shape targeting [Klár et al. 2020a; Ichim et al. 2017]. In Appendix C and Appendix D, we supplement the generality of our approach by providing additional proofs that are extensions of the energies in this section.

5.1 Anisotropic ARAP: Revisited

In this section, we revisit the simplest possible form of the anisotropic version of the ARAP energy [Alexa et al. 2000]

\begin{align} \Psi _{\mathrm{ALX}} &= \frac{\mu }{2} \Vert (\mathbf {F}-\mathbf {R})\mathbf {a}\Vert _F^2 , \end{align}

(35)

\begin{align} &= \frac{\mu }{2} \left(\Vert \mathbf {F}\mathbf {a}\Vert _F^2 + \Vert \mathbf {R}\mathbf {a}\Vert _F^2 - 2\mathbf {a}^T\mathbf {F}^T\mathbf {R}\mathbf {a}\right), \end{align}

(36)

\begin{align} &= \frac{\mu }{2} \left(\mathbf {a}^T\mathbf {F}^T\mathbf {F}\mathbf {a}+ \mathbf {a}\mathbf {R}^T\mathbf {R}\mathbf {a}- 2\mathbf {a}^T\mathbf {F}^T\mathbf {R}\mathbf {a}\right), \end{align}

(37)

\begin{align} &= \frac{\mu }{2}\left(I_5 - 2{I}_4 + 1\right), \end{align}

(38)

which Kim et al. [2019] briefly discussed and showed to be easily expressible in terms of \(I_4\) and \(I_5\). The difference now is that we have an energy that is smooth and fully differentiable everywhere and guarantees robust simulation against failure by using Equation (34) or Equation (33) to address issues of singularity.

Isotropy from Anisotropy. Notably, the constitutive law suggested by Equation (38) is a general model of rotation-invariant anisotropic elasticity, which can be regarded as the foundational parametric description of isotropic (corotational) elasticity [Stomakhin et al. 2012; Lin et al. 2022; Chao et al. 2010; Liu et al. 2008; McAdams et al. 2011b; Sorkine and Alexa 2007; Finnendahl et al. 2023; Alexa et al. 2000]. To see this, let \(\lbrace {\mathbf {u}_i\rbrace }, 1\le i\le d\) be the set of orthonormal basis vectors (\(d = 2, 3\) is the FE dimension) spanning \(\rm I\!R^{d}\) with which we have

\begin{align} \sum _i^d \Psi _{\mathrm{ALX}}(\mathbf {u}_i) &= \sum _i^d I_{5}(\mathbf {u}_i) - 2\sum _i^d I_{4}(\mathbf {u}_i) + \sum _i^d 1, \end{align}

(39)

\begin{align} &= I_2 - 2I_1 + d, \end{align}

(40)

\begin{align} &= \Psi _{i{\rm\small ARAP}} , \end{align}

(41)

as the anisotropic re-presentation of isotropic ARAP energy that is written in terms of the low-order invariants [Smith et al. 2019] and evaluated with analytic rotations [Lin et al. 2022]. This follows from

\begin{equation} \sum _i^d \mathbf {u}_i^T \left(\prod _j^m \mathbf {S}\right) \mathbf {u}_i = \mathrm{tr}(\mathbf {S}^m), \end{equation}

(42)

which states that the sum of \(d\) anisotropic invariants that are of order \(m\) and parameterized by orthonormal vectors \(\mathbf {u}_i\) is an isotropic invariant of the same order.

Proof.

Let \(\mathbf {U}= \begin{bmatrix} \mathbf {u}_{1} & \ldots & \mathbf {u}_{d} \\ \end{bmatrix}\in \rm I\!R^{d\times {d}}\) be an orthogonal matrix with \(\mathbf {U}^T\mathbf {U}= \mathbf {U}\mathbf {U}^T = \mathbf {I}\). Then for \(m=1\), we can write Equation (42) as

\begin{equation} \sum _i^d \mathbf {u}_i^T\mathbf {S}\mathbf {u}_i = \mathrm{tr}(\mathbf {U}^T\mathbf {S}\mathbf {U}) = \mathrm{tr}({\mathbf {U}\mathbf {U}^T}\mathbf {S}) = \mathrm{tr}(\mathbf {S}) = I_1, \end{equation}

(43)

based on the cyclic property of the trace. Since Equation (43) holds for any symmetric matrix \(\mathbf {S}\), we thus also have

\begin{equation} \sum _i^d \mathbf {u}_i^T\mathbf {S}^T\mathbf {S}\mathbf {u}_i = \mathrm{tr}(\mathbf {S}^2) = I_2, \end{equation}

(44)

from which we get the first term in Equation (40). □

It is by generalizing to arbitrarily high-order \(m\) that we obtain Equation (42) as the power sum of the roots (i.e., singular values) of the characteristic polynomial of the deformation gradient \(\mathbf {F}\). This power sum is a symmetric polynomial of order \(m\) following the well-known Newton-Girard identities [Mead 1992] which imply that any isotropic invariant of degree \(m\) may also be written in terms of lower-order invariants (or elementary symmetric polynomials in \(\sigma _i\)) for exploring the space of invariant-expressible energies that are also inversion safe.

Characteristic Symmetry. We designate Equation (35) as having a property we call characteristic symmetry given that each term of Equation (39) (via Equation (42)) resolves to a symmetric polynomial involving the roots of the characteristic polynomial of \(\mathbf {F}\). This property is not preserved with energies like \(\Psi _{AA}\) (Equation (9)) or the Baraff-Witkin stretching model (Section 5.2), which fail to reproduce ARAP. The generality of Equation (35) also means that our model is readily extensible for catering to the various forms of anisotropy (cf. Figure 1). Layering \(d\) energies (e.g., \(d=3\)) as in Equation (39) achieves isotropy, where each energy is associated with the same material parameters. By prescribing a different material parameter for each layered energy, we obtain orthotropy (cf. Figure 3). Finally, a generally anisotropic material is modeled by layering an arbitrary number of anisotropic energies, each prescribed with a unique material parameter and direction.

Fig. 3.

5.2 Cloth

We will now use the preceding discussion to show that a low-order Baraff-Witkin cloth model [Baraff and Witkin 1998; Kim 2020] can be derived using \(I_4\) to yield an invariant-expressible FE formulation (cf. Equation (38)) that is simpler and even faster to optimize. To demonstrate this, we first rename and renumber the invariants found in the work of Kim [2020] to maintain consistency with the equations in preceding sections, \(I_7(\mathbf {a}, \mathbf {b}) = \mathbf {a}^T\mathbf {F}^T\mathbf {F}\mathbf {b}\). The numbering has been shifted so that we can later introduce a new, lower-order invariant, \(I_6\). These are the cross-fiber invariants, penalizing shear with (normalized) fiber directions \(\mathbf {a},\mathbf {b}\in \rm I\!R^{2\times 1}\) to give

\begin{align} \Psi _{BW} &= \frac{1}{2}\left(\left(\sqrt {I_5(\mathbf {a})} - 1\right)^2 + \left(\sqrt {I_5(\mathbf {b})} - 1\right)^2\right) + I_7(\mathbf {a}, \mathbf {b})^2, \end{align}

(45)

as the FE formulation of the Baraff-Witkin cloth model in the work of Kim [2020]. These fiber directions satisfy the orthogonality condition \(\mathbf {a}^T\mathbf {b}= 0\), which if relaxed leads to \((I_7(\mathbf {a}, \mathbf {b}) - \mathbf {a}^T\mathbf {b})^2\) as the generalized shear corresponding to the third term in Equation (45). Notably, the first two (i.e., stretching) terms of Equation (45)

\begin{align} \underbrace{{I_5(\mathbf {a})} + {I_5(\mathbf {b})}}_{I_2} + 2\underbrace{\left(\sqrt {I_5(\mathbf {a})} + \sqrt {I_5(\mathbf {b})}\right)}_{\approx I_1} + 2 \approx {\rm\small ARAP} \end{align}

(46)

form an anisotropic approximation of isotropic ARAP (Equation (40)), which unfortunately indicates that is not possible to express Equation (45) as an isotropic FEM energy.

Lower-Order Form. Using our analysis in Equation (39), we also have

\begin{align} \Psi _{BW46} &= \Vert (\mathbf {F}-\mathbf {R})\mathbf {a}\Vert _F^2 + \Vert (\mathbf {F}-\mathbf {R})\mathbf {b}\Vert _F^2, \end{align}

(47)

\begin{align} &= I_5(\mathbf {a})+I_5(\mathbf {b}) - 2({I}_4(\mathbf {a}) + {I}_4(\mathbf {b})) + 2, \end{align}

(48)

\begin{align} &= I_2 - 2(I_1 - 1), \end{align}

(49)

as the lower-order form of Equation (45), where the PK1 and Hessian follow naturally. The fact that we can reduce anisotropic cloth energy to invariant-expressible ARAP is significant,⁷ but the use of Equation (47) has yet one problem that we must address.

Decoupling Stretch and Shear. The energy \(\Psi _{BW46}\) will simultaneously penalize the stretch and shear components of cloth deformation w.r.t. the fiber directions, which removes our ability to have decoupled control over these components like Equation (45). To resolve this, we introduce

\begin{align} \tilde{\Psi }_{BW461} &= \underbrace{I_4(\mathbf {a})^2+I_4(\mathbf {b})^2+2I_6(\mathbf {a},\mathbf {b})^2}_{I_2}-2 (I_1- 1), \end{align}

(50)

as the equivalent form (proof in Appendix D), where \(I_6(\mathbf {a}, \mathbf {b}) = \mathbf {a}^T\mathbf {F}^T\mathbf {R}\mathbf {b}= \mathbf {a}^T\mathbf {S}\mathbf {b}\) is our novel cross-fiber invariant. The revised energy is then

\begin{align} \Psi _{BW461} &= \frac{\mu _1}{2} \left(I_4(\mathbf {a})^2+I_4(\mathbf {b})^2-2 (I_1 - 1)\right) + \mu _2 I_6(\mathbf {a},\mathbf {b})^2, \end{align}

(51)

\begin{align} &\equiv \frac{\mu _1}{2} \left((I_4(\mathbf {a})-1)^2+(I_4(\mathbf {b})-1)^2\right) + \mu _2 I_6(\mathbf {a},\mathbf {b})^2, \end{align}

(52)

which is the decoupled expression of our low-order cloth model that is derived from isotropic FEM energy in Equation (47). We also have \(\mu _1\) and \(\mu _2\) as the parameters determining the amount by which stretch and shear are penalized, respectively.

Relaxing the orthogonality condition also gives

\begin{align} \Psi _{BW46\ast } &= \frac{\mu _1}{2} \left(I_4(\mathbf {a})^2+I_4(\mathbf {b})^2- 2((\underbrace{{I}_4(\mathbf {a}) +{I}_4(\mathbf {b})}_{\approx I_1}) - 1)\right) \nonumber \nonumber\\ &+ \mu _2 (I_6(\mathbf {a}, \mathbf {b}) - \mathbf {a}^T\mathbf {b})^2, \end{align}

(53)

as the revised decoupled form with generalized shear. The form of Equation (53) follows the fact that Equation (51) is presented under the assumption that the fiber directions \(\mathbf {a}\) and \(\mathbf {b}\) are orthonormal, which allows for some simplification and substitutions that are not possible with Equation (53).

5.3 Shape Targeting

In this section, we revisit the problem of simulating so-called ‘active’ elasticity, where a target shape is given for each FE. The formulation is a natural extension to the fiber-based anisotropic models that we have described in preceding sections, where a multi-dimensional analogue of the anisotropy vector is now used.

Constitutive Model. Klár et al. [2020a] showed that the ‘Phace’ facial modeling energy originating from Ichim et al. [2017] can be recast as an elastic constitutive model (see Equation (1) in the work of Klár et al. [2020b]) by

\begin{equation} \Psi _{ST} = \frac{\mu }{2}\Vert \mathbf {F}- \mathbf {R}\mathbf {S}_t\Vert _F^2 + \frac{\lambda }{2}(J - \det (\mathbf {S}_t))^2, \end{equation}

(54)

where \(\mathbf {S}_t\) is a symmetric 3 \(\times\) 3 matrix that is used to encode a multi-dimensional ‘activation’ in the same fashion as the 1D anisotropic ARAP model described in previous sections. The variables \(\mu\) and \(\lambda\) are material parameters. We examine Equation (54) to show that this too can be further cast into a fully differentiable low-order invariant-expressible energy, revealing a more direct relationship among shape targeting [Ichim et al. 2017; Klár et al. 2020a], anisotropic ARAP energy in Equation (38), and \(I_4\).

Rewriting the Constitutive Model. To show this, we focus on the first term of Equation (54), which is most relevant,⁸ and drop the material parameter \(\mu\) for brevity to give

\begin{align} \Psi _{ST}&=\Vert \mathbf {F}-\mathbf {R}\mathbf {S}_t\Vert _F^2 , \end{align}

(55)

\begin{align} &=\Vert \mathbf {F}\Vert _F^2+\Vert \mathbf {S}_t\Vert _F^2-2\text{tr}(\mathbf {F}^T\mathbf {R}\mathbf {S}_t), \end{align}

(56)

\begin{align} &=\Vert \mathbf {F}\Vert _F^2+\Vert \mathbf {S}_t\Vert _F^2-2\text{tr}(\mathbf {S}\mathbf {S}_t), \end{align}

(57)

\begin{align} &=\Vert \mathbf {F}\Vert _F^2+\Vert \mathbf {S}_t\Vert _F^2-2\text{tr}(\mathbf {S}\mathbf {S}_t)+{d}-{d}+2\text{tr}(\mathbf {S})-2\text{tr}(\mathbf {S}), \end{align}

(58)

\begin{align} &=\Vert \mathbf {F}\Vert _F^2-2\text{tr}(\mathbf {S})+{d}+\Vert \mathbf {S}_t\Vert _F^2-2\text{tr}(\mathbf {S}\mathbf {S}_t)-{d}+2\text{tr}(\mathbf {S}), \end{align}

(59)

\begin{align} &=\Psi _{i{\rm\small ARAP}}+\Vert \mathbf {S}_t\Vert _F^2-2\text{tr}(\mathbf {S}\mathbf {S}_t)-{d}+2\text{tr}(\mathbf {S}), \end{align}

(60)

\begin{align} &=\Psi _{i{\rm\small ARAP}}+2\text{tr}(\mathbf {S}(\mathbf {I}-\mathbf {S}_t))+\Vert \mathbf {S}_t\Vert _F^2-{d}, \end{align}

(61)

as the algebraic expansion, where implicit ARAP [Lin et al. 2022] appears immediately as the isotropic component. We can further reduce Equation (61) by focusing on the second term to give

\begin{align} 2\text{tr}(\mathbf {S}(\mathbf {I}-\mathbf {S}_t))&=2\text{tr} \left(\mathbf {S}\sum _{i=1}^{d}(1-\lambda _i)\mathbf {v}_i\mathbf {v}_i^T\right), \end{align}

(62)

\begin{align} &=\sum _{i=1}^{d}2(1-\lambda _i)\text{tr}(\mathbf {S}\mathbf {v}_i\mathbf {v}_i^T), \end{align}

(63)

\begin{align} &=\sum _{i=1}^32(1-\lambda _i)I_4(\mathbf {v}_i), \end{align}

(64)

since \(\mathbf {S}_t = \sum _{i=1}^{d}\lambda _i\mathbf {v}_i\mathbf {v}_i^T\) is a symmetric matrix with real eigenvalues \(\lambda _i\) and eigenvectors \(\mathbf {v}_i\) from which we get \(\mathbf {I}-\mathbf {S}_t=\sum _{i=1}^{d}(1-\lambda _i)\mathbf {v}_i\mathbf {v}_i^T\). Thus, the final expression of the stretch term in elastic energy associated with an ‘active’ material [Klár et al. 2020a] is

\begin{align} \Psi _{ST}&=\Psi _{i{\rm\small ARAP}}+\sum _{i=1}^{d}2(1-\lambda _i)I_4(\mathbf {v}_i)+\Vert \mathbf {S}_t\Vert _F^2-{d,} \end{align}

(65)

which is the invariant-based expression of Equation (55) where the last two terms are constants since they are independent of the deformation gradient \(\mathbf {F}(\mathbf {x})\).

A comparison of Equations (65) and (55) shows that they achieve the same minimum (rest shape) for the same values of \(\mathbf {F}\). Equation (65), however, offers a more direct interface for specifying the shape target \(\mathbf {S}_t\) via eigenpairs \((\lambda _i, \mathbf {v}_i)\) as a way to induce anisotropic stiffening/softening in a specific direction unlike the algebraic relation of the constitutive model presented in the work of Klár et al. [2020b]. We can likewise apply the solutions in Section 4 to address any potential issues relating to force singularity since Equation (65) contains \(I_4\). Another advantage of our formulation is that simply setting \(\lambda _i = 0\), \(\forall i \gt 1\) will recover the canonical form used in muscle-based models, where it is assumed that the isotropic and anisotropic components of the model \(\Psi = \Psi _{iso} + \Psi _{aniso}\) are additively layered. The added benefit is that there is a clean separation of the qualitative behavior of each term to allow for isolated analysis and control.

6 Results

In this section, we present discussions about our analysis and combine with experiments to characterize the behavior of our energies and performance. We optimized our energies using a projected Newton solver with line search. The Preconditioned Conjugate Gradient (PCG) method is used to solve the linear system arising in each Newton iteration, which we implement with the Eigen library. In our testing for the convergence of the Newton solver, we employed a dynamic simulation with a value of \(10^{-2}L\Delta t\) and a static simulation with a value of \(10^{-4}L\), where \(L\) represents the length of the diagonal of the bounding box of the object being simulated, and \(\Delta t\) is the size of the timestep. Additionally, we set \(\delta =10^{-6}\) as our PCG threshold. A standard diagonal preconditioner is used for PCG. Our results are produced on an Ubuntu system with a 20-core 3.8-GHZ Intel Core I9 CPU and 32 GB of RAM. A summary of our timings is also provided in Table 1 for reference.

Table 1.

Scenario	Elem	#Elems	\(E\)	\(\nu\)	\(\Delta {t}\)	#\(\Delta {t}\)	Newt’ iters	PCG Iters	Time
Figure 3 (compress)	tet	54,985	1e3	0.3	1e-2	50	1.9	906.3	0.41
Figure 3 (pull)	tet	34,992	1e3	0.3	1e-2	50	1.5	874.7	0.38
Figure 4	tet	34,992	1e5	0.0	1e-2	2,000	2.8	411.3	0.71
Figure 5	tet	637,584	1e4	0.0	1e-2	200	26.0	15,663.6	128.76
Figure 6 (compress)	tet	17,840	1e4	0.0	1e-2	100	1.3	490.6	0.16
	tet	17,840	1e4	0.3	1e-2	100	1.4	644.8	0.18
	tet	17,840	1e4	0.49	1e-2	100	1.8	2,509.4	0.31
Figure 6 (pull)	tet	17,840	1e4	0.0	1e-2	200	1.6	544.7	0.18
	tet	17,840	1e4	0.3	1e-2	200	1.8	775.8	0.22
	tet	17,840	1e4	0.49	1e-2	200	3.1	3,685.7	0.44
Figure 9(a)	tet	308,200	1e3	0.0	1e-2	400	3.5	1,738.9	8.54
Figure 9(b)	tet	308,200	1e3	0.0	1e-2	400	3.8	24,010.2	53.02
Figure 9(c)	tet	308,200	1e5	0.0	1e-2	400	1.0	4,980.3	12.61
Figure 13 (interp)	hex	200,739	1e2	0.4	–	–	3.3	77,701.5	1,026.22
Figure 13 (interact)	hex	200,739	1e2	0.4	1e-2	200	1.0	745.9	68.04
Figure 15(b)	tet	292,476	1e2	0.0	–	–	25	85,179	219.69
Figure 15(c)	tet	292,476	1e2	0.0	–	–	67	269,221	685.93

Table 1. Performance Data of Dynamic Volumetric Simulation

Measurements are given as averages over the total number of timesteps. \(E\) refers to Young’s modulus, and \(\nu\) is Poisson’s ratio. \(\Delta {t}\) is the timestep size, and #\(\Delta {t}\) is the number of timesteps. All time measurements are in units of seconds.

Inversion-Safe ARAP. Figure 4 shows the result of reproducing the standard ‘scramble’ benchmark test that is common in FEM literature (see, e.g., the work of Lin et al. [2022], Sin et al. [2013], Smith et al. [2018], and Stomakhin et al. [2012]). The benchmark is simulated with an isotropic energy that is constructed using Equation (39) with \(d=3\) orthogonal anisotropy directions \(\mathbf {u}_i\). We randomly position the vertices of the Armadillo within a cube whose volume is twice that of the Armadillo’s rest volume, except a small subset of vertices on the hands and feet which remain in their original positions. This benchmark demonstrates the overall robustness of our method, showing our ability to withstand inversion and flattening for recovering the original rest configuration gracefully. In contrast, using the anisotropic ARAP energy proposed by Kim et al. [2019] (cf. Equation (9)) in the same manner fails at the very first step, which highlights the significance of our analysis. The result of Figure 4 also provides further evidence for the validity and practical utility of our finding in Equation (42).

Fig. 4.

Figure 5 provides a visual demonstration of another practically challenging scenario, where the Stanford Bunny is inverted such that all of its elements are in a singularity-inducing state. We use this experiment (same energy as in Figure 4) as a way to to provide further insight into the practical utility, robustness, and potential applications of our model under unique and challenging conditions. Our method successfully recovers the original bunny shape when using either of our solutions that we propose to address singularity.

Fig. 5.

We further test our anisotropic rewriting of the isotropic ARAP energy (Equation (39)) with a pulling and compression test, where we also include volume preservation by adding the term \(\tfrac{\lambda }{2}(J-1)^2\) to our energy as in previous works (e.g., those of Lin et al. [2022], Smith et al. [2018], and Stomakhin et al. [2012]). We simulate a cylinder that is pulled in one instance and compressed in another. We vary Poisson’s ratio \(\nu =0,0.3,0.49\) across three distinct values while maintaining Young’s modulus \(E=1\mathrm{e}4\) at a constant. Our results are shown in Figure 6, which demonstrate the range of material behaviors that are achieved with our method while also penalizing volume loss.

Fig. 6.

Singularity-Fixing solutions. To analyze the \(I_4\)-energy profile near singularity, we plot several graphs in Figures 7 and 8 which characterize the original (unfiltered) behavior of \(I_4\) as well as our singularity-fixing solutions discussed in Section 4. The unfiltered energy (cf. Equation (21)) can be observed to clearly shoot off to infinity as \(\alpha \rightarrow 0\). Our first solution (Equation (33)) will ensure that the local \(I_4\)-force magnitudes are smoothly attenuated to zero near singularity by using the mollifier. Figure 8 shows our aberration barrier profile corresponding to our second solution (Equation (34)), which guarantees, by construction, that the material exhibits a very strong reaction to extreme deformations in the direction of singularity. This is a property unique to this solution; using Equation (33) will allow the material to deform to the singularity state \(\alpha \approx 0\) while absorbing a finite amount of energy. In contrast, Equation (34) prevents the material from deforming toward the singularity state by injecting a finite amount of energy. In contrast to the work of Kim et al. [2019], our solutions ensure that singularities are eliminated in practically all configurations of an FE by guaranteeing energy smoothness everywhere and thus robustness (with additional safeguards) against any potential for failure.

Fig. 7.

Fig. 8.

Anisotropic ARAP Revisited. As another validation of our formulation of \(I_4\), we simulate an octopus model resembling the experiment in the work of Kim et al. [2019] to show our ability to likewise model anisotropic stiffening effects. We model the main body of the octopus using isotropic material \(\Psi _{i{\rm\small ARAP}}\) (see Figure 9(a)). We also define fibers along each of the tentacles for modeling localized control of anisotropic stiffness with \(\Psi _{i{\rm\small ARAP}} + \tfrac{\gamma }{2}(I_4(\mathbf {a}) - 1)^2\), where \(\gamma\) is a parameter controlling this stiffness. As shown in Figure 9(b), stiffening of the tentacles leads to a natural increase in rigidity or slackness along their span, while the bending behavior in the other directions remains the same. Conversely, increasing the general the stiffness of the isotropic model representing the octopus body will have a global effect on the model (see Figure 9(c)), unlike the localized changes observed when we stiffen only the tentacle fibers. The benefit of our analysis is that we now have the \(I_4\) Hessian expression, which allows us to directly employ this invariant for simulation.

Fig. 9.

Cloth Simulation. To evaluate our model for cloth simulation with decoupled control over stretch and shear, we set up two scenarios for analyzing stretching and draped motion which are shown in Figures 10 through 12. We compare our \(I_4\)-based model in Equation (51) with the FEM formulation of Baraff-Witkin cloth by Kim [2020] (cf. Equation (45)), which is well known and widely recognized in the field. In both scenarios, our model is able to produce comparable deformation behavior, and with faster performance given the lower-order nature of our energy. Table 2 also shows that our method yields a lower number of Newton and PCG iterations most of the time, which is significant yielding up to 1.28 to \(2.91\times\) speedup over the model by Kim [2020]. We use \(\mathbf {a}= [1, 1]\) and \(\mathbf {b}= [1, -1]\) in Figures 10 and 12 (normalized) as our anisotropy vectors since using the canonical basis vectors reduces the amount of wrinkling observed when compared to the case of any other pair of orthonormal bases. Overall, our model provides a more efficient and general solution for cloth simulation, with comparable quality of results to the state of the art [Kim 2020]. We have compared our results against those of Kim [2020] by using implementations with and without analytic eigensystems for completeness, as we found the approach using explicit/numerical projections of the full element Hessian (prior to global assembly) to be generally faster than analytically projecting the eigensystems of individual energy terms as described in the work of Kim [2020]. This difference in performance may be due to the modifications of the local energy landscape (curvature) caused by the term-by-term analytic projections versus full element Hessian projection. Moreover, the sign and magnitude of each eigenvalue of the (full) Hessian may change dramatically as multiple Hessian terms are incrementally added together, where the severity of this change is dependent on the number of terms.

Table 2.

Scenario	Energy model	Stretch	Shear	Bending	#\(\Delta {t}\)	Newton iters	PCG iters	Time (sec)
Isotropic drop	Equation (49) [Lin et al. 2022]	40	40	1e-3	1,000	1.3	1,270.2	0.67
	Equation (51) Ours	40	80	1e-3	1,000	1.4	1,540.0	0.81
	Equation (45) [Kim 2020]	40	80	1e-3	1,000	1.4	2,830.0	0.95
	Equation (45) [Kim 2020] (analytic)	40	80	1e-3	1,000	1.5	2,876.2	1.00
Anisotropic drop	Ours	40	80	1e-3	1,000	2.5	3,039.3	1.31
Anisotropic drop	Kim [2020]	40	80	1e-3	1,000	2.6	5,396.4	1.73
Isotropic stretch	Equation (49) [Lin et al. 2022]	20	20	1e-3	200	1.1	587.3	0.56
	Equation (51) Ours	20	30	1e-3	200	1.6	946.3	0.65
	Equation (45) [Kim 2020]	20	30	1e-3	200	3.0	3,504.5	1.44
	Equation (45) [Kim 2020] (analytic)	20	30	1e-3	200	1.9	2,422.7	1.00
Anisotropic stretch	Ours	20	20	0	200	2.31	1,451.3	0.37
	Ours	20	80	0	200	2.34	1,532.4	0.42
	[Kim 2020]	20	20	0	200	5.3	5,759.5	1.08
	[Kim 2020]	20	80	0	200	3.7	5,451.1	1.10

Table 2. Performance Data of Dynamic Cloth Simulation with 32,768 Triangle Elements

The timestep size is set to 1e-2 seconds. Measurements are given as an average over the total number of timesteps.

Fig. 10.

Fig. 11.

Fig. 12.

Active Elasticity Using Shape Targets. Building on a template hexahedral mesh, we conduct two demonstrations for evaluating our shape targeting energy in Equation (65). The first example optimizes energy over a sequence of target matrices per element. Here we effectively optimize the deformation process by driving the model to achieve specific given shape transformations over time. In the second example, we interact with a model containing pre-specified shape targets, forcefully manipulating the template by pushing on both cheeks of the retargeted face. Our results are shown in Figure 13, which are implemented following the hexahedral FE formulation described in the work of Irving et al. [2006]. These examples showcase the versatility of our approach in dynamically controlling and optimizing mesh deformations by using shape targets.

Fig. 13.

Parameterization. Figure 14 shows a different application of our analysis, this time for accomplishing the task of parameterizing (i.e., fitting a texture to) a mesh, where we have used the 2D version of Equation (39) with two anisotropy directions and \(\mathbf {F}\in \rm I\!R^{2\times 2}\). The parameterization is initialized with Tutte embedding. In Figure 14(a), we show the result of standard (isotropic) parameterization, where the anisotropy directions correspond to the canonical basis vectors with \(\mathbf {u}_1 = [1, 0]\) and \(\mathbf {u}_2 = [0, 1]\). In Figure 14(b), we show the result we obtain after multiplying all \(\mathbf {u}_1\)-related terms in the energy by a factor of \(10\times\). We further rotate both vectors \(\mathbf {u}_i\) by 45 degrees and show the result in Figure 14(c) as a demonstration of the generality of our fully differentiable and robust expressions of \(I_4\).

Fig. 14.

Generating Caricatured Models from Curvature. We also apply our method to generate exaggerated artistic versions of 3D models by amplifying the characteristic traits based on curvature. Our results are shown in Figure 15. For a given model, we first compute the normal and curvature at each surface vertex. These attributes are then propagated to the interior vertices of the respective tetrahedral mesh by solving the Laplace equation with the surface attributes as the boundary condition such that we can calculate the average curvature value and normal vector for each tetrahedron using information stored at the vertices. We then minimize the energy

\begin{align} \Psi _{\text{C}} = \Psi _{i{\rm\small ARAP}} + \beta cI_4(\mathbf {n}), \end{align}

(66)

where \(\beta\) represents the degree to which curvature characteristics should be emphasized. If \(\beta\) is positive, regions with positive curvature expand, whereas those with negative curvature contract, and vice versa. The variable \(c\) denotes the ratio of element curvature to the object’s average curvature, and \(\mathbf {n}\) is the normal direction. Equation (66) is in fact just one form of Equation (65), where \(2(1-\lambda _1)={{\beta }}c\), \(\mathbf {v}_1=\mathbf {n}\), and \(\lambda _{2,3} = 1\) with \(\mathbf {v}_{2,3}\) representing any unit vectors forming an orthonormal basis with \(\mathbf {n}\).

Fig. 15.

7 Conclusion and Discussion

We have presented a new formulation for simulating anisotropic material that also generalizes to isotropic elasticity in low-order invariant-expressible form, and with robust behavior near singularity-inducing states. This method works by using analytic rotations rather than numerically determined energy terms (e.g., from numerical SVD), which allow us to derive novel analytic and thus fully differentiable expressions for the low-order anisotropic invariants \(I_4\) and \(I_6\). A further key result of our analysis is a novel perspective from which to identify the sources of force singularity to then be able to construct numerous inversion-safe energies that are expressed in terms of these invariants. Crucially, our solutions to addressing issues of singularity eliminate the perilous impact of problematic geometric configurations on the energy and its gradients. With these methods in place, we can construct an overall timestepping energy for low-order anisotropic elastic deformation that can leverage superlinear convergence and robustness of Newton-type stepping.

7.1 Discussion

Real-world materials are generally classified as either isotropic or anisotropic with further distinctions between the three prevalent types of anisotropy, namely transversely isotropic, orthotropic, and generally anisotropic materials. We have demonstrated a low-order formulation of transverse isotropy that is a superset of the other forms, providing a foundation upon which numerous models of deformation can be constructed. This is significant, as it applies to the innumerable real-world materials which exhibit unique (and sometimes uniform) properties when measured along different axes. Our results thus also suggest a potential to apply our analysis for solving problems in areas beyond our most immediate use cases like physics-driven actuated soft bodies [Yang et al. 2022; Srinivasan et al. 2021; Klár et al. 2020a; Lee et al. 2018; Du et al. 2021], biomechanics [Teran et al. 2003; Blemker and Delp 2005; Chagnon et al. 2015; Sifakis et al. 2005], engineering [Bonet and Burton 1998; Weiss et al. 1996], and image processing [Welk et al. 2006], among others.

Overall, our singularity-fixing solutions extend \(I_4\) and \(I_6\) over the critical region spanning the problematic collapsed and inverted regimes to obtain valid forces when \(\alpha \approx 0\). Our modifications can then be applied to arbitrary energies that can be written in their terms as we show in Section 5. The resulting expressions are identical to the physical model most of the time, and they also allow the simulation to continue even if some elements are within the prescribed singularity region. Our extensions also provide at least \(C^1\) continuity around the problematic collapsed and inverted cases, which avoids sudden jumps or oscillations that could effect neighboring elements.

In practice, we apply these singularity-fixing solutions to energies that explicitly contain \(I_4\) and/or \(I_6\) like Equations (48), (53), and (66). This includes instances where the respective energy is equivalent (i.e., can be reduced) to isotropic ARAP⁹ like Equation (48) because individual force terms (\(\partial I_4/\partial \mathbf {x}\) and \(\partial I_6/\partial \mathbf {x}\)) can become singular as \(\alpha \rightarrow 0\).

By using \(g(\mathbf {F}) = c\) with \(c=0\) in our first solution to address singularity, the \(I_4\) and \(I_6\) terms of a local energy will be effectively weakened as an element enters the prescribed region \(\alpha \le \varepsilon\), by attenuating these terms near \(\alpha =0\) to remove singularity with a mollifier at a given deformation limit \(\varepsilon\). Such term-specific weakening is inspired by the behavior of St Venant-Kirchhoff elasticity under uni-axial compression of an element until total collapse [Picinbono et al. 2000; Sifakis and Barbic 2012]. In our particular case, once the critical threshold is reached (\(\alpha = \varepsilon\)), the strength of restorative \(I_4\) force reaches a maximum. Further deformation in the direction of singularity will be met with decreasing resistance—in fact, the restorative force, for example,

\begin{equation*} \lim _{\alpha \rightarrow 0} \frac{\partial \hat{I}_4}{\partial \mathbf {x}} = 0, \end{equation*}

will vanish in the limit. The difference with St Venant-Kirchoff, however, is that our imposed weakening is strictly local as it applies only to \(I_4\) (and \(I_6\)) terms of an energy, and our proxy function \(g(\mathbf {F})\) is active only within elements that are currently configured to be in the small region \(\alpha \le \varepsilon\). For situations where the local energy is written purely in terms of \(I_4\) and/or \(I_6\), neighboring elements can serve to act as external tie-breakers in the limit, providing the load necessary to restore an element if the prescribed local material is rendered too weak.

It follows that alternative choices for \(g(\mathbf {F})\) will of course lead to some other (perhaps more complex) effects than mere weakening, but we reserve proper exploration of such avenues for future work. This is especially relevant for aforementioned energies that are expressed purely in terms of \(I_4\) and/or \(I_6\) where, for example, several (perhaps all) adjacent volumetric or co-dimensional elements are collapsed to a line. Such cases may require an alternative choice of \(g(\mathbf {F})\) to serve as tie-breaker by providing the force necessary to restore local material.

Our alternative singularity-fixing solution in Equation (34) will introduce a function exhibiting rapid growth as \(\alpha \rightarrow 0\) to modify the local constitutive model so as to remove singularity in the style of classic Neohookean elasticity with its near origin barrier [Bonet and Wood 2008]. The difference is that we have a ‘near singularity’ barrier instead. The fact that this solution

defines a generally impassable barrier as \(\alpha \rightarrow 0\) implies that there is technically no mechanism for handling what happens when, accidentally, a simulated element is forced into a (theoretically impossible) singular configuration. In such cases, our aberration barrier energy and stress are undefined, since \(\alpha = 0\). Should such a scenario arise, it is advisable to temporarily replace Equation (34) with Equation (33), which further highlights the importance of exploring alternative choices for \(g(\mathbf {F})\) based on our preceding discussion.

Limitations. An eigenanalysis of \(I_4\) and \(I_6\) (as well as \(V_C\), \({I_8,}\) and \(VI_{C}\)) is still missing, which could lead to a better understanding of materials defined in their terms. Analytic eigensystems are especially relevant due to the implementation flexibility they provide, for example, for parallelization of simulations on the GPU in the style of Huang et al. [2024]. Additionally, our solutions to addressing singularities are robust but of high order (quintic-blend and log-barrier functions). Therefore, resolving singularity will generally require more Newton iterations than the more stable configurations of an FE. For example, the iterations can be up to \(2\times\) more when using the barrier solution, and up to \(10\times\) more with the quintic blend solution in the worst case.

Acknowledgments

We thank the reviewers for their detailed feedback, van Zon [2024] for the Octopus model in used Figure 9, and Qin et al. [2023] for the face model in Figure 13.

Footnotes

Refer, for example, to the technical report of McAdams et al. [2011a] for one example, which applies a modified Jacobi iteration method to compute SVD.

A hyperelastic material is one that is characterized by the property that its strain energy is independent of prior deformation history. In other words, energy is measured as a function of both the current and reference configuration of geometry only.

The potential energy function is also known as the ‘constitutive model.’

⁴

\(f\) here corresponds to the standard textbook quadratic formula, with the plus sign.

⁵

Stated differently, ‘inverting’ a co-dimensional element (e.g., such that \(\sigma _{1} = -\sigma _2\)) is equivalent to a transformation by a full \(\pi\)-rotation in 3D space.

⁶

We do not consider the first lower-order invariant \(I_1 \equiv \mathrm{tr}(\mathbf {R}^T\mathbf {F})\) here because the norm \(\Vert \partial ^2 I_1/\partial \mathbf {F}^2\Vert _F \equiv \Vert \partial \mathbf {R}/\partial \mathbf {F}\Vert _F\) of its Hessian tends toward infinity as \(\alpha \rightarrow 0\), which significantly reduces the magnitude of search directions to negatively impact solver convergence rates.

⁷

This relationship to FEM has generally been unclear [Kim 2020].

⁸

The second term can be expressed following Smith et al. [2019], where \(\det (\mathbf {S}_t)\) evaluates to a constant because \(\mathbf {S}_t\) is independent of \(\mathbf {F}\).

⁹

Isotropic ARAP energy is itself free from issues of singularities, as shown in Appendix C.

¹⁰

A matrix \(\mathbf {M}\) is said to be skew symmetric if it satisfies \(\mathbf {M}+ \mathbf {M}^T = 0\).

¹¹

The product of two symmetric matrices is not necessarily a symmetric matrix.

¹²

Note too that \(I_6(\mathbf {a},\mathbf {b})=\mathbf {a}^T\mathbf {S}\mathbf {b}=\mathbf {b}^T\mathbf {S}\mathbf {a}\) since \(\mathbf {S}\) is symmetric, and \(I_6(\mathbf {a},\mathbf {b})^2\) can be rewritten as \(I_6(\mathbf {a},\mathbf {b})^2=\mathbf {a}^T\mathbf {S}\mathbf {b}\mathbf {b}^T\mathbf {S}\mathbf {a}=\mathbf {b}^T\mathbf {S}\mathbf {a}\mathbf {a}^T\mathbf {S}\mathbf {b}\).

Supplementary Material

tog-23-0120-File003 (tog-23-0120-file003.mp4)

Supplementary material

Download
27.13 MB

A Partial Derivatives of I₄ and I₆

Here we provide the expressions for the partial derivatives that are needed to evaluate the PK1 and Hessian of our low-order formulations for \(I_4\) and \(I_6\).

Notation. For brevity, we will adopt the following shorthand notation,

\begin{align*} f_l &:= \frac{\partial f}{\partial i}, &\quad f_{lm} &:= \frac{\partial ^2 f}{\partial i\partial j}, &\quad f_{lmn} &:= \frac{\partial ^3 f}{\partial i\partial j\partial k},\\ p_{l} &:= \frac{\partial I_4}{\partial i} &\quad p_{lm} &:= \frac{\partial ^2 I_4}{\partial i \partial j}, \\ q_{l} &:= \frac{\partial I_6}{\partial i}, &\quad q_{lm} &:= \frac{\partial ^2 I_6}{\partial i \partial j}, \end{align*}

for the partial derivatives of the root function \(f\), \(I_4,\) and \(I_6\), respectively. The shorthand subscript identifiers \(l,m,n \in \lbrace 1, 2, 3, J, 4, 5, 6, 7\rbrace\) correspond to the invariants \(i,j,k \in \lbrace I_C, II_C, III_C, J, IV_C, V_C, {I_8}, {VI_{C}} \rbrace ,\) in the same sequential order, to give, for example, \(f_3\), \(p_J,\) and \(q_4\) as the shorthand forms for \(\tfrac{\partial f}{\partial III_C}, \tfrac{\partial I_4}{\partial J},\ \text{and}\ \tfrac{\partial I_6}{\partial IV_C},\) respectively.

A.1 \(I_4\) Scalar Expressions When \(\mathbf {F}\in \rm I\!R^{3\times 3}\)

The first-order partial derivatives of \(I_4\) w.r.t. the invariants \(\tfrac{\partial I_4}{\partial i}\) when \(\mathbf {F}\) is a \(3\times 3\) matrix are

\begin{align} p_{1}&=2 f_{11} IV_C+4 f_{21}V_C+ f_{J1}{I_8}, \end{align}

(67)

\begin{align} p_{2}&=2f_{12}IV_C+4f_{22}V_C+f_{J2}{I_8}, \end{align}

(68)

\begin{align} p_{J}&=2f_{1J}IV_C+4f_{2J}V_C+f_{JJ}{I_8}, \end{align}

(69)

\begin{align} p_{4}&=2f_{1}, \end{align}

(70)

\begin{align} p_{5}&=4f_{2}, \end{align}

(71)

\begin{align} p_{6}&=f_{J}. \end{align}

(72)

Similarly, the second-order partial derivatives \(\tfrac{\partial ^2 I_4}{\partial i \partial j}\) are as follows. For \(\tfrac{\partial ^2 I_4}{\partial I_C \partial {i}},\) we have

\begin{align} p_{11}&=2f_{111}IV_C+4f_{211}V_C+f_{J11}{I_8}, \end{align}

(73)

\begin{align} p_{12}&=2f_{112}IV_C+4f_{212}V_C+f_{J12}{I_8}, \end{align}

(74)

\begin{align} p_{1J}&=2f_{11J}IV_C+4f_{21J}V_C+f_{J1J}{I_8}, \end{align}

(75)

\begin{align} p_{14}&=2f_{11}, \end{align}

(76)

\begin{align} p_{15}&=4f_{12}, \end{align}

(77)

\begin{align} p_{16}&=f_{1J}, \end{align}

(78)

for \(\tfrac{\partial ^2 I_4}{\partial II_C \partial {i}}\)

\begin{align} p_{21}&=2f_{121}IV_C+4f_{221}V_C+f_{J21}{I_8}, \end{align}

(79)

\begin{align} p_{22}&=2f_{122}IV_C+4f_{222}V_C+f_{J22}{I_8}, \end{align}

(80)

\begin{align} p_{2J}&=2f_{12J}IV_C+4f_{22J}V_C+f_{J2J}{I_8}, \end{align}

(81)

\begin{align} p_{24}&=2f_{21}, \end{align}

(82)

\begin{align} p_{25}&=4f_{22}, \end{align}

(83)

\begin{align} p_{26}&=f_{2J}, \end{align}

(84)

for \(\tfrac{\partial ^2 I_4}{\partial J \partial {i}}\)

\begin{align} p_{J1}&=2f_{1J1}IV_C+4f_{2J1}V_C+f_{JJ1}{I_8}, \end{align}

(85)

\begin{align} p_{J2}&=2f_{1J2}IV_C+4f_{2J2}V_C+f_{JJ2}{I_8}, \end{align}

(86)

\begin{align} p_{JJ}&=2f_{1JJ}IV_C+4f_{2JJ}V_C+f_{JJJ}{I_8}, \end{align}

(87)

\begin{align} p_{J4}&=2f_{J1}, \end{align}

(88)

\begin{align} p_{J5}&=4f_{J2}, \end{align}

(89)

\begin{align} p_{J6}&=f_{JJ}, \end{align}

(90)

for \(\tfrac{\partial ^2 I_4}{\partial IV_C \partial {i}}\)

\begin{align} p_{41}&=2f_{11}, \end{align}

(91)

\begin{align} p_{42}&=2f_{12}, \end{align}

(92)

\begin{align} p_{4J}&=2f_{1J}, \end{align}

(93)

\begin{align} p_{44} &= p_{45} = p_{46} =0, \end{align}

(94)

for \(\tfrac{\partial ^2 I_4}{\partial V_C \partial {i}}\)

\begin{align} p_{51}&=4f_{21}, \end{align}

(95)

\begin{align} p_{52}&=4f_{22}, \end{align}

(96)

\begin{align} p_{5J}&=4f_{2J}, \end{align}

(97)

\begin{align} p_{54}&=p_{55}=p_{56}=0, \end{align}

(98)

and finally, for \(\tfrac{\partial ^2 I_4}{\partial {I_8} \partial {i}},\) we have

\begin{align} p_{61}&=f_{J1}, \end{align}

(99)

\begin{align} p_{62}&=f_{J2}, \end{align}

(100)

\begin{align} p_{6J}&=f_{JJ}, \end{align}

(101)

\begin{align} p_{64}&=p_{65}=p_{66}=0. \end{align}

(102)

Computing Equation (67) through Equation (102) requires derivatives of the root function \(f\). The first-order partial derivatives

\begin{align*} f_1 = \frac{2f^2+2I_C}{\alpha }, \quad f_2= \frac{-2}{\alpha },\quad f_J = \frac{8f}{\alpha }, \end{align*}

with \(\alpha\) (see Equation (27)) are from Lin et al. [2022], to whom we also refer the reader for the second-order partial derivatives. We also provide the third-order partial derivatives of this \(f\) in Table 3.

Table 3.

Term(s)		Expression
\(f_{111}\)	\(=\)	\(\frac{4f_1^2+4ff_{11}-\left(24ff_1^2+12f^2f_{11}-8f_1-4I_Cf_{11}\right)f_1-2\left(12f^2f_1-4f-4I_Cf_1\right)f_{11}}{\alpha }\)
\(f_{112}=f_{121}=f_{211}\)	\(=\)	\(\frac{4f_1f_2+4ff_{12}-\left(24ff_1f_2+12f^2f_{12}-4f_2-4I_Cf_{12}\right)f_1-\left(12f^2f_1-4f-4I_Cf_1\right)f_{12}-\left(12f^2f_2-4I_Cf_2\right)f_{11}}{\alpha }\)
\(f_{11J}=f_{1J1}=f_{J11}\)	\(=\)	\(\frac{4f_Jf_1+4ff_{1J}-\left(24ff_1f_J+12f^2f_{1J}-4f_J-4I_Cf_{1J}\right)f_1-\left(12f^2f_1-4f-4I_Cf_1\right)f_{1J}-\left(12f^2f_J-4I_Cf_J-8\right)f_{11}}{\alpha }\)
\(f_{122}=f_{212}=f_{221}\)	\(=\)	\(\frac{4f_2^2+4ff_{22}-\left(24ff_2^2+12f^2f_{22}-4I_Cf_{22}\right)f_1-2\left(12f^2f_2-4I_Cf_2\right)f_{12}}{\alpha }\)
\(f_{12J}=f_{1J2}=f_{21J}=f_{2J1}=f_{J12}=f_{J21}\)	\(=\)	\(\frac{4f_Jf_2+4ff_{2J}-\left(24ff_2f_J+12f^2f_{2J}-4I_Cf_{2J}\right)f_1-\left(12f^2f_2-4I_Cf_2\right)f_{1J}-\left(12f^2f_J-4I_Cf_J-8\right)f_{12}}{\alpha }\)
\(f_{1JJ}=f_{J1J}=f_{JJ1}\)	\(=\)	\(\frac{4f_J^2+4ff_{JJ}-\left(24ff_J^2+12f^2f_{JJ}-4I_Cf_{JJ}\right)f_1-2\left(12f^2f_J-4I_Cf_J-8\right)f_{1J}}{\alpha }\)
\(f_{222}\)	\(=\)	\(\frac{-\left(24ff_2^2+12f^2f_{22}-4I_Cf_{22}\right)f_2-2\left(12f^2f_2-4I_Cf_2\right)f_{22}}{\alpha }\)
\(f_{22J}=f_{2J2}=f_{J22}\)	\(=\)	\(\frac{-\left(24ff_2f_J+12f^2f_{2J}-4I_Cf_{2J}\right)f_2-\left(12f^2f_2-4I_Cf_2\right)f_{2J}-\left(12f^2f_J-4I_Cf_J-8\right)f_{22}}{\alpha }\)
\(f_{2JJ}=f_{JJ2}=f_{J2J}\)	\(=\)	\(\frac{-\left(24ff_J^2+12f^2f_{JJ}-4I_Cf_{JJ}\right)f_2-2\left(12f^2f_J-4I_Cf_J-8\right)f_{2J}}{\alpha }\)
\(f_{JJJ}\)	\(=\)	\(\frac{8f_{JJ}-\left(24ff_J^2+12f^2f_{JJ}-4I_Cf_{JJ}\right)f_J-2\left(12f^2f_J-4I_Cf_J-8\right)f_{JJ}}{\alpha }\)

Table 3. Expressions for the Third-Order Partial Derivatives of the Quartic Root Function \(f\) w.r.t. the Invariants

A.2 \(I_4\) Scalar Expressions for 2D Elements

Here we provide the expressions required to evaluate the PK1 and Hessian for 2D elements.

For 2D planar elements, where \(\mathbf {F}\) is a \(2\times 2\) matrix (cf. Section 3.2), we have

\begin{align} p_1&=2f_{11}IV_C+f_{J1}{I_8}, \end{align}

(103)

\begin{align} p_J&=2f_{1J}IV_C+f_{JJ}{I_8}, \end{align}

(104)

\begin{align} p_4&=2f_{1}, \end{align}

(105)

\begin{align} p_6&=f_{J}, \end{align}

(106)

and

\begin{align} f_1&=\frac{1}{2f}, \end{align}

(107)

\begin{align} f_J&=\frac{1}{f}. \end{align}

(108)

Similarly, we have

\begin{align} p_1&=2f_{11}IV_C+2f_{J1}{VI_{C}}, \end{align}

(109)

\begin{align} p_3&=2f_{13}IV_C+2f_{33}{VI_{C}}, \end{align}

(110)

\begin{align} p_4&=2f_{1}, \end{align}

(111)

\begin{align} p_7&=2f_{3}, \end{align}

(112)

and

\begin{align} f_1&=\frac{1}{2f}, \end{align}

(113)

\begin{align} f_3&=\frac{1}{2f\sqrt {III_C}}, \end{align}

(114)

for co-dimensional elements, where \(\mathbf {F}\) is a \(3\times 2\) matrix (cf. Section 3.3). We also evaluate these expressions with

\begin{align} f_{11}&=\frac{-f_1^2}{f}, \end{align}

(115)

\begin{align} f_{1J}&=f_{J1}=\frac{-f_1f_J}{f}, \end{align}

(116)

\begin{align} f_{JJ} &= \frac{-f_J^2}{f}, \end{align}

(117)

\begin{align} f_{13}&=f_{31}=\frac{-f_1f_3}{f}, \end{align}

(118)

\begin{align} f_{33} &= \frac{-\frac{1}{2III_C^{\frac{3}{2}}}-2f_3^2}{2f}, \end{align}

(119)

as second-order partial derivatives of \(f\).

The second-order partial derivatives of \(I_4\) are then

\begin{align} p_{11}&=2f_{111}IV_C+f_{J11}{I_8}, \end{align}

(120)

\begin{align} p_{J1}=p_{1J}&=2f_{11J}IV_C+f_{J1J}{I_8}, \end{align}

(121)

\begin{align} p_{JJ}&=2f_{1JJ}IV_C+f_{JJJ}{I_8}, \end{align}

(122)

\begin{align} p_{41}=p_{14}&=2f_{11}, \end{align}

(123)

\begin{align} p_{61}=p_{16}&=f_{J1}, \end{align}

(124)

\begin{align} p_{4J}=p_{J4}&=2f_{1J}, \end{align}

(125)

\begin{align} p_{6J}=p_{J6}&=f_{JJ}, \end{align}

(126)

for planar elements in 2D. In the case of co-dimensional elements, we have

\begin{align} p_{11}&=2f_{111}IV_C+2f_{311}{VI_{C}}, \end{align}

(127)

\begin{align} p_{31}=p_{13}&=2f_{113}IV_C+2f_{313}{VI_{C}}, \end{align}

(128)

\begin{align} p_{33}&=2f_{133}IV_C+2f_{333}{VI_{C}}, \end{align}

(129)

\begin{align} p_{71}=p_{17}&=2f_{31}, \end{align}

(130)

\begin{align} p_{4J}=p_{J4}&=2f_{13}, \end{align}

(131)

\begin{align} p_{73}=p_{37}&=2f_{33}. \end{align}

(132)

where \(p_{41}=p_{14}\) have the same expression (in terms of \(f\)) as Equation (123). Lastly, we have

\begin{align} f_{111}&=\frac{-3f_1f_{11}}{f}, \end{align}

(133)

\begin{align} f_{11J}=f_{1J1}=f_{J11}&=\frac{-2f_1f_{1J}-f_Jf_{11}}{f}, \end{align}

(134)

\begin{align} f_{1JJ}=f_{J1J}=f_{JJ1}&=\frac{-2f_Jf_{1J}-f_1f_{JJ}}{f}, \end{align}

(135)

\begin{align} f_{JJJ}&=\frac{-3f_Jf_{JJ}}{f}, \end{align}

(136)

\begin{align} f_{113}=f_{131}=f_{311}&=\frac{-2f_1f_{13}-f_3f_{11}}{f}, \end{align}

(137)

\begin{align} f_{133}=f_{313}=f_{331}&=\frac{-2f_3f_{13}-f_1f_{33}}{f}, \end{align}

(138)

\begin{align} f_{333}&=\frac{\frac{3}{4III_C^{\frac{5}{2}}}-6f_3f_{33}}{2f}, \end{align}

(139)

as the third-order partial derivatives of \(f\) that are required to evaluate Equation (120) through Equation (132).

A.3 Matrix Derivatives

The isotropic matrix \(\tfrac{\partial i}{\partial \mathbf {F}}\) and higher-order tensor terms \(\tfrac{\partial ^2 i}{\partial \mathbf {F}^2}\) (where \(i \in \lbrace I_C, II_C, J \rbrace\)) can be found, for example, in the works of Smith et al. [2019] or Kim and Eberle [2022], to whom we refer the reader. We also have

\begin{align} \frac{\partial IV_C}{\partial \mathbf {F}} &= 2\mathbf {F}\mathbf {A}, \quad \frac{\partial ^2 IV_C}{\partial \mathbf {F}^2} = 2\frac{\partial \mathbf {F}}{\partial \mathbf {F}}:\mathbf {A}, \end{align}

(140)

\begin{align} \frac{\partial V_C}{\partial \mathbf {F}} &= 2\frac{\partial \mathbf {C}}{\partial \mathbf {F}}: \mathbf {F}^T\mathbf {F}\mathbf {A}, \end{align}

(141)

\begin{align} \frac{\partial ^2 V_C}{\partial \mathbf {F}^2} &= 2\frac{\partial \mathbf {C}}{\partial \mathbf {F}}:\left(\frac{\partial \mathbf {F}}{\partial \mathbf {F}}:\mathbf {A}\right):\frac{\partial \mathbf {C}}{\partial \mathbf {F}} + 2\frac{\partial ^2\mathbf {C}}{\partial \mathbf {F}^2}:\mathbf {F}^T\mathbf {F}\mathbf {A}, \end{align}

(142)

\begin{align} \frac{\partial {I_8}}{\partial \mathbf {F}} &= \frac{\partial J}{\partial \mathbf {F}}\text{tr}(\mathbf {A}), \quad \frac{\partial ^2 {I_8}}{\partial \mathbf {F}^2} = \frac{\partial ^2 J}{\partial \mathbf {F}^2}\text{tr}(\mathbf {A}), \end{align}

(143)

\begin{align} \frac{\partial {VI_{C}}}{\partial \mathbf {F}} &= \frac{\partial III_C}{\partial \mathbf {F}}\text{tr}(\mathbf {A}), \quad \frac{\partial ^2 {VI_{C}}}{\partial \mathbf {F}^2} = \frac{\partial ^2 III_C}{\partial \mathbf {F}^2}\text{tr}(\mathbf {A}), \end{align}

(144)

as the matrix and fourth-order tensor terms required to evaluate the PK1 and Hessian of \(I_4\), where \(\tfrac{\partial \mathbf {C}}{\partial \mathbf {F}_{ij}} = \tfrac{\partial \mathbf {F}^T}{\partial \mathbf {F}_{ij}}\mathbf {F}+\mathbf {F}^T\tfrac{\partial \mathbf {F}}{\partial \mathbf {F}_{ij}}\), \(\tfrac{\partial ^2 \mathbf {C}}{\partial \mathbf {F}_{ij}\partial \mathbf {F}_{kl}} = \tfrac{\partial \mathbf {F}^T}{\partial \mathbf {F}_{ij}}\tfrac{\partial \mathbf {F}^T}{\partial \mathbf {F}_{kl}}+\tfrac{\partial \mathbf {F}^T}{\partial \mathbf {F}_{kl}}\tfrac{\partial \mathbf {F}}{\partial \mathbf {F}_{ij}}\), and \(\tfrac{\partial \mathbf {F}}{\partial \mathbf {F}_{ij}}\) is a constant matrix with non-zero value 1 only at the i-th row, j-th column.

A.4 The Case of I₆

All scalar, matrix, and tensor terms that are associated with \(I_6\) (i.e., \(q_l\), \(q_{lm}\)) will have the same form/structure as the expressions given for \(I_4\) in preceding sections. The primary distinction is that the anisotropic invariants \(IV_C, V_C, {I_8,}\) and \({VI_{C}}\) will now be evaluated as cross-fiber invariants with two unit vectors \(\mathbf {a}\) and \(\mathbf {b}\) rather than one. For example, the invariant \(V_C(\mathbf {a}) = \mathbf {a}^T\mathbf {C}^T\mathbf {C}\mathbf {a}\) will become \(V_C(\mathbf {a},\mathbf {b}) = \mathbf {b}^T\mathbf {C}^T\mathbf {C}\mathbf {a}\) to give

\begin{equation*} \frac{\partial V_C(\mathbf {a},\mathbf {b})}{\partial \mathbf {F}} = 2\frac{\partial \mathbf {C}}{\partial \mathbf {F}}:\mathbf {F}^T\mathbf {F}\mathbf {D}, \end{equation*}

and

\begin{equation*} \frac{\partial ^2 V_C}{\partial \mathbf {F}^2} =2\frac{\partial \mathbf {C}}{\partial \mathbf {F}}:\left(\frac{\partial \mathbf {F}}{\partial \mathbf {F}}:{\mathbf {D}}\right):\frac{\partial \mathbf {C}}{\partial \mathbf {F}} + 2\frac{\partial ^2\mathbf {C}}{\partial \mathbf {F}^2}:\mathbf {F}^T\mathbf {F}\mathbf {D}, \end{equation*}

as the respective PK1 and Hessian that are derived using the chain rule. Comparing with Equation (142), the matrix \(\mathbf {A}= \mathbf {a}\mathbf {a}^T\) is simply replaced with \(\mathbf {D}=(\mathbf {a}\mathbf {b}^T+\mathbf {b}\mathbf {a}^T)/2\). The same pattern holds for the other cross-fiber invariants.

B PK1 and Hessian of α

In this section, we provide the derivatives of \(\alpha\), which are required to resolve singularity.

To leverage Equation (34) or Equation (33) in the Newton method, it is necessary to compute the first- and second-order partial derivatives of \(\alpha\) w.r.t. \(\mathbf {F}\). When \(\mathbf {F}\) is a \({3\times 3}\) matrix (Section 3.1), these are given by

\begin{align} \frac{\partial \alpha }{\partial \mathbf {F}} &= 4\biggr [(3f^2-I_C)\frac{\partial f}{\partial \mathbf {F}}-\biggr (f\frac{\partial I_C}{\partial \mathbf {F}}+2\frac{\partial J }{\partial \mathbf {F}}\biggr)\biggr ], \end{align}

(145)

\begin{align} \frac{\partial ^2 \alpha }{\partial \mathbf {F}^2} &= 4\biggr [6f\frac{\partial f }{\partial \mathbf {F}}\otimes \frac{\partial f}{\partial \mathbf {F}}+3f^2\frac{\partial ^2 f}{\partial \mathbf {F}^2}- \biggr (\frac{\partial f}{\partial \mathbf {F}}\otimes \frac{\partial I_C}{\partial \mathbf {F}} +f\frac{\partial ^2 I_C}{\partial \mathbf {F}^2} \nonumber \nonumber\\ &+\frac{\partial I_C}{\partial \mathbf {F}}\otimes \frac{\partial f}{\partial \mathbf {F}}+I_C\frac{\partial ^2 f}{\partial \mathbf {F}^2}+2\frac{\partial ^2 J}{\partial \mathbf {F}^2}\biggr) \biggr ]. \end{align}

(146)

When \(\mathbf {F}\) is a \({2\times 2}\) matrix (Section 3.2), we have

\begin{align} \frac{\partial \alpha }{\partial \mathbf {F}} &= \frac{\partial f}{\partial \mathbf {F}} = \mathbf {R}, \end{align}

(147)

\begin{align} \frac{\partial ^2 \alpha }{\partial \mathbf {F}^2} &= \frac{\partial ^2 f}{\partial \mathbf {F}^2} = \frac{\partial \mathbf {R}}{\partial \mathbf {F}}, \end{align}

(148)

which is just the rotation \(\mathbf {R}\) and its (analytic) gradient \(\partial \mathbf {R}/\partial \mathbf {F}\) as the respective PK1 and the Hessian. Finally, when \(\mathbf {F}\) is a \({3\times 2}\) matrix (Section 3.3), we have

\begin{align} \frac{\partial \alpha }{\partial \mathbf {F}} &= \sqrt {III_C}\frac{\partial f}{\partial \mathbf {F}}+f\frac{1}{2\sqrt {III_C}}\frac{\partial III_C}{\partial \mathbf {F}}, \end{align}

(149)

\begin{align} \frac{\partial ^2 \alpha }{\partial \mathbf {F}^2} &= \sqrt {III_C}\frac{\partial ^2 f}{\partial \mathbf {F}^2}+\frac{1}{2\sqrt {III_C}}\frac{\partial III_C}{\partial \mathbf {F}} \otimes \frac{\partial f}{\partial \mathbf {F}} + \frac{1}{2\sqrt {III_C}}\frac{\partial f}{\partial \mathbf {F}} \otimes \frac{\partial III_C}{\partial \mathbf {F}} \nonumber \nonumber\\ & +f\frac{-1}{4III_C^{\frac{3}{2}}}\frac{\partial III_C}{\partial \mathbf {F}} \otimes \frac{\partial III_C}{\partial \mathbf {F}}+f\frac{1}{2\sqrt {III_C}} \frac{\partial ^2 III_C}{\partial \mathbf {F}^2}, \end{align}

(150)

as the expressions required to leverage Equation (34) or Equation (33) in a Newton-type method. Higher-order (tensor) derivatives follow naturally, but these are beyond the scope of our analysis.

C Isotropy and the Elimination of Singularity

Our analysis in Section 4 showed another perspective on how the rotation gradient \(\partial \mathbf {R}/\partial \mathbf {F}\) is the source of singularity in the energy \(\Psi _4\) (cf. Equation (8)), where entries of the PK1 matrix \(\partial \Psi _4/\partial \mathbf {F}\) are scaled inversely proportional to \(\alpha\) which exacerbates anisotropic force magnitudes near problematic configurations of an FE. We supplement our analysis with the results of this section to show how the PK1 of any isotropic energy that is written purely in terms of \(I_1\) is always singularity free. In doing so, we arrive at the main result of this section, which is a constructive proof for the equation \(\mathbf {R}= \partial I_1/\partial \mathbf {F}\), with a further revelation that any non-isotropic model (i.e., anisotropic and orthotropic models) will suffer from singularities.

Elimination of Singularity in Isotropic Energies. The presence or absence of singularity can be demonstrated using Equation (43), from which we can obtain the PK1 in Equation (8) by

\begin{align} \frac{\partial I_1}{\partial \mathbf {F}}= \frac{\partial \mathrm{tr}(\mathbf {F}^T\mathbf {R})}{\partial \mathbf {F}} &=\sum _i^d\frac{\partial I_4(\mathbf {u}_i) }{\partial \mathbf {F}}, \end{align}

(151)

\begin{align} &= \frac{\partial \mathbf {R}}{\partial \mathbf {F}} : \mathbf {F}\left(\sum _i^d \mathbf {u}_i\mathbf {u}_i^T\right) + \mathbf {R}\left(\sum _i^d \mathbf {u}_i\mathbf {u}_i^T\right), \end{align}

(152)

\begin{align} &= \frac{\partial \mathbf {R}}{\partial \mathbf {F}} : \mathbf {F}+ \mathbf {R}, \end{align}

(153)

where \(\mathbf {I}= \sum _i^d \mathbf {u}_i\mathbf {u}_i^T\) is the identity matrix. Singularity is eliminated if and only if \(\partial I_1/\partial \mathbf {F}= \mathbf {R}\), which implies that the first term of Equation (153) must be the zero matrix.

The singularity is eliminated precisely because the summation \(\sum _i\mathbf {u}_i\mathbf {u}_i^T\) (cf. Equation (152)) gives the identity matrix \(\mathbf {I}\) iff there are \(d\) orthonormal vectors \(\mathbf {u}_i\). Otherwise, we arrive, for example, at

\begin{equation} \frac{\partial \mathbf {R}}{\partial \mathbf {F}} : \mathbf {F}\mathbf {A}= \frac{\partial \mathbf {R}}{\partial \mathbf {F}} : \mathbf {R}\mathbf {S}\mathbf {A}= \mathbf {R}^T\frac{\partial \mathbf {R}}{\partial \mathbf {F}} : \mathbf {S}\mathbf {A}\ne 0, \end{equation}

(159)

in the case of anisotropy with \(\mathbf {A}= \mathbf {a}\mathbf {a}^T\) for some vector \(\mathbf {a}\) of unit length (cf. Equation (8)). Moreover, Equation (159) is non-zero because \(\mathbf {S}\mathbf {A}\) is not a symmetric matrix.¹¹

D Decoupling Stretch and Shear

This section provides supplementary proofs for showing that the isotropic invariant \(I_2\) may be expressed in terms of \(I_4\) and \(I_6\), which we use to arrive at Equation (51) and its generalization to 3D.

2D Case. Suppose that the low-order FEM formulation of the Baraff-Witkin model is indeed

\begin{align} \Psi _{i{\rm\small ARAP}}&=(I_4(\mathbf {a})-1)^2+(I_4(\mathbf {b})-1)^2+2I_6(\mathbf {a},\mathbf {b})^2, \end{align}

(160)

as stated in Equation (50) to give

\begin{align} \Psi _{i{\rm\small ARAP}} &\equiv (I_4(\mathbf {a})^2-2I_4(\mathbf {a})+1)+(I_4(\mathbf {b})^2-2I_4(\mathbf {b})+1)+2I_6(\mathbf {a},\mathbf {b})^2 , \end{align}

(161)

\begin{align} &=I_4(\mathbf {a})^2+I_4(\mathbf {b})^2+2I_6(\mathbf {a},\mathbf {b})^2-2(I_4(\mathbf {a})+I_4(\mathbf {b}))+2, \end{align}

(162)

\begin{align} &=I_4(\mathbf {a})^2+I_4(\mathbf {b})^2+2I_6(\mathbf {a},\mathbf {b})^2-2I_1+2, \end{align}

(163)

where Equation (163) has been simplified using \(I_4(\mathbf {a})+I_4(\mathbf {b})=I_1\) (from Equation (42)). Then, Equation (160) is the anisotropic and cross-fiber-invariant-based expression of isotropic ARAP energy iff

\begin{equation} {I_4(\mathbf {a})^2+I_4(\mathbf {b})^2+2I_6(\mathbf {a},\mathbf {b})^2 + 2} = {I_2} \end{equation}

(164)

holds since Equation (163) would then become equivalent to Equation (40).

Proof.

Expanding out the quadratic terms of Equation (164), we have

\begin{align} &I_2 = I_4(\mathbf {a})^2+I_4(\mathbf {b})^2+2I_6(\mathbf {a},\mathbf {b})^2 , \end{align}

(165)

\begin{align} &=\mathbf {a}^T\mathbf {S}\mathbf {a}\mathbf {a}^T\mathbf {S}\mathbf {a}+\mathbf {b}^T\mathbf {S}\mathbf {b}\mathbf {b}^T\mathbf {S}\mathbf {b}+\mathbf {a}^T\mathbf {S}\mathbf {b}\mathbf {b}^T\mathbf {S}\mathbf {a}+\mathbf {b}^T\mathbf {S}\mathbf {a}\mathbf {a}^T\mathbf {S}\mathbf {b}, \end{align}

(166)

\begin{align} &=\text{tr}(\mathbf {a}^T\mathbf {S}\mathbf {a}\mathbf {a}^T\mathbf {S}\mathbf {a})+\text{tr}(\mathbf {b}^T\mathbf {S}\mathbf {b}\mathbf {b}^T\mathbf {S}\mathbf {b})+\text{tr}(\mathbf {a}^T\mathbf {S}\mathbf {b}\mathbf {b}^T\mathbf {S}\mathbf {a})+\text{tr}(\mathbf {b}^T\mathbf {S}\mathbf {a}\mathbf {a}^T\mathbf {S}\mathbf {b}), \end{align}

(167)

\begin{align} &=\text{tr}(\mathbf {S}\mathbf {a}\mathbf {a}^T\mathbf {S}\mathbf {a}\mathbf {a}^T)+\text{tr}(\mathbf {S}\mathbf {b}\mathbf {b}^T\mathbf {S}\mathbf {b}\mathbf {b}^T)+\text{tr}(\mathbf {S}\mathbf {b}\mathbf {b}^T\mathbf {S}\mathbf {a}\mathbf {a}^T)+\text{tr}(\mathbf {S}\mathbf {a}\mathbf {a}^T\mathbf {S}\mathbf {b}\mathbf {b}^T), \end{align}

(168)

\begin{align} &=\text{tr}(\mathbf {S}\mathbf {a}\mathbf {a}^T\mathbf {S}\mathbf {a}\mathbf {a}^T+\mathbf {S}\mathbf {b}\mathbf {b}^T\mathbf {S}\mathbf {b}\mathbf {b}^T+\mathbf {S}\mathbf {b}\mathbf {b}^T\mathbf {S}\mathbf {a}\mathbf {a}^T+\mathbf {S}\mathbf {a}\mathbf {a}^T\mathbf {S}\mathbf {b}\mathbf {b}^T), \end{align}

(169)

\begin{align} &=\text{tr}(\mathbf {S}(\underbrace{\mathbf {a}\mathbf {a}^T\mathbf {S}\mathbf {a}\mathbf {a}^T+\mathbf {b}\mathbf {b}^T\mathbf {S}\mathbf {b}\mathbf {b}^T+\mathbf {b}\mathbf {b}^T\mathbf {S}\mathbf {a}\mathbf {a}^T+\mathbf {a}\mathbf {a}^T\mathbf {S}\mathbf {b}\mathbf {b}^T}_{\mathbf {S}_{\mathbf {a}\mathbf {b}}})), \end{align}

(170)

which we get by using the cyclic property of the trace in Equations (167) and (168).¹² Using the symbol \(\mathbf {S}_{\mathbf {a}\mathbf {b}}\) to represent the terms of the innermost parentheses of Equation (170), we have

\begin{align} \mathbf {S}_{\mathbf {a}\mathbf {b}}&=\mathbf {a}\mathbf {a}^T\mathbf {S}\mathbf {a}\mathbf {a}^T+\mathbf {b}\mathbf {b}^T\mathbf {S}\mathbf {a}\mathbf {a}^T+\mathbf {a}\mathbf {a}^T\mathbf {S}\mathbf {b}\mathbf {b}^T+\mathbf {b}\mathbf {b}^T\mathbf {S}\mathbf {b}\mathbf {b}^T, \end{align}

(171)

\begin{align} &= (\mathbf {a}\mathbf {a}^T\mathbf {S}+\mathbf {b}\mathbf {b}^T\mathbf {S})(\mathbf {a}\mathbf {a}^T+\mathbf {b}\mathbf {b}^T), \end{align}

(172)

\begin{align} &=(\mathbf {a}\mathbf {a}^T\mathbf {S}+\mathbf {b}\mathbf {b}^T\mathbf {S})(\mathbf {a}\mathbf {a}^T+\mathbf {b}\mathbf {b}^T), \end{align}

(173)

\begin{align} &=(\mathbf {a}\mathbf {a}^T+\mathbf {b}\mathbf {b}^T)\mathbf {S}(\mathbf {a}\mathbf {a}^T+\mathbf {b}\mathbf {b}^T), \end{align}

(174)

\begin{align} &=\mathbf {U}^T\mathbf {S}\mathbf {U} \end{align}

(175)

\begin{align} &=\mathbf {S}, \end{align}

(176)

where \(\ \mathbf {U}\ =\ (\mathbf {a}\mathbf {a}^T+\mathbf {b}\mathbf {b}^T)\) is an orthonormal matrix to give

\begin{align} \text{tr}(\mathbf {S}({\mathbf {S}_{\mathbf {a}\mathbf {b}}})) &= \text{tr}(\mathbf {S}^2) = I_2, \end{align}

(177)

as the reduced form of Equation (170). □

Equation (177) thus proves that Equation (163) is indeed the anisotropic and cross-fiber-invariant-based expression of isotropic ARAP energy.

3D Case. A similar approach can be applied in 3D if we suppose that

\begin{align} \Psi _{i{\rm\small ARAP}}&=(I_4(\mathbf {a})-1)^2+(I_4(\mathbf {b})-1)^2+(I_4(\mathbf {c})-1)^2\nonumber \nonumber\\ &+2I_6(\mathbf {a},\mathbf {b})^2+2I_6(\mathbf {a},\mathbf {c})^2+2I_6(\mathbf {b},\mathbf {c})^2 \end{align}

(178)

is the anisotropic and cross-fiber-invariant-based expression of isotropic ARAP energy, where \(\mathbf {a}, \mathbf {b}, \mathbf {c}\in \rm I\!R^{3\times 1}\) are three unit-length vectors that we assume to be orthogonal.

Proof.

Expanding out the terms in Equation (178), we get

\begin{align} \Psi _{i{\rm\small ARAP}}&\equiv I_4(\mathbf {a})^2+I_4(\mathbf {b})^2+I_4(\mathbf {c})^2+2I_6(\mathbf {a},\mathbf {b})^2+2I_6(\mathbf {a},\mathbf {c})^2 \nonumber \nonumber\\ &\quad + 2I_6(\mathbf {b},\mathbf {c})^2-2(I_4(\mathbf {a})+I_4(\mathbf {b})+I_4(\mathbf {c}))+3, \end{align}

(179)

\begin{align} &=I_4(\mathbf {a})^2+I_4(\mathbf {b})^2+I_4(\mathbf {c})^2+2I_6(\mathbf {a},\mathbf {b})^2+2I_6(\mathbf {a},\mathbf {c})^2 \nonumber \nonumber\\ &\quad +2I_6(\mathbf {b},\mathbf {c})^2-2I_1+3, \end{align}

(180)

which is simplified using \(I_4(\mathbf {a})+I_4(\mathbf {b})+I_4(\mathbf {c})=I_1\) (from Equation (42)). We likewise gather and expand all quadratic terms in Equation (180) to posit

\begin{align} I_2 &= I_4(\mathbf {a})^2+I_4(\mathbf {b})^2+I_4(\mathbf {c})^2+2I_6(\mathbf {a},\mathbf {b})^2+2I_6(\mathbf {a},\mathbf {c})^2+2I_6(\mathbf {b},\mathbf {c})^2, \end{align}

(181)

\begin{align} &=\text{tr}(\mathbf {S}(\mathbf {S}_{\mathbf {a}\mathbf {b}\mathbf {c}})), \end{align}

(182)

where

\begin{align} \mathbf {S}_{\mathbf {a}\mathbf {b}\mathbf {c}} &= \mathbf {a}\mathbf {a}^T\mathbf {S}\mathbf {a}\mathbf {a}^T+\mathbf {b}\mathbf {b}^T\mathbf {S}\mathbf {b}\mathbf {b}^T+\mathbf {c}\mathbf {c}^T\mathbf {S}\mathbf {c}\mathbf {c}^T+ \mathbf {b}\mathbf {b}^T\mathbf {S}\mathbf {a}\mathbf {a}^T+\mathbf {a}\mathbf {a}^T\mathbf {S}\mathbf {b}\mathbf {b}^T \nonumber \nonumber\\ &\quad +\mathbf {c}\mathbf {c}^T\mathbf {S}\mathbf {a}\mathbf {a}^T+\mathbf {a}\mathbf {a}^T\mathbf {S}\mathbf {c}\mathbf {c}^T+ \mathbf {b}\mathbf {b}^T\mathbf {S}\mathbf {c}\mathbf {c}^T+\mathbf {c}\mathbf {c}^T\mathbf {S}\mathbf {b}\mathbf {b}^T , \end{align}

(183)

\begin{align} &= \mathbf {U}^T\mathbf {S}\mathbf {U}, \end{align}

(184)

\begin{align} &= \mathbf {S}, \end{align}

(185)

with \(\mathbf {U}= (\mathbf {a}\mathbf {a}^T+\mathbf {b}\mathbf {b}^T+\mathbf {c}\mathbf {c}^T)\) being an orthonormal matrix. □

From Equation (182), we arrive at a form resembling Equation (177) to show that Equation (178) is the anisotropic and cross-fiber-invariant-based expression of isotropic ARAP energy in 3D. Equation (178) thus also presents an avenue for a decoupling strategy similar to Equation (51) but in 3D volumetric elements.

References

[1]

Marc Alexa, Daniel Cohen-Or, and David Levin. 2000. As-rigid-as-possible shape interpolation. In Proceedings of the 27th Annual Conference on Computer Graphics and Interactive Techniques (SIGGRAPH’00). 157–164. DOI:

Digital Library

Google Scholar

[2]

David Baraff and Andrew Witkin. 1998. Large steps in cloth simulation. In Proceedings of the 25th Annual Conference on Computer Graphics and Interactive Techniques (SIGGRAPH’98). ACM, New York, NY, USA, 43–54. DOI:

Digital Library

Google Scholar

[3]

Silvia S. Blemker and Scott L. Delp. 2005. Three-dimensional representation of complex muscle architectures and geometries. Annals of Biomedical Engineering 33, 5 (May 2005), 661–673. DOI:

Crossref

Google Scholar

[4]

Silvia S. Blemker, Peter M. Pinsky, and Scott L. Delp. 2005. A 3D model of muscle reveals the causes of nonuniform strains in the biceps brachii.Journal of Biomechanics 38 4 (2005), 657–65.

Crossref

Google Scholar

[5]

J. Bonet and A. J. Burton. 1998. A simple orthotropic, transversely isotropic hyperelastic constitutive equation for large strain computations. Computer Methods in Applied Mechanics and Engineering 162, 1 (1998), 151–164. DOI:

Crossref

Google Scholar

[6]

Javier Bonet, Antonio J. Gil, and Rogelio Ortigosa. 2015. A computational framework for polyconvex large strain elasticity. Computer Methods in Applied Mechanics and Engineering 283 (2015), 1061–1094.

Crossref

Google Scholar

[7]

Javier Bonet, Antonio J. Gil, and Rogelio Ortigosa. 2016. On a tensor cross product based formulation of large strain solid mechanics. International Journal of Solids and Structures 84 (2016), 49–63.

Crossref

Google Scholar

[8]

Javier Bonet and Richard D. Wood. 2008. Nonlinear Continuum Mechanics for Finite Element Analysis (2nd ed.). Cambridge University Press. DOI:

Crossref

Google Scholar

[9]

Grégory Chagnon, Marie Rebouah, and Denis Favier. 2015. Hyperelastic energy densities for soft biological tissues: A review. Journal of Elasticity 120 (2015), 129–160.

Crossref

Google Scholar

[10]

Isaac Chao, Ulrich Pinkall, Patrick Sanan, and Peter Schröder. 2010. A simple geometric model for elastic deformations. ACM Transactions on Graphics 29, 4 (July 2010), Article 38, 6 pages. DOI:

Digital Library

Google Scholar

[11]

Tao Du, Kui Wu, Pingchuan Ma, Sebastien Wah, Andrew Spielberg, Daniela Rus, and Wojciech Matusik. 2021. DiffPD: Differentiable projective dynamics. ACM Transactions on Graphics 41, 2 (Nov. 2021), Article 13, 21 pages.

Google Scholar

[12]

Ugo Finnendahl, Matthias Schwartz, and Marc Alexa. 2023. ARAP revisited discretizing the elastic energy using intrinsic Voronoi cells. Computer Graphics Forum 42, 6 (2023), e14790. DOI:

Crossref

Google Scholar

[13]

Zhan Gao, Theodore Kim, Doug L. James, and Jaydev P. Desai. 2009. Semi-automated soft-tissue acquisition and modeling for surgical simulation. In Proceedings of the 5th Annual IEEE International Conference on Automation Science and Engineering(CASE’09). IEEE, 268–273.

Digital Library

Google Scholar

[14]

Nicholas J. Higham. 2008. Functions of Matrices. Society for Industrial and Applied Mathematics. DOI:

Crossref

Google Scholar

[15]

Kemeng Huang, Floyd M. Chitalu, Huancheng Lin, and Taku Komura. 2024. GIPC: Fast and stable Gauss-Newton optimization of IPC barrier energy. ACM Transactions on Graphics 43, 2 (March 2024), Article 23, 18 pages. DOI:

Digital Library

Google Scholar

[16]

Alexandru-Eugen Ichim, Petr Kadleček, Ladislav Kavan, and Mark Pauly. 2017. Phace: Physics-based face modeling and animation. ACM Transactions on Graphics 36, 4 (July 2017), Article 153, 14 pages.

Digital Library

Google Scholar

[17]

G. Irving, J. Teran, and R. Fedkiw. 2004. Invertible finite elements for robust simulation of large deformation. In Proceedings of the 2004 ACM SIGGRAPH/Eurographics Symposium on Computer Animation(SCA’04). 131–140. DOI:

Digital Library

Google Scholar

[18]

G. Irving, J. Teran, and R. Fedkiw. 2006. Tetrahedral and hexahedral invertible finite elements. Graphical Models 68, 2 (2006), 66–89. DOI:

Digital Library

Google Scholar

[19]

Theodore Kim. 2020. A finite element formulation of Baraff-Witkin cloth. In Proceedings of the ACM SIGGRAPH/Eurographics Symposium on Computer Animation(SCA’20). Article 16, 9 pages.

Digital Library

Google Scholar

[20]

Theodore Kim, Fernando De Goes, and Hayley Iben. 2019. Anisotropic elasticity for inversion-safety and element rehabilitation. ACM Transactions on Graphics 38, 4 (July 2019), Article 69, 15 pages.

Digital Library

Google Scholar

[21]

Theodore Kim and David Eberle. 2022. Dynamic deformables: Implementation and production practicalities (now with code!). In ACM SIGGRAPH 2022 Courses(SIGGRAPH’22). ACM, New York, NY, USA, Article 7, 259 pages.

Digital Library

Google Scholar

[22]

Gergely Klár, Andrew Moffat, Ken Museth, and Eftychios Sifakis. 2020a. Shape targeting: A versatile active elasticity constitutive model. In ACM SIGGRAPH 2020 Talks(SIGGRAPH’20). ACM, New York, NY, USA, Article 59, 2 pages.

Digital Library

Google Scholar

[23]

Gergely Klár, Andrew Moffat, Ken Museth, and Eftychios Sifakis. 2020b. Shape targeting: A versatile active elasticity constitutive model (supplemental material). In ACM SIGGRAPH 2020 Talks(SIGGRAPH’20). ACM, New York, NY, USA, Article 59, 2 pages.

Digital Library

Google Scholar

[24]

Seunghwan Lee, Ri Yu, Jungnam Park, Mridul Aanjaneya, Eftychios Sifakis, and Jehee Lee. 2018. Dexterous manipulation and control with volumetric muscles. ACM Transactions on Graphics 37, 4 (July 2018), Article 57, 13 pages.

Digital Library

Google Scholar

[25]

Minchen Li, Zachary Ferguson, Teseo Schneider, Timothy Langlois, Denis Zorin, Daniele Panozzo, Chenfanfu Jiang, and Danny M. Kaufman. 2020. Incremental potential contact: Intersection- and inversion-free large deformation dynamics. ACM Transactions on Graphics 39, 4 (2020), Article 49, 20 pages.

Digital Library

Google Scholar

[26]

Huancheng Lin, Floyd M. Chitalu, and Taku Komura. 2022. Isotropic ARAP energy using Cauchy-Green invariants. ACM Transactions on Graphics 41, 6 (Nov. 2022), Article 275, 14 pages.

Digital Library

Google Scholar

[27]

Ligang Liu, Lei Zhang, Yin Xu, Craig Gotsman, and Steven J. Gortler. 2008. A local/global approach to mesh parameterization. In Proceedings of the Symposium on Geometry Processing(SGP’08). 1495–1504.

Digital Library

Google Scholar

[28]

Jerrold E. Marsden and Thomas J. R. Hughes. 1994. Mathematical Foundations of Elasticity. Dover Publications.

Google Scholar

[29]

Aleka McAdams, Andrew Selle, Rasmus Tamstorf, Joseph Teran, and Eftychios Sifakis. 2011a. Computing the Singular Value Decomposition of 3x3 Matrices with Minimal Branching and Elementary Floating Point Operations. Technical Report. University of Wisconsin, Madison. http://graphics.cs.wisc.edu/Papers/2011/MSTTS11

Google Scholar

[30]

Aleka McAdams, Yongning Zhu, Andrew Selle, Mark Empey, Rasmus Tamstorf, Joseph Teran, and Eftychios Sifakis. 2011b. Efficient elasticity for character skinning with contact and collisions. ACM Transactions on Graphics 30, 4 (July 2011), Article 37, 12 pages. DOI:

Digital Library

Google Scholar

[31]

D. G. Mead. 1992. Newton’s identities. American Mathematical Monthly 99, 8 (1992), 749–751.

Crossref

Google Scholar

[32]

Théodore Papadopoulo and Manolis I. A. Lourakis. 2000. Estimating the Jacobian of the singular value decomposition: Theory and applications. In Computer Vision—ECCV 2000. Springer, Berlin, Germany, 554–570.

Crossref

Google Scholar

[33]

Ken Perlin. 2002. Improving noise. ACM Transactions on Graphics 21, 3 (July 2002), 681–682. DOI:

Digital Library

Google Scholar

[34]

Guillaume Picinbono, Herve Delingette, and Nicholas Ayache. 2000. Real-time large displacement elasticity for surgery simulation: Non-linear tensor-mass model. In Proceedings of the 3rd International Conference on Medical Image Computing and Computer-Assisted Intervention (MICCAI’00). 643–652.

Crossref

Google Scholar

[35]

Roman Poya, Rogelio Ortigosa, and Theodore Kim. 2023. Geometric optimisation via spectral shifting. ACM Transactions on Graphics 42, 3 (April 2023), Article 29, 15 pages.

Digital Library

Google Scholar

[36]

Dafei Qin, Jun Saito, Noam Aigerman, Thibault Groueix, and Taku Komura. 2023. Neural face rigging for animating and retargeting facial meshes in the wild. In ACM SIGGRAPH 2023 Conference Proceedings(SIGGRAPH’23). ACM, New York, NY, USA, Article 68, 11 pages. DOI:

Digital Library

Google Scholar

[37]

Christian Schüller, Ladislav Kavan, Daniele Panozzo, and Olga Sorkine-Hornung. 2013. Locally injective mappings. Computer Graphics Forum 32, 5 (2013), 125–135.

Digital Library

Google Scholar

[38]

Eftychios Sifakis and Jernej Barbic. 2012. FEM simulation of 3D deformable solids: A practitioner’s guide to theory, discretization and model reduction. In ACM SIGGRAPH 2012 Courses(SIGGRAPH’12). ACM, New York, NY, USA, Article 20, 50 pages.

Digital Library

Google Scholar

[39]

Eftychios Sifakis, Igor Neverov, and Ronald Fedkiw. 2005. Automatic determination of facial muscle activations from sparse motion capture marker data. ACM Transactions on Graphics 24, 3 (July 2005), 417–425. DOI:

Digital Library

Google Scholar

[40]

F. S. Sin, D. Schroeder, and J. Barbič. 2013. Vega: Non-linear FEM deformable object simulator. Computer Graphics Forum 32, 1 (2013), 36–48.

Crossref

Google Scholar

[41]

Breannan Smith, Fernando De Goes, and Theodore Kim. 2018. Stable neo-Hookean flesh simulation. ACM Transactions on Graphics 37, 2 (March 2018), Article 12, 15 pages.

Digital Library

Google Scholar

[42]

Breannan Smith, Fernando De Goes, and Theodore Kim. 2019. Analytic eigensystems for isotropic distortion energies. ACM Transactions on Graphics 38, 1 (Feb. 2019), Article 3, 15 pages.

Digital Library

Google Scholar

[43]

Olga Sorkine and Marc Alexa. 2007. As-rigid-as-possible surface modeling. In Proceedings of the 5th Eurographics Symposium on Geometry Processing(SGP’07). 109–116.

Google Scholar

[44]

Olga Sorkine-Hornung and Michael Rabinovich. 2016. Least-Squares Rigid Motion Using SVD. Technical Note. ETH Zurich.

Google Scholar

[45]

Sangeetha Grama Srinivasan, Qisi Wang, Junior Rojas, Gergely Klár, Ladislav Kavan, and Eftychios Sifakis. 2021. Learning active quasistatic physics-based models from data. ACM Transactions on Graphics 40, 4 (July 2021), Article 129, 14 pages.

Digital Library

Google Scholar

[46]

Alexey Stomakhin, Russell Howes, Craig Schroeder, and Joseph M. Teran. 2012. Energetically consistent invertible elasticity. In Proceedings of the ACM SIGGRAPH/Eurographics Symposium on Computer Animation(SCA’12). 25–32.

Digital Library

Google Scholar

[47]

J. Teran, S. Blemker, V. Ng Thow Hing, and R. Fedkiw. 2003. Finite volume methods for the simulation of skeletal muscle. In Proceedings of the 2003 ACM SIGGRAPH/Eurographics Symposium on Computer Animation(SCA’03). 68–74.

Digital Library

Google Scholar

[48]

Christopher D. Twigg and Zoran Kačić-Alesić. 2010. Point cloud glue: Constraining simulations using the procrustes transform. In Proceedings of the 2010 ACM SIGGRAPH/Eurographics Symposium on Computer Animation (SCA’10). 45–54. DOI:

Crossref

Google Scholar

[49]

Dirk van Zon. 2024. Octopus Rig. Retrieved June 2, 2024 from https://sketchfab.com/3d-models/octopus-rig-48f40cb3dfd5473fbb200cf1dcd7b944

Google Scholar

[50]

Jeffrey A. Weiss, Bradley N. Maker, and Sanjay Govindjee. 1996. Finite element implementation of incompressible, transversely isotropic hyperelasticity. Computer Methods in Applied Mechanics and Engineering 135, 1 (1996), 107–128. DOI:

Crossref

Google Scholar

[51]

Martin Welk, Joachim Weickert, and Gabriele Steidl. 2006. From tensor-driven diffusion to anisotropic wavelet shrinkage. In Proceedings of the 9th European Conference on Computer Vision—Volume Part I(ECCV’06). 391–403. DOI:

Digital Library

Google Scholar

[52]

Lingchen Yang, Byungsoo Kim, Gaspard Zoss, Baran Gözcü, Markus Gross, and Barbara Solenthaler. 2022. Implicit neural representation for physics-driven actuated soft bodies. ACM Transactions on Graphics 41, 4 (July 2022), Article 122, 10 pages.

Digital Library

Google Scholar

Index Terms

Analytic rotation-invariant modelling of anisotropic finite elements

Index terms have been assigned to the content through auto-classification.

Recommendations

Anisotropic finite strain viscoelasticity based on the Sidoroff multiplicative decomposition and logarithmic strains

In this paper a purely phenomenological formulation and finite element numerical implementation for quasi-incompressible transversely isotropic and orthotropic materials is presented. The stored energy is composed of distinct anisotropic equilibrated ...
Orthotropic rotation-free thin shell elements

A method to simulate orthotropic behaviour in thin shell finite elements is proposed. The approach is based on the transformation of shape function derivatives, resulting in a new orthogonal basis aligned to a specified preferred direction for all ...
Insight into 3-node triangular shell finite elements: the effects of element isotropy and mesh patterns

In this paper, we study the convergence characteristics of some 3-node triangular shell finite elements. We review the formulations of three different isotropic 3-node elements and one non-isotropic 3-node element. We analyze a clamped plate problem and ...

Comments

Information & Contributors

Information

Published In

cover image ACM Transactions on Graphics

ACM Transactions on Graphics Volume 43, Issue 5

October 2024

207 pages

EISSN:1557-7368

DOI:10.1145/3613708

Editor:
Carol O'Sullivan
Trinity College Dublin, Ireland

Issue’s Table of Contents

This work is licensed under a Creative Commons Attribution International 4.0 License.

Publisher

Association for Computing Machinery

New York, NY, United States

Publication History

Published: 09 August 2024

Online AM: 28 May 2024

Accepted: 14 May 2024

Revised: 07 May 2024

Received: 18 October 2023

Published in TOG Volume 43, Issue 5

Check for updates

Author Tags

Qualifiers

Research-article

Funding Sources

RGC grant
Innovation and Technology Commission of the HKSAR Government under the InnoHK initiative (TransGP project)
JC STEM Lab of Robotics for Soft Materials funded by the Hong Kong Jockey Club Charities Trust

Contributors

Other Metrics

View Article Metrics

Bibliometrics & Citations

Bibliometrics

Article Metrics

0
Total Citations
1,737
Total Downloads

Downloads (Last 12 months)1,737
Downloads (Last 6 weeks)177

Reflects downloads up to 25 Feb 2025

Other Metrics

View Author Metrics

Citations

View Options

View options

PDF

View or Download as a PDF file.

PDF

eReader

View online with eReader.

eReader

Login options

Check if you have access through your login credentials or your institution to get full access on this article.

Abstract

1 Introduction

2 Background

2.1 The Derivatives of Arbitrary I4 Energies

2.2 Inversion-Safe Anisotropic ARAP Energy

2.3 Challenges with \(\Psi_AA\)-Energy Minimization

3 Expressing I4 in Proper Closed form

3.1 Volumetric FEs in 3D

3.2 Planar FEs in 2D

3.3 Co-dimensional Elements

4 Fixing the Singularity Problem

4.1 The Sources of Singularity

4.2 Two Solutions for Fixing Singularity

5 Applications and Generalizations

5.1 Anisotropic ARAP: Revisited

5.2 Cloth

5.3 Shape Targeting

6 Results

7 Conclusion and Discussion

7.1 Discussion

Acknowledgments

Footnotes

Supplementary Material

A Partial Derivatives of I4 and I6

A.1 \(I_4\) Scalar Expressions When \(\mathbf {F}\in \rm I\!R^{3\times 3}\)

A.2 \(I_4\) Scalar Expressions for 2D Elements

A.3 Matrix Derivatives

A.4 The Case of I6

B PK1 and Hessian of α

C Isotropy and the Elimination of Singularity

D Decoupling Stretch and Shear

References

Index Terms

Recommendations

Anisotropic finite strain viscoelasticity based on the Sidoroff multiplicative decomposition and logarithmic strains

Orthotropic rotation-free thin shell elements

Insight into 3-node triangular shell finite elements: the effects of element isotropy and mesh patterns

Comments

Information

Published In

Publisher

Publication History

Check for updates

Author Tags

Qualifiers

Funding Sources

Contributors

Other Metrics

Bibliometrics

Article Metrics

Other Metrics

Citations

View options

PDF

eReader

Login options

Full Access

Share

Share this Publication link

Share on social media

Affiliations

2.1 The Derivatives of Arbitrary I₄ Energies

3 Expressing I₄ in Proper Closed form

A Partial Derivatives of I₄ and I₆

A.4 The Case of I₆