Stable Neo-Hookean Flesh Simulation

12
Stable Neo-Hookean Flesh Simulation
BREANNAN SMITH, FERNANDO DE GOES, and THEODORE KIM, Pixar Animation Studios
Fig. 1. Left: Thirteen skeletal bones drive a hexahedral lattice with 45,809 elements and 156,078 degrees of freedom. Center: A quasi-static simulation with
our new Neo-Hookean model and a Poisson’s ratio of ν = 0.488. Wrinkles and bulges emerge from our model’s excellent volume-preserving properties. An
average time step with our model took 13.7 Newton iterations, 5,860 Conjugate Gradient (CG) iterations, and 25.6 seconds. Right: The same simulation
with corotational elasticity and ν = 0.488. The model fails to preserve volume and instead collapses the trapezius and forms a spurious fold around the
shoulder blade. The artifacts persist across all values of ν . An average time step with this model took 17.9 Newton iterations, 16,183 CG iterations, and 46.6
seconds. ©Disney/Pixar.
Nonlinear hyperelastic energies play a key role in capturing the fleshy appearance
of virtual characters. Real-world, volume-preserving biological
tissues have Poisson’s ratios near 1/2, but numerical simulation within this
regime is notoriously challenging. In order to robustly capture these visual
characteristics, we present a novel version of Neo-Hookean elasticity.
Our model maintains the fleshy appearance of the Neo-Hookean model,
exhibits superior volume preservation, and is robust to extreme kinematic
rotations and inversions. We obtain closed-form expressions for the eigenvalues
and eigenvectors of all of the system’s components, which allows us
to directly project the Hessian to semipositive definiteness, and also leads
to insights into the numerical behavior of the material. These findings also
inform the design of more sophisticated hyperelastic models, which we
explore by applying our analysis to Fung and Arruda-Boyce elasticity. We
provide extensive comparisons against existing material models.
CCS Concepts: • Computing methodologies → Physical simulation;
Additional Key Words and Phrases: Physically-based simulation, elasticity
ACM Reference format:
Breannan Smith, Fernando de Goes, and Theodore Kim. 2018. Stable NeoHookean
Flesh Simulation. ACM Trans. Graph. 37, 2, Article 12 (March
2018), 15 pages.
https://doi.org/10.1145/3180491
Authors’ addresses: B. Smith, F. de Goes, and T. Kim, 1200 Park Avenue, Emeryville,
CA 94608; emails: {breannan, fernando, tkim}@pixar.com.
Permission to make digital or hard copies of all or part of this work for personal or
classroom use is granted without fee provided that copies are not made or distributed
for profit or commercial advantage and that copies bear this notice and the full citation
on the first page. Copyrights for components of this work owned by others than
the author(s) must be honored. Abstracting with credit is permitted. To copy otherwise,
or republish, to post on servers or to redistribute to lists, requires prior specific
permission and/or a fee. Request permissions from permissions@acm.org.
2018 Copyright is held by the owner/author(s). Publication rights licensed to ACM.
ACM 0730-0301/2018/03-ART12 $15.00
https://doi.org/10.1145/3180491
1 INTRODUCTION
The elastic energy used to model deformations determines the visual
quality of a simulation, so it must be chosen carefully for virtual
characters and virtual humans. The defining quality of biological
tissues such as muscle and fat is volume preservation, which is
reflected in their high Poisson’s ratios ν ∈ [0.45, 0.5) (Greaves et al.
2011) and the corresponding visual features that emerge (Figure 2).
However, this near-incompressible regime is notoriously difficult
to simulate robustly and accurately. For instance, the popular corotational
model for elasticity (Chao et al. 2010; McAdams et al. 2011;
Müller et al. 2002; Nesme et al. 2009) ostensibly supports materials
within this range, but it linearizes the volume term in a way that
compromises both its volumetric properties and its visual results
(Figure 1, right). Alternatively, force-filtering approaches (Irving
et al. 2004; Teran et al. 2005) inherit the visual quality of the NeoHookean
model but introduce additional user parameters.
In this work, we propose a novel Neo-Hookean energy that retains
the rich, volumetric character of biological materials and does
not need any filter parameters. While motivated by biological tissues
where ν ≈ 0.5, our model behaves well across a wide range
of Poisson’s ratios (ν ∈ [0, 0.5)). Our model exhibits better volume
preservation than equivalent filtered models and does not break
down under the large kinematic rotations that often occur in quasistatic
simulations. The computational performance of our method
is competitive with, and often better than, existing methods. We
know of no other model that possesses all of these advantages.
We achieve these properties by performing an analysis that
completely characterizes the Hessian of the energy and results in
closed-form expressions for the eigenvectors and eigenvalues of
each of its terms. The eigendecomposition of the overall system
can then be written as a set of compact expressions, which we
ACM Transactions on Graphics, Vol. 37, No. 2, Article 12. Publication date: March 2018.
12:2 • B. Smith et al.
Fig. 2. Volume preservation is the key to convincing flesh simulations. We simulate an elbow flexion over a range of Poisson’s ratios ν with our stable
Neo-Hookean model. Smaller values of ν resemble a skinning, and it is not until ν ≈ 1/2 that we obtain significant volumetric bulging in the forearm.
©Disney/Pixar.
use to isolate the sources of indefiniteness and project them back
to semidefiniteness. As a consequence, we can safely apply our
energy to Newton-type implicit integration schemes that use
conjugate gradient-based linear solvers. While this analysis was
motivated by our Neo-Hookean model, it is sufficiently generic
that it also supports any model formulated in terms of the first
and third invariants of the deformation gradient. As preliminary
examples, we show how it can be extended to Fung (2013) and
Arruda-Boyce (1993) elasticity. We conclude with detailed comparisons
to existing hyperelastic energies. In particular, we are
able to phrase the popular corotational model as a linearization of
the Neo-Hookean energy (Section 5.1) and position our model on
a spectrum of successive approximations to this energy.
2 BACKGROUND
Before presenting our material model, we provide a brief review
of existing hyperelastic energies and numerical simulators. In the
following, we represent scalars as unbolded letters (x, J), vectors
as bold lowercase letters (x), and matrices as bold uppercase letters
(X). Greek letters (Ψ, α) represent scalars.
2.1 Deformation Gradient
The fundamental measure of nonlinear deformation is the deformation
gradient F. Here, we label its columns with vectors and its
entries with scalars:
F =
⎡
⎢
⎢
⎢
⎢
⎢
⎣
f0 f1 f2
⎤
⎥
⎥
⎥
⎥
⎥
⎦
=
⎡
⎢
⎢
⎢
⎢
⎢
⎣
f0 f3 f6
f1 f4 f7
f2 f5 f8
⎤
⎥
⎥
⎥
⎥
⎥
⎦
. (1)
Table 1 summarizes the derived quantities of F we will use. One
factorization requires special attention: the rotation variant SVD
(Higham 2008; Irving et al. 2004; Twigg and Kačić-Alesić 2010) of
F, which we write as
F = UΣVT and Σ =
⎡
⎢
⎢
⎢
⎢
⎢
⎣
σ0 0 0
0 σ1 0
0 0 σ2
⎤
⎥
⎥
⎥
⎥
⎥
⎦
. (2)
Unlike the standard SVD convention, the rotation variant moves
reflections to Σ, so U and V are rotations with det(U) = det(V) = 1,
and Σ is allowed to have a negative entry.
2.2 Hyperelastic Energies
The elastic behavior of a deformable body can be specified in terms
of a hyperelastic energy density Ψ, such as the corotational (CR)
Table 1. Quantities Derived from the Deformation Gradient F
Symbol Definition
F = RS Polar decomposition
J = det(F) Relative volume change
C = FT F Right Cauchy-Green tensor
IC = tr (C) First right Cauchy-Green invariant
material (McAdams et al. 2011):
ΨCR = μ||F − R||2
F +
λ
2
tr2 (S − I) . (3)
Here, μ and λ are the Lamé constants, R and S form the polar decomposition
of F (Table 1), and I is a 3 × 3 identity matrix. Forces
can be derived from this energy by several means, but a popular
approach uses the first Piola-Kirchhoff (PK1) stress tensor, which
we denote P(F). The PK1 for the corotational material is
PCR(F) = ∂Ψ
∂F = R [μ(S − I) + λ tr(S − I)I] . (4)
After integrating over a volume V and changing variables to vertex
displacements u, we obtain a matrix,V P(F) ∂F
∂u , whose columns
are vertex forces. The term ∂F
∂u varies based on the element type,
e.g., tetrahedra or hexahedra, and basis type, e.g., linear or quadratic
(Bargteil and Cohen 2014). Other texts describe this derivation
(Müller et al. 2008; Sifakis and Barbic 2012), so we will deal
directly with Ψ and P(F) and defer to these references for further
details.
2.3 Robustness Challenges
A key drawback of nonlinear hyperelastic energies is that they
usually contain singularities. Early work on deformable simulation
(Hirota et al. 2001) detected these “illegal” states and carefully
avoided them using a backtracking line search. These searches are
challenging to perform robustly, because a NaN in a single degree
of freedom can quickly grow to pollute larger regions of the
simulation. This backtracking approach also precludes simulations
where the solution contains inverted elements, which often occurs
in production environments.
Irving et al. (2004) described an explicit method for simulating
nonlinear materials that instead used the rotation variant SVD
(Equation (2)) to write a filtered stress tensor P(Σ) defined in F’s
principal stretch space Σ. For specific values of Σ where P(Σ)
was known to behave badly, the tensor was clamped to constant
ACM Transactions on Graphics, Vol. 37, No. 2, Article 12. Publication date: March 2018.
Stable Neo-Hookean Flesh Simulation • 12:3
or linear functions. Thus, no backtracking was needed to avoid
illegal states. In order to use this method, the user identifies badly
behaved regions of P(Σ) and sets a filtering threshold. Notably,
the usual convention for the rotation variant SVD is to push reflections
into the smallest singular value, but subsequent work has
shown that this can introduce subtle problems, including discontinuities
and nonoptimal inversion recovery directions (Georgii
and Westermann 2008; Schmedding and Teschner 2008). The work
of Civit-Flores and Susín (2014) provides an excellent overview.
Implicit integration of the force-filtering approach was described
by Teran et al. (2005). Since the fourth-order tensor ∂P(F)
∂F
may become indefinite, a separate derivation was needed to guarantee
its positive definiteness. The tensor was diagonalized into
one 3 × 3 and three 2 × 2 matrices, and eigenvalues were clamped
to a threshold. This approach has the limitation that Ψ, P(F) and
∂P(F)
∂F can drift out of sync. This issue was alleviated by Stomakhin
et al. (2012) by applying the filter directly to Ψ. Xu et al. (2015) proposed
an alternative where users design force curves in principal
stretch space. We instead present an energy that is agnostic to any
reflection convention, even in the presence of inverted elements.
Alternative deformation schemes such as position-based dynamics
(Bender et al. 2015; Umetani et al. 2014) and projective dynamics
(Bouaziz et al. 2014) have recently emerged. One advantage
of these schemes is that the Hessian ∂P(F)
∂F does not need to
be constructed explicitly. However, projective dynamics still need
to specify filtering thresholds and reflection conventions (Narain
et al. 2016; Wang and Yang 2016) to avoid numerical singularities.
The advantage that our model requires neither carries over
directly.
3 ENERGY FORMULATION
Following the preceding discussion, we design a new hyperelastic
energy that is stable in four important ways:
• Inversion stability: The energy is singularity-free and does
not need any filters or thresholds.
• Reflection stability: The energy is well behaved regardless
of the reflection convention used in the SVD.
• Rest stability: Under zero load, an element’s rest shape is
preserved.
• Meta-stability under degeneracy: Under point, line, and
plane degeneracies, the forces are defined up to rotation.
We start by examining existing hyperelastic energies that go by
the name “Neo-Hookean.”
3.1 Existing Neo-Hookean Energies
The most common version of Neo-Hookean elasticity is (see,
e.g., Bonet and Wood (2008))
ΨNeo = μ
2
(IC − 3) − μ log J +
λ
2
(log J)
2. (5)
Here, IC and J are defined in Table 1. However, many other versions
appear in the literature:
ΨA = μ
2
(IC − 3) − μ log J +
λ
2
(J − 1)
2 [Ogden 1997],
ΨB = μ
2
(J−2/3
IC − 3) +
λ
2
(J − 1)
2 [Bower 2009],
ΨC = μ
2
(J−2/3
IC − 3) +
λ
2
(J − 1) [Wang and Yang 2016].
We examine these energies under the Valanis-Landel hypothesis
(Xu et al. 2015), which posits that many hyperelastic energies can
be separated into length (1D), area (2D), and volume (3D) components.
The above models contain length and volume terms, but no
area terms. Adding an area term yields a Mooney-Rivlin model.
1D Length Term: Mooney (1940) originally proposed the
energy
ΨM = μ
2
(IC − 3), (6)
which was later dubbed the “Neo-Hookean” energy by Rivlin
(1948). When unconstrained, this energy achieves its minimum
when the element has collapsed to zero volume, i.e., when IC = 0
and ΨM = −3. Mooney additionally imposed the hard constraint
that J = 1, so the energy is instead minimized at the volumepreserving
configuration that is the closest to the stretch space
origin. Even without this constraint, we note that ΨM is well behaved
under inversion. The energy relative to a zero-volume configuration
is always well defined irrespective of an element’s current
state.
In contrast, the ΨB and ΨC energies use a modified Neo-Hookean
term μ/2(J−2/3
IC − 3) (Rivlin 1948). The J−2/3 term factors off the
nearest squared isotropic stretch from IC, leaving only the squared
distortional part of the deformation. It also introduces numerical
problems by growing without bound under compression, i.e., as
J → 0, and becoming undefined at J = 0. Therefore, we prefer to
use a ΨM length term in lieu of a modified energy.
3D Volume Term: The remaining terms in the above energies
are volume-preserving penalty terms. The volume terms from
ΨNeo clearly present numerical difficulties:
ΨNeo,volume = −μ log J +
λ
2
(log J)
2. (7)
In addition to growing unbounded as J → 0, they become infeasible
for J < 0. These difficulties motivated the design of the forcefiltering
approach (Irving et al. 2004). Any energy with a logarithm,
including ΨA,volume = −μ log J + λ/2(J − 1)
2, will also have these
problems.
In contrast, ΨB,volume = λ/2(J − 1)
2 does not have these issues
and is bounded, well defined, and invertible. Such terms have been
used (Martin et al. 2011; Teschner et al. 2004) to sidestep the need
for any inversion handling. A promising candidate is then
ΨD = μ
2
(IC − 3) +
λ
2
(J − 1)
2, (8)
but it has a key drawback that we must address.
3.2 Rest Stabilization
While ΨD possesses inversion stability, it lacks what we will call
rest stability. In the solid mechanics literature (e.g., Bonet and
Wood (2008), Section 6.4.3), it is often mentioned that the hyperelastic
energy should vanish at identity. This occurs with ΨD, but it
is more important that the PK1 resolve to zero, as this is the true
indicator that the energy has been extremized. Otherwise, forces
appear when the body is at rest and overwrite the intended rest
state with a different, parameter-dependent state. This is a known
ACM Transactions on Graphics, Vol. 37, No. 2, Article 12. Publication date: March 2018.
12:4 • B. Smith et al.
side effect of certain barrier functions (Schüller et al. 2013). One solution
is to increase the volume penalty constant λ dramatically in
order to reduce the visual appearance of the artifact (e.g., Blemker
et al. (2005) suggests ∼1,000μ). We avoid this approach, because it
precludes the artifact-free simulation of materials with lower Poisson’s
ratios. Scaling λ also increases the indefiniteness of the Hessian
(see Section 4.6) and should be done judiciously.
Intuitively, the problem is that the effective rest state no
longer coincides with J = 1. Instead, it coincides with a smaller
J that causes an element to slightly deflate. Consequently, we ask
whether we can inflate each element so that, when it deflates, it
settles to rest at J = 1. To answer this question, we consider the
PK1 of ΨD:
PD(F) = μF + λ ∂J
∂F (J − 1). (9)
An energy is rest stable if PD(I) = 0 but PD(I) is nonzero:
PD(I) = μI + λ ∂ det(I)
∂F (det(I) − 1) = μI. (10)
To recover rest stability, we modify (J − 1)
2 to shift the root from
1 to α; i.e., we inflate the element to a volume greater than one:
⎧⎪⎪⎪⎪
⎨
⎪⎪⎪⎪
⎩
ΨE = μ
2
(IC − 3) +
λ
2
(J − α)
2
PE(F) = μF + λ ∂J
∂F (J − α).
Solving for an α that satisfies PE(I) = 0 yields α = 1 + μ/λ:
ΨE = μ
2
(IC − 3) +
λ
2

J − 1 − μ
λ
2
. (11)
We can additionally expand the quadratic to obtain
ΨE = μ
2
(IC − 3) − μ(J − 1) +
λ
2
(J − 1)
2 +
 μ
λ
2
. (12)
Since constants disappear under differentiation, this material
model is functionally equivalent to
ΨE = μ
2
(IC − 3) − μ(J − 1) +
λ
2
(J − 1)
2. (13)
This is strikingly similar to the original ΨNeo from Equation (5),
except the log J terms have been replaced with (J −1). We make
two important observations. First, ΨNeo is rest stable, but not
singularity-free. Second, (J −1) is the first term in the Taylor series
of log J at J =1. Essentially, we have performed a singularityremoving
Taylor truncation of ΨNeo about a volume-preserving
state, resulting in an inversion and rest-stable material model. We
can also verify that our energy ΨE is reflection stable. The IC term
involves the squared singular values of F, so any negation convention
is irrelevant. Since the J term is the product of the singular
values, the sign convention is again irrelevant. Similar reasoning
can be used for PK1, and we show in Section 4.5.1 that the Hessian
is agnostic to the convention as well.
3.3 Meta-Stability under Degeneracy
The energy now has inversion, reflection, and rest stability. However,
its behavior must be examined under degeneracy, i.e., when
an element has been crushed to a plane, line, or point. This examination
can also be viewed as a Drucker stability analysis (Bower
2009).
Fig. 3. To stress-test our stable Neo-Hookean (ν = 0.49) model, we randomly
scramble a cube’s vertices within a space of twice its volume. The
cube recovers its rest pose, highlighting our model’s robustness under extreme,
inverted configurations.
An extensive treatment of all three degeneracies is given in the
supplemental materials, but we will sketch the main results here.
For the case where an element is crushed to a plane, the energy
is stable. Material forces appear along the normal direction to the
plane and work to restore the original shape.
When the element is crushed to a line, the energy is meta-stable.
The forces resolve to zero, but the Hessian is negative definite.
For a material that is invertible, this meta-stable state is correct.
If an unconstrained element has been crushed to exactly a line,
the shape it should return to is underdetermined, because the configuration
is equidistant from all possible rotations of the original
shape about the line. However, negative definiteness guarantees
that any perturbation (e.g., momentum) will select a rotation, and
the element will self-restore. The alternative to meta-stability is a
singularity at this configuration, which is clearly undesirable.
When an element is crushed to a point, i.e., F = 0, the forces
disappear and the Hessian is positive definite, which means the
configuration will remain stable even in the presence of perturbations.
This behavior is undesirable, though the basin of attraction
is vanishingly small for Poisson’s ratios that correspond to biological
tissue (ν ≥ 0.45). For completeness, we introduce one more
term to eliminate the basin across a wide range of ν.
We add a regularized origin barrier log (IC + δ ) that inserts a
peak of negative definiteness about F = 0 but uses δ to smooth
away the logarithmic singularity. In the supplemental material, we
show that δ = 1 is the value that exactly produces semidefiniteness
at F = 0. Our final energy is then
Ψnew = μ
2
(IC − 3) +
λ
2
(J − α)
2 − μ
2 log (IC + 1). (14)
As we derive in the supplemental material, the rest-stability term
now shifts to α = 1 + μ/λ − (μ/4) λ.
3.4 Lamé Reparameterization
With this energy, we now consider the reparameterization of the
material parameters μ and λ in order to maintain consistency with
Hooke’s law. To this end, we shift the values of μ and λ so that our
model reproduces the PK1 of linear elasticity, which is of the form
P(F) = 2μLaméϵ + λLamé tr (ϵ )I, (15)
where ϵ =1/2(F + FT ) − I is the linearized strain tensor, and the coefficients
μLamé and λLamé are the Lamé parameters in linear elasticity.
The linearization of the stress derived from Equation (14)
ACM Transactions on Graphics, Vol. 37, No. 2, Article 12. Publication date: March 2018.
Stable Neo-Hookean Flesh Simulation • 12:5
becomes consistent with this expression if we set μ = 4/3μLamé and
λ = λLamé + 5/6μLamé. The expression for Poisson’s ratio shifts to
ν = λ − (5/8) μ
2(λ + (1/8) μ)
. (16)
Note that the same methodology can be applied to any of the previous
energies. In particular, the reparameterization of Equation (13)
yields μ = μLamé, λ = λLamé + μLamé, and ν = (λ−μ)/(2λ).
Using this new expression for Poisson’s ratio (Equation (16)), we
can analyze our energy to confirm that it does not introduce any
spurious minima over the range ν ∈ [0, 0.5). The details are again
given in the supplemental materials, and we will sketch the results
here. We are able to show symbolically that the critical points of
Ψnew (including maxima and saddle points) only appear under uniform
scaling. Using techniques similar to that employed for metastability,
we analyze this deformation mode and find that spurious
minima only appear when λ/μ < 0.1525. This corresponds to
the auxetic regime of ν ≤ −0.85, so our model is stable for any
ν ∈ [0, 0.5).
4 ENERGY EIGENANALYSIS
We show that it is possible to perform a complete eigenanalysis
of Equation (14). In particular, we construct closed-form eigenvalues
and eigenvectors of each of the terms in the Hessian, and
obtain a compact expression for the eigendecomposition of their
sum. These expressions will be used to construct positive semidefinite
versions of the Hessian, and also to develop a qualitative
understanding of the energy. We begin by specifying a notation
for fourth-order tensors.
4.1 Tensor Notation
There are many nth order tensor notations (e.g. (Kolda and Bader
2009; Simmonds 2012)), but we will specialize to 4th order tensors.
Similar to Golub and Van Loan (2012), we define vectorization
vec(·) as column-wise flattening of a matrix into a vector:
A =

a c
b d 

vec(A) =
⎡
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎣
a
b
c
d
⎤
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎦
. (17)
Since vec(A) is a vector, we could denote it as a. However, we
will add the symbol ˇ· to indicate that it is a vectorized matrix,
i.e. vec(A) = aˇ. We denote 4th order tensors in matrix-of-matrices
form using blackboard bold:
A =
⎡
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎣

a c
b d 
 i k
j l 


e д
f h 
 m o
n p 

⎤
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎦
=

[A00] [A01]
[A10] [A11]


.
When vec(·) is applied, two unfoldings reorder the 4th order tensor
into a 2nd order matrix:
vec(A) =
⎡
⎢
⎢
⎢
⎢
⎢
⎣
vec(A00) vec(A10) vec(A01) vec(A11)
⎤
⎥
⎥
⎥
⎥
⎥
⎦
=
⎡
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎣
aeim
bf jn
cдk o
dhl p
⎤
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎦
= Aˇ .
To emphasize that a matrix is a vectorized (“flattened”) tensor, we
again use the ˇ· notation, as shown in the Aˇ above. This vectorization
convention differs from Teran et al. (2005), which clusters
the diagonal entries of the submatrices in A. Our convention will
allow us to write several expressions in terms of cross products.
4.2 First Piola-Kirchhoff Stress (PK1)
We start from the PK1 for Equation (14),
P(F) = μ

1 − 1
IC + 1

F + λ(J − α)
∂J
∂F
, (18)
where α = 1 + μ/λ − (μ/4) λ. We omit the subscript, as we only consider
one model in this section. Using the column-wise notation
for F (Equation (1)) and the identity J = f0 · (f1 × f2), we write ∂J
∂F
(a.k.a. the cofactor matrix) as cross products:
∂J
∂F =
⎡
⎢
⎢
⎢
⎢
⎢
⎣
f1 × f2 f2 × f0 f0 × f1
⎤
⎥
⎥
⎥
⎥
⎥
⎦
. (19)
This is a convenient shorthand for computing ∂J
∂F , and will be useful
when analyzing ∂2 J
∂F2 .
4.3 The Energy Hessian Terms
Using the scalar notation for F (see Equation (1)), we can write the
Hessian of the energy in fourth-order matrix-of-matrices form:
∂2Ψ
∂F2 = ∂P(F)
∂F =
⎡
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎣

∂P(F)
∂f0
  ∂P(F)
∂f3
  ∂P(F)
∂f6


∂P(F)
∂f1
  ∂P(F)
∂f4
  ∂P(F)
∂f7


∂P(F)
∂f2
  ∂P(F)
∂f5
  ∂P(F)
∂f8

⎤
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎦
, (20)
where each entry is defined as
∂P(F)
∂fi
= μ

1 − 1
IC + 1
 ∂F
∂fi 	
Ti
+ μ
2
(IC + 1)2 Ffi
	
Mi
+ λ ∂J
∂F
∂J
∂fi 	
Gi
+λ(J − α) ∂2 J
∂F∂fi 	
Hi
,
(21)
and each ∂P(F)
∂fi ∈ R3×3. Since we will examine each of the terms
in detail, we have tagged them as Ti , Mi , Gi , and Hi , respectively
the Tikhonov, Mu, volume Gradient, and volume Hessian terms.
The complete Hessian is then their sum:
∂2Ψ
∂F2 = μ

1 − 1
IC + 1

T + μ
2
(IC + 1)2 M(F)
+ λG(F) + λ(J − α)H(F).
(22)
We will examine the eigensystem of each term, and then build the
system that results from their sum.
ACM Transactions on Graphics, Vol. 37, No. 2, Article 12. Publication date: March 2018.
12:6 • B. Smith et al.
4.4 The Tikhonov, Mu, and Gradient Terms
We start with the Tikhonov term. Since ∂F
∂fi is zero except in the
ith entry, this becomes an identity matrix under vectorization:
vec (T) = Tˇ = I, (23)
where I denotes a 9 × 9 identity matrix. This matrix is full rank,
positive definite, and independent of the values in F, and serves as
a diagonal (Tikhonov) regularizer for the rest of the energy.
For the Mu term, we vectorize F into a nine-vector, vec(F) = ˇ
f,
which allows the tensor to be written as an outer product:
vec (M) = Mˇ = ˇ
fˇ
f
T. (24)
This rank-one matrix has a single nonzero eigenvalue:
ˇ
f2
2 ≡ ||F||2
F ≡ 

σ2
0 + σ2
1 + σ2
2

, (25)
where ·F denotes the Frobenius norm and {σi} the singular values
from Equation (2). The eigenvector is ˇ
f up to normalization: ˇ
f/ˇ
f. The eigenvalue is always nonnegative, and if F contains a
large stretch, the eigenvalue will be large as well.
The gradient term has a similar structure. Using the vectorization
vec( ∂J
∂F ) = gˇ, the tensor becomes an outer product:
vec (G) = Gˇ = gˇgˇT. (26)
Again, there is a single nonzero, nonnegative eigenvalue:
gˇ2
2 ≡





∂J
∂F





2
F
≡ 
(σ0σ1)
2 + (σ0σ2)
2 + (σ1σ2)
2

. (27)
The eigenvector is gˇ up to normalization: gˇ/gˇ.
4.5 The Volume Hessian
The Hi term is more involved. Vectorizing Hi reveals the structure
vec(H) = Hˇ =
⎡
⎢
⎢
⎢
⎢
⎢
⎢
⎣
0 −f
2 f
1
f
2 0 −f
0
−f
1 f
0 0
⎤
⎥
⎥
⎥
⎥
⎥
⎥
⎦
, (28)
where the symbol · denotes a cross-product matrix:
x =
⎡
⎢
⎢
⎢
⎢
⎢
⎣
0 −x2 x1
x2 0 −x0
−x1 x0 0
⎤
⎥
⎥
⎥
⎥
⎥
⎦
. (29)
The Hˇ matrix forms a self-similar cross-product matrix. Equation
(28) contains cross-product matrices, and its blocks form
a macro-structure that is a cross-product matrix. Similar crossproduct-like
structures will continue to appear later.
4.5.1 Volume Hessian Eigenvalues. The eigenvalues of Hˇ can be
factored into two characteristic polynomials:
ε3 − tr (C)ε − 2J = 0 (30)
ε3 − tr (C)ε2 +
1
2


tr2 (C) − tr (C2)

ε − det(C) = 0, (31)
where C is from Table 1, and ε denotes the eigenvalues of Hˇ . The
usual λ notation for eigenvalues has been put aside and replaced
by ε because λ is already reserved as a material constant.
Equation (31) is simpler to solve, as it corresponds to the characteristic
polynomial of C. Given its roots (εα ,εβ ,εγ ), six of the
Fig. 4. The set of volume-preserving configurations in 2D and 3D stretch
space are shown in red. The energy’s saddle-point structure is visible in
2D.
eigenvalues of Hˇ are (±
√εα , ±√εβ , ±√εγ ). The square roots of
eigenvalues can be restated using the singular values of F:
ε3 = σ0 ε6 = −σ0
ε4 = σ1 ε7 = −σ1
ε5 = σ2 ε8 = −σ2.
(32)
If the SVD of F has been computed, then six of the eigenvalues of
Hˇ are already known.
For the remaining three eigenvalues, Equation (30) is a depressed
cubic whose roots can be obtained directly:
εk = 2
IC
3 cos 
1
3

arccos 
3J
IC
 3
IC

+ 2πk

k = 0, 1, 2. (33)
We now have all of the eigenvalues of Hˇ . The negations in Equation
(32) and the IC and J terms in Equation (33) render the reflection
convention applied to the singular values irrelevant.
There are always negative eigenvalues in the volume Hessian.
Three of the six values ε3...8 are less than or equal to zero, and
the cosine function in ε0...2 guarantees that another one or two
eigenvalues are negative. Moreover, we have determined that the
other eigenvalues in ∂2Ψ
∂F2 are nonnegative, so the volume Hessian
is the only possible source of negative eigenvalues.
We investigate further by examining the geometry of J = det(F)
in principal stretch space (Ogden 1997; Stomakhin et al. 2012; Xu
et al. 2015). The 2D case is easier to visualize, and the intuition
carries over to 3D. As shown in red (Figure 4(a)), the set of volumepreserving
configurations forms two distinct, disconnected sets.
This reflects the nonconvexity of J = det(F) (Liu et al. 2016), and
its saddle-point structure is visible in the contours. The red curve
closely resembles the “primary contour” of Stomakhin et al. (2012).
The red points are energetically equivalent, so points outside
this curve have no preferred Newton direction, even when the energy
is squared, e.g., (J − 1)
2. Additional regularization is needed
to resolve this ambiguity, which is the role played by the other
terms in ∂2Ψ
∂F2 . The situation worsens in 3D, as the set of volumepreserving
configurations forms four components (Figure 4(b)).
The intuition remains that saddle points exist everywhere and
need to be regularized.
The structure we see in Equations (32) and (33) extends the understanding
of the findings in Teran et al. (2005). There, they factor
ACM Transactions on Graphics, Vol. 37, No. 2, Article 12. Publication date: March 2018.
Stable Neo-Hookean Flesh Simulation • 12:7
the Hessian into one 3 × 3 and three 2 × 2 eigensystems. This corresponds
exactly to the three entangled roots we found for Equation
(33) and the three pair-wise roots we found for Equation (32).
This underscores the fact that our current analysis can be applied
to any system that uses J = det(F) and is not specific to the energy
Ψvolume = (J − 1)
2.
4.5.2 Volume Hessian Eigenvectors. The eigenvectors of Hˇ can
be written in terms of Qˇ from the eigendecomposition Hˇ = Qˇ ΛQˇ T.
However, we prefer the tensor form:
Q =
⎡
⎢
⎢
⎢
⎢
⎢
⎣
[Q0] [Q3] [Q6]
[Q1] [Q4] [Q7]
[Q2] [Q5] [Q8]
⎤
⎥
⎥
⎥
⎥
⎥
⎦
. (34)
Each eigenvector is then a 3 × 3 matrix instead of a nine-vector.
We start with the eigenvectors corresponding to Equation (31).
The eigenvector corresponding to ε3 yields an elegant structure:
Q3 = 1
√
2
UD3VT (35)
D3 =
⎡
⎢
⎢
⎢
⎢
⎢
⎣
0 00
0 01
0 −1 0
⎤
⎥
⎥
⎥
⎥
⎥
⎦
, (36)
where U and V are from the SVD of F (Equation (2)) and the
1/
√
2 is a normalization factor. This eigenvector has a “pseudocross-product”
structure that we have not seen anywhere else. The
middle matrix D3 is a cross-product matrix for the basis vector
e0 = [−1 0 0]T . In fact, for the special case where U=V, Equation
(35) becomes a standard cross-product on the first column of
U, i.e., Q3 = 1/
√
2u0. Equation (35) has a more general form where
U and V differ, and corresponds to an infinitesimal rotation along
the σ0 axis in stretch space.
The complementary eigenvector for ε6, has a similar structure:
Q6 = 1
√
2
UD6VT (37)
D6 =
⎡
⎢
⎢
⎢
⎢
⎢
⎣
000
001
010
⎤
⎥
⎥
⎥
⎥
⎥
⎦
. (38)
Note that the negative has been removed from the last row. D3
has been multiplied by a reflection, making Q6 a reflected pseudocross-product.
This corresponds to an infinitesimal rotation along
the antipodal axis to D3. The eigenvectors for Q4, Q5, Q7, and Q8
follow analogously as the cross-products for e1 = [0 − 1 0]T and
e2 = [0 0 − 1]T . We explicitly list all six eigenvectors in Appendix
A.
The eigenvectors for the depressed cubic (Equation (30)) are also
defined by U and V but contain a conventional diagonal:
Qk = 1
qk
UDkVT k = 0...2 (39)
Dk =
⎡
⎢
⎢
⎢
⎢
⎢
⎣
σ0σ2 + σ1εk 0 0
0 σ1σ2 + σ0εk 0
0 0 ε2
k − σ2
2
⎤
⎥
⎥
⎥
⎥
⎥
⎦
, (40)
where qk = Dk F is a normalization factor. A Mathematica notebook
Verify_Eigenvectors.nb that symbolically confirms this
eigensystem is included in the supplemental materials.
4.6 The Complete Eigensystem
With these closed-form expressions for the individual terms,
we can analyze the complete system, which we denote Aˇ for
Assembled. The vectorized version of Equation (22) is
Aˇ = μ

1 − 1
IC + 1

I + μ
2
(IC + 1)2
ˇ
fˇ
f
T + λgˇgˇT + λ(J − α)Hˇ . (41)
Deriving expressions for the the eigenvalues of a sum of matrices is
generally nontrivial (Knutson and Tao 2001). However, our energy
has a special structure that allows us to obtain simple expressions
for this system.
Since the μ(1 − 1
IC +1 )I term is diagonal, it has nine identical
eigenvalues that can be directly added to the eigenvalues of any existing
system. The main difficulty arises from computing the sum
μ 2
(IC +1)2 ˇ
fˇ
fT + λgˇgˇT + λ(J − α)Hˇ . However, both ˇ
fˇ
fT and gˇgˇT are
rank-one matrices, which allows the problem to be viewed as a
rank-two update to the eigensystem for λ(J − α)Hˇ .
The problem can be simplified using deflation (Bunch et al.
1978). To this end, we show in Appendix B that both gˇ and ˇ
f
are orthogonal to the six pseudo-cross-product eigenvectors of Hˇ .
As a consequence, the eigenvectors of Hˇ and their corresponding
eigenvalues remain unchanged, and the rank-two update has now
been simplified to the effect solely on the depressed cubic (Equations
(33) and (39)). The system is sufficiently small that direct
symbolic examination becomes tractable. The eigenvalues of the
rank-two updated matrix are the roots of
ε¯
3 + c2ε¯
2 + c1ε¯ + c0 = 0 (42)
c2 = −


gˇ



2
2
ρ − ICη
c1 = −(1 + 2Jρ)IC − 6Jη +




gˇ



2
2
IC − 9J 2

ρη
c0 = −(2 + 3Jρ)J +

I
2
C − 4



gˇ



2
2

η + 2J

I
2
C − 3



gˇ



2
2

ρη,
where η and ρ are defined as
η = 2μ
(IC + 1)
2 

λ (J − 1) − 3
4 μ
 ρ = λ
λ (J − 1) − 3
4 μ
. (43)
The cubic is no longer depressed, and while we could apply the
cubic formula to obtain a closed-form expression, we found that
solving for the roots using the Jenkins and Traub (1970) rpoly algorithm
was more practical. The form of the eigenvectors remains
the same, Qk = 1
qk UDkVT , but the diagonal entries differ,
Dk =
⎡
⎢
⎢
⎢
⎢
⎢
⎣
α0 0 0
0 α1 0
0 0 α2
⎤
⎥
⎥
⎥
⎥
⎥
⎦
, (44)
ACM Transactions on Graphics, Vol. 37, No. 2, Article 12. Publication date: March 2018.
12:8 • B. Smith et al.
and instead take the form
α0 = ε¯k (σ1 + σ0σ2η + Jσ1ρ)
+ σ0σ2 + σ1 (σ2
0 − σ2
1 + σ2
2 )η + Jσ0σ2ρ
+ σ0 (σ2
0 − σ2
1 )σ2 (σ2
1 − σ2
2 )ρη
α1 = ε¯k (σ0 + σ1σ2η + Jσ0ρ)
+ σ1σ2 − σ0 (σ2
0 − σ2
1 − σ2
2 )η + Jσ1σ2ρ
− σ1 (σ2
0 − σ2
1 )σ2 (σ2
0 − σ2
2 )ρη
α2 = ε¯
2
k − ε¯k (σ2
0 + σ2
1 )(η + σ2
2 ρ)
− σ2
2 − 2Jη − 2Jσ2
2 ρ + ((σ2
0 − σ2
1 )σ2)
2ρη.
Final eigensystem: With these expressions for the rank-two
updated system, we can add the μT = μ (1 − 1/(IC + 1)) regularization
term and arrive at the final eigenvalues. The first three are
εk = λ(J − α)ε¯k + μT k = 0...2, (45)
where ε¯k are the roots of Equation (42). The other six are
ε3 = λ(J − α)σ0 + μT ε6 = −λ(J − α)σ0 + μT (46)
ε4 = λ(J − α)σ1 + μT ε7 = −λ(J − α)σ1 + μT (47)
ε5 = λ(J − α)σ2 + μT ε8 = −λ(J − α)σ2 + μT . (48)
The first three eigenvectors are specified by Equations (39) and
(44), and the last six remain the ones described in Appendix A.
5 DISCUSSION AND EXTENSION
In this section, we compare our energy to corotational elasticity
and show the generality of our analysis by applying it to Fung and
Arruda-Boyce hyperelasticity.
5.1 Corotational as Linearized Neo-Hookean
In Section 3.2, we described how our energy could be viewed
as a singularity-removing Taylor truncation of the original NeoHookean
energy ΨNeo. In order to establish a common basis for
comparison, we show next how linearizing ΨNeo yields corotational
elasticity to within a constant. Starting from the energy in
Equation (5),
ΨNeo = μ
2
(IC − 3) − μ log J +
λ
2
(log J)
2,
we can employ the linearization log J ≈tr (S − I) to establish that
μ
2
(IC − 3) − μ log J ≈ μ
2
||F − R||2
F ,
λ
2
(log J)
2 ≈ λ
2
tr2 (S − I) ,
and then conclude that (see proof in Appendix C)
ΨCR = μ
2
||F − R||2
F +
λ
2
tr2 (S − I) .
Note that ΨCR contains a μ/2 instead of the μ from Equation
(3), because the linearization introduces a factor of two that
should be accounted for using a Lamé reparameterization similar
to Section 3.4. The constant factor is otherwise irrelevant to our
discussion.
5.2 Comparison to Our Energy
From the discussions in Section 3.2 and Section 5.1, it is clear
that our model applies the approximation log J ≈ (J − 1), while the
corotational model applies the linearization log J ≈ tr (S − I). This
is a critical difference, because when the linearization tr (S − I) is
used, the actual volume J = σ0σ1σ2 is no longer computed. Instead,
it is assumed that (σ0σ1σ2 − 1) ≈ (σ0 + σ1 + σ2 − 3), which only
holds for small deformations where {σi} are near 1.
The unacceptable severity of the
linearization becomes clear when
we visualize its energy in the inset
on the right. The volume preservation
term λ
2 tr2 (S − I) is minimized
when tr Σ = 3. The corresponding
contour is shown in red and includes
inverted configurations, as
can be seen when the line strays
outside the upper right quadrant.
Thus, the energy can report that
volume has been perfectly preserved when in fact elements have
been crushed to zero volume or inverted. The red contour again resembles
Figure 1 from Stomakhin et al. (2012). In that paper, they
propose a “fixed” corotational model:
ΨFixCR = μ||F − R||2
F +
λ
2
(J − 1)
2. (49)
This model is halfway between our model and the original corotational
model. As with our model, the log J ≈ (J − 1) approximation
is applied to the λ
2 (log J)
2 term, but the problematic linearization
log J ≈ tr (S − I) is still applied to the μ log J term. While ΨFixCR is
more stable than ΨCR, we show in Figure 5 that it introduces visual
artifacts under large deformation.
Finally, we compare to the St. Venant-Kirchhoff (StVK) energy:
ΨStVK = μ||E||2
F +
λ
2
tr2 (E). (50)
Unlike the corotational energy, the linear (F − R) has been replaced
with a quadratic Green’s strain E = 1
2 (FT F − I). The true
volume J is still not computed anywhere, so all of the volume
linearization problems from the corotational model are present.
Irving et al. (2004) observed that this energy contains spuriously
stable inverted rest and flattened states, which makes the standard
StVK model unattractive for large deformation. Additional
compression resistance terms (Kikuuwe et al. 2009) have been
proposed:
ΨStVK,Kikuuwe = μ||E||2
F +
λ
2
tr2 (E) + κ
12  1 − J
6
3
, (51)
where the new term is clamped to zero when J ≥ 1. While this term
helps to preserve volume and remove the spurious stable configurations,
it deactivates during extension. The artifacts that arise
from the linearized volume term then persist in this regime. We
show some of these artifacts in Figure 10 and include the model in
a more extensive version of Figure 5 in the supplemental material.
5.3 Extension to Other Energies
Our eigenanalysis is sufficiently general that it can be used to analyze
any strain energy that employs the IC and J strain invariants
ACM Transactions on Graphics, Vol. 37, No. 2, Article 12. Publication date: March 2018.
Stable Neo-Hookean Flesh Simulation • 12:9
Fig. 5. A cylinder with 306,406 hexahedra and 965,004 degrees of freedom is stretched by a factor of 3.4 with ν = 0.49. The corotational model fails entirely,
and the “fixed” corotational model grows unnatural, sharp fins. Our stable Neo-Hookean model retains the rubber-like appearance of the filtered NeoHookean
model, exhibits slightly better volume preservation, and does not require any parameter tuning.
(see Section 7 for comments on IIC). We have found it straightforward
to extend to increased nonlinearities in IC.
General Model: Consider a strain energy density expressed in
terms of the IC and J invariants that leads to a system of the form:
Aˇ = γI I +γF ˇ
fˇ
f
T +γG gˇgˇT +γH Hˇ . (52)
The eigensystem can be obtained by setting η = γF /γH and ρ =
γG /γH in Equation (42), and μT = γI and λ(J − α) = γH in Equations
(45) through (48). We now examine two specific extensions
of our analysis.
Fung Hardening: The exponential hardening of Fung-like
models (Fung 2013; Pan et al. 2015; Wang and Yang 2016) is a secondary
feature of many biological tissues (Kautzman et al. 2012),
so we propose a stabilized Fung model:
ΨFung = μ0
2 (IC − 3) +
λ
2 (J − α)
2 + γ
2

e
μ1
2 (IC −3) − 1

, (53)
where the rightmost term is the novel one. Following the approach
of Section 3.2, the energy can be stabilized by setting α =
1 + μ0+γ μ1
λ . The PK1 then expands to
PFung(F) = μ0F + λ ∂J
∂F (J − α) +γ μ1e
μ1
2 (IC −3)
F. (54)
The unfolded, assembled Hessian can then be written as
Aˇ Fung =

μ0 +γ μ1e
μ1
2 (IC −3)

I + λgˇgˇT
+ λ(J − α)Hˇ + 2γ μ2
1e
μ1
2 (IC −3)ˇ
fˇ
f
T .
(55)
The following constants change in the eigenanalysis:
η = 2γ μ2
1 e
μ1
2 (IC −3)
λ(J−α ) ρ = 1
J−α μT =

μ0 +γ μ1e
μ1
2 (IC −3)

. (56)
Arruda-Boyce: The hardening term can also be represented with
a polynomial using the model of Arruda and Boyce (1993):
ΨAB =

n
i=1
ai (I
i
C − 3i
). (57)
We use the three-term expansion for illustrative purposes, shift the
constants to maintain consistency with the Neo-Hookean parameters,
and arrive at the following energy:
ΨAB = μ
2
(IC − 3) + μβ1
4 (I
2
C − 9) + μβ2
6 (I
3
C − 27) +
λ
2 (J − α)
2.
(58)
The energy is rest stable when α = 1 + μ
λ (1 + 3β1 + 9β2), leading
to a PK1 of
PAB(F) = μ


1 + β1IC + β2I
2
C

F + λ ∂J
∂F (J − α) (59)
and a Hessian of
Aˇ AB = μ


1 + β1IC + β2I
2
C

I + λgˇgˇT
+ λ(J − α)Hˇ + 2μ (β1 + 2β2IC ) ˇ
fˇ
f
T . (60)
This system leads to the following constants in the eigenanalysis:
η = 2μ (β1+2β2IC )
λ(J−α ) ρ = 1
J−α μT = μ


1 + β1IC + β2I
2
C

. (61)
6 RESULTS
We evaluate the robustness of our model across a variety of scenarios.
All quasi-static examples use a (potentially nonmanifold
(Mitchell et al. 2015)) hexahedral mesh with eight-point quadrature.
To show the agnosticism of our energy to the topology and
integration scheme, we also include a dynamic example using linear
tetrahedra. Quasi-static simulations are performed with a standard
Newton solver augmented with a line search, and linear systems
are solved using ViennaCL’s Conjugate Gradient implementation
running on OpenCL (Rupp et al. 2010). All Newton solves
were run until the absolute L2-norm of the force residual was less
than 1e−2. After normalizing by the number of degrees of freedom,
we are essentially solving to single precision. All simulations were
performed with a 2.3GHz Intel Xeon E5-2699 on eight cores, 118GB
of RAM, and an NVIDIA Quadro M6000 GPU.
Stretch Test. We stretch a cylinder composed of 306,406
hexahedra to 3.4 times its original length by translating hardconstrained
vertices on two opposing faces. We compare to three
models:
• The corotational energy ΨCR from Irving et al. (2004)
• The “fixed” corotational ΨFixCR of Stomakhin et al. (2012)
• The original Neo-Hookean energy ΨNeo, but “filtered” using
Teran et al. (2005).
In all cases, we call the original implementations in the PhysBAM
codebase (Dubey et al. 2011), which one of the original authors has
generously provided for us. The default force-filtering threshold
σf = 0.25 and Hessian projection threshold ϵp = −1e−4 from the
code were used. We show results for an extreme Poisson’s ratio of
ν = 0.49 in Figure 5. A wide range of additional ν, as well as comparisons
to theC2 model from Stomakhin et al. (2012) and the StVK
ACM Transactions on Graphics, Vol. 37, No. 2, Article 12. Publication date: March 2018.
12:10 • B. Smith et al.
Fig. 6. In the left graph, our model exhibits the best volume preservation under stretching. In the right graph, our CG iteration count is competitive with
the other models, and under large stretches, our model consistently converges faster. Our model took an average of 100.64s per frame, filtered Neo-Hookean
took 537.72s, and “fixed” corotational took 144.07s.
model from Kikuuwe et al. (2009), are shown in the supplemental
materials.
We first examine volume preservation. In Figure 6, the corotational
model loses over 83.0% of the original volume; linearizing
the volume penalty is clearly inadequate for this case. The
“fixed” corotational model only gains 4.5% volume, but the final
shape contains highly unnatural crinkling. An analogous artifact is
shown in Figure 6 of Stomakhin et al. (2012), and our analysis suggests
that it stems from the linearization hidden in the μ
2 ||F − R||2
F
term. Filtered Neo-Hookean gains 4.9% volume, slightly more than
the other model. Finally, our stable Neo-Hookean model gains 4.3%
volume, which is the smallest of any of the models.
In order to compare computational performance, we use the total
number of conjugate gradient (CG) iterations performed by the
Newton solver at each time step. This is more instructive than the
total Newton iterations in this case, because many of the solves in
Figure 5 completed in fewer than three iterations. It also removes
constant factors associated with low-level optimization that has
been applied to various components of the code, e.g., the PhysBAM
classes.
The “fixed” corotational model requires the most CG iterations,
roughly 2.3 times the total CG iterations of our model for the
stretch in Figure 5. This may be because at larger deformations,
the inconsistent approximations applied to the log J term begin to
conflict. On average, our stable Neo-Hookean model takes slightly
more CG iterations than filtered Neo-Hookean, with 1.2 times the
total CG iterations; the computational costs are virtually identical.
Our model contains none of the artifacts from the “fixed” corotational
model and preserves volume better than the filtered NeoHookean
model.
The CG iteration counts for the original corotational model are
not shown because they are an order of magnitude larger than the
other models. Converging to the degenerate geometry in Figure 5
is computationally expensive. We also attempted to simulate the
C2 continuous Neo-Hookean model from Stomakhin et al. (2012).
While the model converged fine for lower ν (see supplemental materials),
we were unable to locate settings that did not diverge before
completing the test for ν = 0.49.
Twist Test: We take a unit cube and subject it to large kinematic
rotations by hard-constraining one face of the cube and rotating it
twice by π/2 radians. Production quasi-static flesh simulations often
undergo rotations of this magnitude, so it is critical that models
remain robust under these conditions. The cube is discretized into
153 hexahedra, and ν = 0.49. Our model finds the correct twisted
shape, where the solves respectively took 12 and 16 Newton iterations,
and 3,278 and 4,956 CG iterations (Figure 9).
The “fixed” corotational model is also able to robustly resolve
both rotations in 11 and 15 Newton iterations and 2,940 and 4,348
CG iterations. The corotational model resolves the π/2 rotation
but then collapses. For the π rotation, 223 Newton iterations were
needed; resolving collapsed configurations is again costly. Collisions
were disabled, so the results involve purely the material
model. Filtered Neo-Hookean yields unusable results. The solve
stalls when the line search is enabled and diverges when it is
disabled.
Scramble. In Figure 3, we perform a “scramble” test similar to
those in Teran et al. (2005) and Stomakin et al. (2012) by randomly
placing the vertices of a unit cube within a cube of twice the rest
volume. To separate the behavior of our model from that of massregularization,
we perform a quasi-static solve. We pin four cube
corners in order to constrain the rigid modes. As shown in the
supplemental video, our model successfully recovers the rest pose.
Torso. In Figure 1, we perform a quasi-static flesh simulation
of arms and a torso that are driven kinematically by a set of internal
bones. The bones are transformed by the skinning of a production
rig, which often results in nonrigid transforms and intersecting
bones. Like the corotational model, our model is robust to
these noisy inputs. However, our model produces realistic bulges
and folds, while the corotational model produces artifacts and collapses
the shoulder entirely. Increasing ν only worsens the artifact;
the parameter accomplishes the opposite of its stated purpose. In
Figure 10, we compare our model to an StVK model with a compression
resistance term (Kikuuwe et al. 2009). The StVK model
fails to capture the fleshy profile of our model.
Collisions are enabled in this simulation and handled with a reference
map-based penalty force (Hirota et al. 2001; Irving et al.
ACM Transactions on Graphics, Vol. 37, No. 2, Article 12. Publication date: March 2018.
Stable Neo-Hookean Flesh Simulation • 12:11
Fig. 7. Top: Our model captures the bulge in the palm at the base of the fingers. Bottom: Our model captures the bulge in the webbing between the pointer
finger and the thumb. At ν = 0.2, corotational elasticity resembles contact-aware skinning. Increasing to ν = 0.488 only introduces artifacts and fails to
capture the features of our model. The simulation is driven by 24 kinematic bones in the hand and the forearm, and the mesh consists of 41,050 hexahedra
and 142,995 degrees of freedom. Stable Neo-Hookean took an average of 15.9 Newton iterations, 9,486.7 CG iterations, and 53.8 seconds. With corotational
at ν = 0.2, the average is 13.1 Newton iterations, 2,471.9 CG iterations, and 44.0 seconds. With corotational at ν = 0.488, the simulation did not complete.
©Disney/Pixar.
2004). Starting from a surface point sampling, a spatial grid acceleration
structure (McAdams et al. 2011) is used to locate sample
points in a tetrahedralization of the hexahedral mesh. If a sample
point collides with a tetrahedron, it is mapped to the reference
domain, where the closest surface point is computed. The surface
point and the original surface sample are mapped back to the deformed
domain and define a penalty spring. Due to the two linear
maps, this target position is not always the closest point on the deformed
surface mesh. If clean and conforming contacts are desired,
an additional local optimization is performed on the deformed surface
mesh that iteratively searches for closer points. This process
typically converges in a single iteration.
Hand. In Figures 7 and 8, we simulate the flesh on a hand
and forearm. As with the torso, bones drive the simulation,
but their transforms come from a skinning that is noisy and
contains difficult stretches and intersections. Similar to McAdams
et al. (2011), the simulation both cleans up skinning artifacts
around difficult joints and adds physics-based details that were
missing in the original rig.
The corotational model consistently collapses regions where
volume should be preserved, such as wrinkles in the palm of the
hand or folds in the webbing between the thumb and pointer
finger. As with the torso example, increasing ν can worsen artifacts
and in some cases causes the Newton solve to diverge.
Dynamics. In Figure 11, we show our material model undergoing
dynamic, extreme collisions. Despite the majority of the model
being squashed to almost zero volume, our elastic energy allows it
to recover. The simulation uses a tetrahedral mesh and shows that
our model is agnostic to element type and quadrature scheme.
7 CONCLUSIONS AND FUTURE WORK
We have designed a new hyperelastic energy that requires no
special machinery to remain stable under large deformations,
ACM Transactions on Graphics, Vol. 37, No. 2, Article 12. Publication date: March 2018.
12:12 • B. Smith et al.
Fig. 8. Top: The thumb at its original rest pose. Top, middle: When
pinched, the flesh between the thumb and pointer bulges with our model
set to ν = 0.488. Bottom, middle: With the corotational model, there is
no feature at ν = 0.2. Bottom: Only slight bulging appears in the corotational
model at ν = 0.488. The settings and timings are the same as Figure
7. ©Disney/Pixar.
including inversions. We are able to derive closed-form expressions
for the eigendecompositions of all its components. The sum
of these components also possesses a special structure that allows
us to write simple expressions for the eigensystem. Together,
these findings allow our energy to be used in Newton-based
implicit schemes that use conjugate gradient solvers.
We have established a hierarchy of Neo-Hookean approximations.
The corotational elasticity model is the simplest, as it applies
Fig. 9. The top face of the cube undergoes two successive π/2 rotations.
Stable Neo-Hookean and “fixed” corotational resolve the twist robustly.
Corotational has trouble with the second rotation, and filtered NeoHookean
is unable to resolve the shape. The cube is discretized with 153
hexahedral elements and the material models are simulated with ν = 0.49.
Fig. 10. Left: A bicep under an imposed elbow flexion simulated with our
stable Neo-Hookean model (ν = 0.488). An average time step took 13.7
Newton iterations, 5,860 conjugate gradient (CG) iterations, and 25.6 seconds.
Right: The same pose with StVK and a compression resistance term
(ν = 0.488, κ = 8 × 108). The StVK model fails to capture the fleshy silhouette
we observe with our model, and instead leaves a gap between the
bicep and forearm. An average time step with this model took 32 Newton
iterations, 14,446 CG iterations, and 140.3 seconds. ©Disney/Pixar.
the log J ≈ tr(S − I) linearization to both Neo-Hookean terms,
followed by the model of Stomakhin et al. (2012), which only applies
it to one. Our energy instead consistently uses log J ≈ (J − 1)
everywhere. This critical difference allows volume to be measured
accurately under large deformation, and improves both the
stability and quality of the simulation. The log J ≈ (J − 1) approximation
is a first-order truncation, so one direction for future
work is to determine the benefits of higher-order truncations.
In Section 4, we performed a complete eigenanalysis of the first
stress invariant IC and the volume measure J. Our analysis of J
also constitutes an analysis of the third stress invariant, IIIC = J 2.
These results can be used to analyze any energy that can be written
in terms of IC and IIIC, as we show with the Fung and ArrudaBoyce
models in Section 5.3. An eigenanalysis of the second invariant,
IIC = tr(CT C), is still missing, and could lead to a better
understanding of materials such as the Mooney-Rivlin and Verdona
models.
ACM Transactions on Graphics, Vol. 37, No. 2, Article 12. Publication date: March 2018.
Stable Neo-Hookean Flesh Simulation • 12:13
Fig. 11. We dynamically drop an octopus containing 394,268 tetrahedra and 278,991 degrees of freedom with ν = 0.46 onto three cylinders. We squash it
with another cylinder, but it recovers as soon as the cylinder is removed. For clarity, the cylinders are visualized with partial transparency. On the right is
a closeup of the most squashed state. The simulation took 10 semi-implicit (Baraff and Witkin 1998) substeps per frame to resolve the collisions, and each
frame took an average of 20.8 seconds. ©Disney/Pixar.
Finally, a better understanding of how Hessian indefiniteness
impacts Newton convergence is needed. Allowing a slight amount
of indefiniteness to enter into the global Hessian can sometimes
accelerate Newton convergence, so a better understanding of this
phenomenon could yield significant practical benefits.
APPENDIXES
A VOLUME HESSIAN EIGENVECTORS
The first three eigenvectors, ε0...2, are specified by Equations (39)
and (40). The remaining six, ε3...8, are of the form
Qk = 1
√
2
UDkVT. (62)
The only differences are in the Di matrices, which are
D3 =
⎡
⎢
⎢
⎢
⎢
⎢
⎣
0 00
0 01
0 −1 0
⎤
⎥
⎥
⎥
⎥
⎥
⎦
D6 =
⎡
⎢
⎢
⎢
⎢
⎢
⎣
000
001
010
⎤
⎥
⎥
⎥
⎥
⎥
⎦
D4 =
⎡
⎢
⎢
⎢
⎢
⎢
⎣
001
000
−100
⎤
⎥
⎥
⎥
⎥
⎥
⎦
D7 =
⎡
⎢
⎢
⎢
⎢
⎢
⎣
001
000
100
⎤
⎥
⎥
⎥
⎥
⎥
⎦
D5 =
⎡
⎢
⎢
⎢
⎢
⎢
⎣
010
−100
000
⎤
⎥
⎥
⎥
⎥
⎥
⎦
D8 =
⎡
⎢
⎢
⎢
⎢
⎢
⎣
010
100
000
⎤
⎥
⎥
⎥
⎥
⎥
⎦
.
B ORTHOGONALITY TO PSEUDO-CROSS-PRODUCTS
Theorem B.1. If B ∈ R3×3 can be diagonalized using the U and
V from F, then vec(B) = ˇ
b is orthogonal to all six pseudo-crossproduct
eigenvectors in Equation (62).
Proof. We denote the eigenvectors in flattened form as
vec(Qk ) = qˇk . If we can show that qˇT
k
ˇ
b = 0, then we are done.
First, we observe that qˇT
k
ˇ
b = tr(QT
k B). Next, we use the assumption
that B can be diagonalized: B = UΣBVT . Since Qk = UDkVT ,
we have tr(QT
k B) = tr(DT
k ΣB) = 0. Finally, we observe that Dk are
hollow matrices; i.e., they are zero along their diagonals. The trace
of a diagonal times a hollow matrix is zero. Thus, qˇT
k
ˇ
b ≡ tr(QT
k B) ≡
tr(DT
k ΣB) = 0.
Corollary B.2. The vector gˇ is orthogonal to all six pseudocross-product
eigenvectors.
Proof. The vector gˇ is generated from the matrix ∂J
∂F ,
i.e., vec( ∂J
∂F ) = gˇ. This matrix can be diagonalized as
∂J
∂F = U
⎡
⎢
⎢
⎢
⎢
⎢
⎣
σ1σ2 0 0
0 σ0σ2 0
0 0 σ0σ1
⎤
⎥
⎥
⎥
⎥
⎥
⎦
VT. (63)
By the preceding theorem, gˇ is orthogonal.
C COROTATIONAL AS LINEARIZED NEO-HOOKEAN
Lemma C.1. tr (S − I) is a linearization of log J.
Proof. Given the polar decomposition F = RS (Table 1), we have
J =det F=det S. By applying the identity det S=etr(log S)
, we obtain
log J =tr (log S). We then compute the Taylor expansion about S =
I,
log S = (S − I) − 1
2
(S − I)
2 +
1
3
(S − I)
3 ..., (64)
and employ its linear term to conclude: log J ≈tr (S − I).
Theorem C.2. μ/2||F−R||2
F is a linearization of μ/2(IC −3)−
μ log J.
Proof. This follows from an algebraic manipulation:
μ
2 (IC − 3) − μ log J = μ
2


F2
F − tr I
− μ log J
(applying Lemma C.1) ≈ μ
2


F2
F − tr I
− μtr (S − I)
= μ
2


F2
F + tr I − 2trS
= μ
2


F2
F + R2
F − 2tr(RT F)

= μ
2 ||F − R||2
F .
Therefore: μ
2 (IC − 3) − μ log J ≈ μ
2 ||F − R||2
F .
Corollary C.3. λ/2 tr2 (S − I) is a linearization of λ/2(log J)
2.
Proof. We can reuse Lemma C.1 and the linearization appears
immediately.
D EIGENSYSTEM WHEN ν ≥ 0.4
If fleshy material with ν ≥ 0.4 is being simulated, the regularized
origin barrier can be discarded, and the simpler eigensystem from
that results from Equation (13) can be used instead. The Hessian is
then
Aˇ = μI + λgˇgˇT + λ(J − α)Hˇ , (65)
ACM Transactions on Graphics, Vol. 37, No. 2, Article 12. Publication date: March 2018.
12:14 • B. Smith et al.
which is only a rank-one updated system. The deflation argument
from Section 4.6 still applies, and the rank-one updated cubic is
instead
ε¯
3 − ρ


gˇ



2
2
ε¯
2 − (1 + 2ρ J)||F||2
F ε¯ − (2 + 3ρ J)J = 0, (66)
where ρ = (J − α)
−1. The diagonal entries in the eigenvector Qk = 1
qk UDkVT can be written as
Dk =
⎡
⎢
⎢
⎢
⎢
⎢
⎣
α (σ0σ2 + σ1ε¯k ) 0 0
0 α (σ1σ2 + σ0ε¯k ) 0
0 0 ε¯
2
k − βσ2
2
⎤
⎥
⎥
⎥
⎥
⎥
⎦
α = 1 + ρ J
β = α + ρ(J + (σ2
0 + σ2
1 )ε¯i ).
(67)
Note that for the case of ρ = 0, Equation (40) is recovered. The final
eigenvalues are now the same as in Equations (45) to (48), except
that the roots of Equation (66) are used instead, and μT = μ.
ACKNOWLEDGMENTS
We would like to thank David Eberle and Audrey Wong for their
last-minute help getting the tetrahedral octopus simulation from
Figure 11 pushed through Presto and Fizt. We would also like to
thank Jonathan Page for providing both the bone geometry and
finger animation in Figure 7, Bill Wise for providing the geometry
and animation for the torso example in Figure 1, and Rick Sayre for
approving the use of the torso asset. Thanks also to Tony DeRose
and Doug James for their feedback on early drafts of this work.
REFERENCES
Ellen M. Arruda and Mary C. Boyce. 1993. A three-dimensional constitutive model
for the large stretch behavior of rubber elastic materials. Journal of the Mechanics
and Physics of Solids 41, 2 (1993), 389–412.
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). 43–54.
Adam W. Bargteil and Elaine Cohen. 2014. Animation of deformable bodies with quadratic
bézier finite elements. ACM Trans. Graph. 33, 3, Article 27 (June 2014), 10
pages.
Jan Bender, Matthias Müller, and Miles Macklin. 2015. Position-based simulation
methods in computer graphics. In Eurographics Tutorials. Eurographics
Association.
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. J. Biomechanics 38,
4 (2005), 657–665.
Javier Bonet and Richard D. Wood. 2008. Nonlinear Continuum Mechanics for Finite
Element Analysis. Cambridge University Press.
Sofien Bouaziz, Sebastian Martin, Tiantian Liu, Ladislav Kavan, and Mark Pauly.
2014. Projective dynamics: Fusing constraint projections for fast simulation. ACM
Trans. Graph. 33, 4, Article 154 (July 2014), 11 pages.
Allan F. Bower. 2009. Applied Mechanics of Solids. CRC Press.
James R. Bunch, Christopher P. Nielsen, and Danny C. Sorensen. 1978. Rank-one modification
of the symmetric eigenproblem. Numer. Math. 31, 1 (1978), 31–48.
Isaac Chao, Ulrich Pinkall, Patrick Sanan, and Peter Schröder. 2010. A simple geometric
model for elastic deformations. ACM Trans. Graph. 29, 4, Article 38 (2010),
6 pages.
Oscar Civit-Flores and Antonio Susín. 2014. Robust treatment of degenerate elements
in interactive corotational FEM simulations. Comput. Graph. Forum 33, 6 (2014),
298–309.
Pradeep Dubey, Pat Hanrahan, Ronald Fedkiw, Michael Lentine, and Craig Schroeder.
2011. PhysBAM: Physically based simulation. In ACM SIGGRAPH Courses. Article
10, 22 pages.
Yuan-cheng Fung. 2013. Biomechanics: Mechanical Properties of Living Tissues.
Springer Science & Business Media.
Joachim Georgii and Rüdiger Westermann. 2008. Corotated finite elements made fast
and stable. In Proceedings of the Fifth Workshop on Virtual Reality Interaction and
Physical Simulation (VRIPHYS’08). 11–19.
Gene H. Golub and Charles F. Van Loan. 2012. Matrix Computations. Vol. 3. JHU Press.
George Neville Greaves, A. L. Greer, R. S. Lakes, and T. Rouxel. 2011. Poisson’s ratio
and modern materials. Nat. Mater. 10, 11 (2011), 823–837.
Nicholas J. Higham. 2008. Functions of Matrices: Theory and Computation. SIAM.
Gentaro Hirota, Susan Fisher, A. State, Chris Lee, and Henry Fuchs. 2001. An implicit
finite element method for elastic solids in contact. In Proceedings of the Fourteenth
Conference on Computer Animation (Computer Animation’01). 136–254.
Geoffrey Irving, Joseph Teran, and Ronald Fedkiw. 2004. Invertible finite elements for
robust simulation of large deformation. In ACM SIGGRAPH/Eurographics Symposium
on Computer Animation (SCA’04). 131–140.
Michael A. Jenkins and Joseph F. Traub. 1970. A three-stage variable-shift iteration for
polynomial zeros and its relation to generalized rayleigh iteration. Numer. Math.
14, 3 (1970), 252–263.
Ryan Kautzman, Jiayi Chong, and Patrick Coleman. 2012. Stable, art-directable skin
and flesh using biphasic materials. In ACM SIGGRAPH Talks.
Ryo Kikuuwe, Hiroaki Tabuchi, and Motoji Yamamoto. 2009. An edge-based computationally
efficient formulation of Saint Venant-Kirchhoff tetrahedral finite elements.
ACM Trans. Graph. 28, 1, Article 8 (2009), 13 pages.
Allen Knutson and Terence Tao. 2001. Honeycombs and sums of hermitian matrices.
Notices Amer. Math. Soc 48, 2 (2001), 175–186.
Tamara G. Kolda and Brett W. Bader. 2009. Tensor decompositions and applications.
SIAM Rev. 51, 3 (2009), 455–500.
Tiantian Liu, Ming Gao, Lifeng Zhu, Eftychios Sifakis, and Ladislav Kavan. 2016.
Fast and robust inversion-free shape manipulation. Comput. Graph. Forum 35, 2
(2016), 1–11.
Sebastian Martin, Bernhard Thomaszewski, Eitan Grinspun, and Markus Gross. 2011.
Example-based elastic materials. ACM Trans. Graph. 30, 4, Article 72 (July 2011),
8 pages.
Aleka McAdams, Yongning Zhu, Andrew Selle, Mark Empey, Rasmus Tamstorf,
Joseph Teran, and Eftychios Sifakis. 2011. Efficient elasticity for character skinning
with contact and collisions. ACM Trans. Graph. 30, 4, Article 37 (July 2011),
12 pages.
Nathan Mitchell, Mridul Aanjaneya, Rajsekhar Setaluri, and Eftychios Sifakis. 2015.
Non-manifold level sets: A multivalued implicit surface representation with applications
to self-collision processing. ACM Trans. Graph. 34, 6, Article 247 (Oct.
2015), 9 pages.
Melvin Mooney. 1940. A theory of large elastic deformation. J. Appl. Phys. 11, 9 (1940),
582–592.
Matthias Müller, Julie Dorsey, Leonard McMillan, Robert Jagnow, and Barbara Cutler.
2002. Stable real-time deformations. In ACM SIGGRAPH/Eurographics Symposium
on Computer Animation (SCA’02). 49–54.
Matthias Müller, Jos Stam, Doug James, and Nils Thürey. 2008. Real time physics. In
ACM SIGGRAPH Classes. Article 88, 90 pages.
Rahul Narain, Matthew Overby, and George E. Brown. 2016. ADMM ⊇ projective
dynamics: Fast simulation of general constitutive models. In Proceedings of the
ACM SIGGRAPH/Eurographics Symposium on Computer Animation (SCA’16). 21–
28.
Matthieu Nesme, Paul G. Kry, Lenka Jeřábková, and François Faure. 2009. Preserving
topology and elasticity for embedded deformable models. ACM Trans. Graph. 28,
3, Article 52 (July 2009), 9 pages.
Raymond W. Ogden. 1997. Non-linear Elastic Deformations. Dover Publications.
Zherong Pan, Hujun Bao, and Jin Huang. 2015. Subspace dynamic simulation using
rotation-strain coordinates. ACM Trans. Graph. 34, 6, Article 242 (2015),
12 pages.
Ronald Rivlin. 1948. Large elastic deformations of isotropic materials. IV. Further developments
of the general theory. Philos. Trans. Roy. Soc. London A: Math. Phys.
Eng. Sci. 241, 835 (1948), 379–397.
Karl Rupp, Florian Rudolf, and Josef Weinbub. 2010. ViennaCL - A high level linear
algebra library for GPUs and multi-core CPUs. In International Workshop on GPUs
and Scientific Applications. 51–56.
Ruediger Schmedding and Matthias Teschner. 2008. Inversion handling for stable deformable
modeling. Vis. Comput. 24, 7 (2008), 625–633.
Christian Schüller, Ladislav Kavan, Daniele Panozzo, and Olga Sorkine-Hornung.
2013. Locally injective mappings. In Eurographics/ACMSIGGRAPH Symposium on
Geometry Proceedings (SGP’13). 125–135.
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
Courses. Article 20, 50 pages.
James G. Simmonds. 2012. A Brief on Tensor Analysis. Springer Science & Business
Media.
Alexey Stomakhin, Russell Howes, Craig Schroeder, and Joseph M. Teran. 2012. Energetically
consistent invertible elasticity. In ACM SIGGRAPH/Eurographics Symposium
on Computer Animation (SCA’12). 25–32.
Joseph Teran, Eftychios Sifakis, Geoffrey Irving, and Ronald Fedkiw. 2005. Robust
quasistatic finite elements and flesh simulation. In ACM SIGGRAPH/Eurographics
Symposium on Computer Animation (SCA’05). 181–190.
ACM Transactions on Graphics, Vol. 37, No. 2, Article 12. Publication date: March 2018.
Stable Neo-Hookean Flesh Simulation • 12:15
Matthias Teschner, Bruno Heidelberger, Matthias Muller, and Markus Gross. 2004. A
versatile and robust model for geometrically complex deformable solids. In Proceedings
of Computer Graphics International (CGI’04). 312–319.
Christopher D. Twigg and Zoran Kačić-Alesić. 2010. Point cloud glue: Constraining
simulations using the procrustes transform. In ACM SIGGRAPH/Eurographics
Symposium on Computer Animation (SCA’10). 45–54.
Nobuyuki Umetani, Ryan Schmidt, and Jos Stam. 2014. Position-based elastic rods. In
ACM SIGGRAPH/Eurographics Symposium on Computer Animation (SCA’14). 21–
30.
Huamin Wang and Yin Yang. 2016. Descent methods for elastic body simulation on
the GPU. ACM Trans. Graph. 35, 6, Article 212 (Nov. 2016), 10 pages.
Hongyi Xu, Funshing Sin, Yufeng Zhu, and Jernej Barbič. 2015. Nonlinear material
design using principal stretches. ACM Trans. Graph. 34, 4, Article 75 (July 2015),
11 pages.
Received September 2017; accepted November 2017
ACM Transactions on Graphics, Vol. 37, No. 2, Article 12. Publication date: March 2018.