Efficient Deformable Body Simulation using Stiffness

Efficient Deformable Body Simulation using Stiffness-Warped Nonlinear Finite Elements Wai-Hin Wayne Ngan∗ University of British Columbia

John E. Lloyd† University of British Columbia

Figure 1: A finite element sheet made from stiffness-warped linear tetrahedra collides with a sphere (first three frames). The same simulation can be achieved with greater fidelity, fewer elements and similar compute speed using nonlinear elements such as quadratic tetrahedra (last frame).

Abstract Finite element methods are commonly used in the graphics and virtual reality communities for simulating the behavior of deformable objects. For simplicity and low computational cost, practitioners often use a small-deformation linear elastic model in conjunction with tetrahedral elements, with a “stiffness warping” method sometimes employed to correct the distortions which can arise from large deformations. In this paper, we show how stiffness warping can be easily extended to more complex nonlinear elements, particularly hexahedra and quadratic tetrahedra. This involves identifying a “characteristic tetrahedron” within the element and using that to determine the warping transformation. We then show that nonlinear elements can provide advantages in both fidelity and/or computational speed (several fold for high accuracy models) as compared with linear tetrahedra. CR Categories: I.6.8 [Simulation and Modeling]: Types of Simulation—Animation, I.3.5 [Computer Graphics]: Computational Geometry and Object Modeling—Physically based modeling I.3.7 [Computer Graphics]: Three-Dimensional Graphics and Realism—Virtual reality Keywords: physically based animation, finite element method, simulation, deformation, stiffness warping

1

Introduction

Finite element methods are popular as a general purpose way to simulate the behavior of deformable bodies in computer graphics ∗ e-mail:

wayne927 at interchange.ubc.ca at cs.ubc.ca

† e-mail:lloyd

and virtual reality applications. For computational efficiency, tetrahedral elements are commonly used in combination with a strainstrain linear elastic model. Although such models can exhibit significant distortions in the presence of large deformations, this can be corrected with the recently introduced technique of “stiffness warping” [Müller et al. 2002; Müller and Gross 2004]. In this paper, we describe an easy way to extend stiffness warping to nonlinear elements, specifically hexahedrons and quadratic tetrahedrons (henceforth known as quadratic tets), and study the resulting speed/performance trade-offs. Our primary contributions are to ⋄ Describe how stiffness warping can be easily and effectively applied to nonlinear elements; ⋄ Show that with this, nonlinear elements can offer significant speed/performance advantages. For high accuracy models, we have found that nonlinear elements can improve computational efficiency by several times. For lower resolution models, nonlinear elements can offer better fidelity for the same computational speed. Nonlinear elements are also quite easy to implement, with the only added complexity being the precomputation of the stiffness matrices, which requires numerical integration (Appendix A). Moreover, nonlinear elements can sometimes offer fidelity that simply cannot be achieved with tetrahedra (as when simulating incompressible materials, where tetrahedra are prone to a numeric stiffness known as locking [Hughes 2000]). A note of clarification: nonlinear in this paper refers to the element interpolating functions. It does not refer to the material constitutive model, which is still small-strain linear with a warping correction.

2 Related Work Efficient simulation of deformable bodies has received much attention in the graphics and virtual reality communities. Mass spring systems, introduced by [Terzopoulos et al. 1987], have been popular because of their simplicity, and are widely used in cloth simulation [Bridson et al. 2003; Choi and Ko 2002], but can be difficult to tune for physically accurate behavior (though there have been attempts to address this [Hong et al. 2006]). Our experience has been that the limiting computational factor is generally the matrix solve

associated with implicit integration (equation 9 below), which depends more on node count and connectivity than whether springs or finite elements are being used. Linear small-strain finite element methods can be very efficient, especially if it is possible to fully or partially pre-invert the associated matrices; such techniques have been popular for surgical simulation [Goksel et al. 2005; Cotin et al. 1999]. However, they can exhibit significant errors in the case of large deformations. This problem can be addressed by “stiffness warping” methods, described below, in which the elemental rotations are separated from the strain calculations [Müller et al. 2002; Müller and Gross 2004; Nesme et al. 2005; Irving et al. 2006]. Fast simulation can also be achieved by coordinate reduction techniques in which the object’s deformability is expressed by a limited number of basis functions. The simplest of these approaches is linear modal reduction [Pentland and Williams 1989], which is very fast but still suffers from small-strain distortional errors; this has been addressed by extensions to large displacement Green strain [Barbiˇc and James 2005] or by incorporation of stiffness warping [Choi and Ko 2005]. A reduced deformation basis can also be produced from boundary element methods [James and Pai 2003]. For biomedical applications, there has been considerable effort directed at producing realistic models of deformable tissues, including incompressibility, using either projection techniques [Irving et al. 2007] or mixed-formulation methods involving higher order elements [Roth et al. 1998]. Finally, the computational advantages of quadratic tets have also been explored in [Mezger and Straßer 2006], though using Green strain instead of a stiffness warping approach.

3

01 − ν BB ν λB B ν0 E= ν B B 0 0

Small-Strain Linear Elasticity

Elastic body deformations are generally described in terms of a mapping x = φ(X) from undeformed material coordinates X = (X, Y, Z) to deformed spatial coordinates x = (x, y, z). The derivative of this map with respect to X is known as the deformation gradient F: F = ∂x/∂X. Coordinate displacements u = (ux , uy , uz ) ≡ x − X describe the deviation of points from their rest position. The displacement gradient ∂u/∂X satisfies ∂u/∂X = ∂x/∂X − I = F − I and is used to define the Cauchy strain tensor ε: ε=

1 2

∂u ∂X

+

∂u ∂X

T !

=

1 F + FT − I. 2

(1)

This is a symmetric quantity with six independent components, which can be arranged as a vector ε¯:

∂uy ∂ux ∂uz ∂uy ∂uz ∂ux ∂uy ∂uz ∂ux ε¯ = , , , + , + , + ∂x ∂y ∂z ∂y ∂x ∂z ∂x ∂z ∂y (2) For a linear elastic material, strain is linearly related to the stress σ, which can be expressed in matrix form as

λ≡

3.2

Accepted for poster presentation at I3D 2008

(3)

1 2

0 0 0 −ν 0 0

1 2

0 0 0 0 −ν 0

1 2

0 0 0 0 0 −ν

1 C C C C C C A

Eν . (1 + ν)(1 − 2ν)

Discretization using Finite Elements

Finite element modeling involves decomposing a body into a set of discrete volumetric elements, each associated with a set of boundary nodes and shape functions that interpolate a quantity within the element from the quantity’s nodal values. In particular, x = φ(X) can be approximated as x=

n X

xi Ni (X)

i=1

where Ni are the shape functions (one per node), xi are the deformed coordinate values at node i, and n is the number of nodes. The displacement mapping becomes n X

n X

i=1

i=1

(xi − Xi )Ni (X) =

ui Ni (X)

(4)

with ui giving the displacement at node i. The displacement gradient is then n X ∂Ni (X) ∂u = ui ∂X ∂X i=1 Instead of X, shape functions are often expressed in terms of natural coordinates S which are more appropriate to the element’s geometry. In that case, X ∂Ni (S) ∂S ∂u ui = . ∂X ∂S ∂X i=1 n

∂S/∂X can be determined from the inverse of the element Jacobian J, which is defined as follows: ∂X = J≡ ∂S

n X

∂Ni (S) Xi ∂S i=1

!

(5)

If the number of coordinates in S exceeds three (as with quadratic tets, described below), then some coordinates are redundant and the Jacobian J must be augmented with this redundancy information in order to make it square and invertible. Arranging the nodal displacements ui into a single vector ue , (4) can be expressed as u(X) = H ue where

0N (X) H≡ 0 1

σ ¯ = E ε¯.

ν ν 1−ν 0 0 0

where λ is the first Lame coefficient, defined as

Background

ν 1−ν ν 0 0 0

u=x−X=

We give a brief overview of finite element modeling with smallstrain linear elasticity. More details are available in various books on the subject [Bonet and Wood 1997; Hughes 2000].

3.1

For isotropic materials, E is characterized by only two parameters, Young’s modulus E and Poisson’s ratio ν, and takes the form

0

0 N1 (X) 0

0 0 N1 (X)

... ... ...

Nn (X) 0 0

0 Nn (X) 0

0 0 Nn (X)

1 A.

(-1, 1, 1) 8

(-1, -1, 1) 5

(1, -1, 1)

6

(1, 1, 1) 7

(-1, 1, -1)

3

4 1 (-1, -1, -1)

(1, 1, -1)

2 (1, -1, -1)

Figure 3: A hexahedral element showing the natural coordinates (S1 , S2 , S3 ) of its eight nodes. Figure 2: Small-strain distortion of an elastic sheet modeled using quadratic tets (left), corrected using the warping described in this paper (right).

ε¯(X) = B ue ,

∂ ∂x

∂ ∂y

0 ∂ ∂x

0

∂z

∂ ∂z

0

The technique of “stiffness warping” [Müller and Gross 2004] corrects for this problem in the case of tetrahedral elements by removing elemental rotation from the deformation gradient F. This is done by computing a polar decomposition such that

1

0

B0 B B0 B≡B B ∂ B ∂y ∂

0 0C C

∂ C ∂z C

H. 0C C

(6) F = RP

∂ A ∂x ∂ ∂y

From the principle of virtual work, it can be shown that BT maps the stress components σ ¯ back to nodal forces, and so B can be used in conjunction with (3) to determine the element stiffness matrix Ke relating nodal displacements to the nodal forces fe : Z

fe = Ke ue ,

BT EB dV.

Ke =

(7)

Ke is formed by taking a volume integral over the element. For nonlinear elements, this must be done numerically (Appendix A).

Dynamical Simulation

A second order ODE giving the dynamics of the entire system can now be formed. Letting x, x0 , and u ≡ x − x0 represent the positions, rest positions, and deformations of all nodes, we have M¨ x + Cx˙ + Kx = fx − f0

To preserve stability with reasonable time steps, the ODE (8) is usually solved using an implicit integrator. Letting xi and x˙ i denote the positions and velocities at step i, and h the time step size, then a first-order implicit Euler scheme involves solving for x˙ i+1 from

M + hC + h2 K x˙ i+1 = Mx˙ i − h Kxi + f0 − fx . (9) The left side matrix will be symmetric positive definite (SPD) as long as C is symmetric positive semi-definite (SPSD) (which is true for Raleigh damping), permitting efficient solution methods such as Cholesky factorization or conjugate gradient iteration.


fe = Re Ke RTe x − Re Ke x0 .

(11)

This is then used to reassemble and replace K and f0 in (8) with K′ and f0′ . The overhead of doing this is not prohibitive but it does mean that K cannot be preinverted.

4 Computing Element Stiffness Matrices 4.1

Tetrahedrons

(8)

where K is the global stiffness matrix formed by assembling the element stiffness matrices Ke of (7), M and C are the mass and damping matrices, fx are the external forces, and f0 ≡ −Kx0 . In this paper, we use Raleigh damping, so that C = αM + βK for scalars α and β.

(10)

where R is an orthogonal rotation matrix and P is symmetric. Strains are then computed from P − I instead of (1). To maintain consistent behavior in the presence of element inversions, it is important to ensure that det(R) = 1 [Irving et al. 2006]. Given R, one can transform the current node coordinates into an unrotated coordinate system, compute the nodal forces, and rotate back. This involves forming a block diagonal matrix Re of R entries, and replacing (7) with

V

3.3

Stiffness Warping

Small-strain linear models exhibit significant errors whenever an element rotates significantly from its initial position (Figure 2).

Then from (2) we can calculate the strain ε¯: 0

3.4

The shape functions for the four node tetrahedron are linear and can be expressed in material coordinates as Ni = 1

X

Y

Z ai

where ai is a vector of coefficients. If the coefficient vectors are arranged as the columns of a matrix A, then we can determine A from the nodal rest positions Xi by noting that 0

1 B1 A= 1 1

X1 X2 X3 X4

Y1 Y2 Y3 Y4

1−1

Z1 Z2 C Z3 A Z4

.

The last three elements of each row of A equal ∂Ni /∂X, from which we can form the B matrix and then compute Ke from (7). Since B is constant, this can be done by multiplying BT E B by the elements rest volume.

(1, 0, 0, 0) 4 (0.5, 0.5, 0, 0) 8

(0.5, 0, 0, 0.5) 10

(0.5, 0, 0.5, 0) 9 7 (0, 0.5, 0, 0.5)

1 (0, 1, 0, 0)

3 6

(0, 0, 0, 1)

(0, 0, 0.5, 0.5)

5 (0, 0.5, 0.5, 0) 2 (0, 0, 1, 0)

Figure 4: A quadratic tetrahedral element showing the natural coordinates (S1 , S2 , S3 , S4 ) of its 10 nodes.

4.2

Figure 5: The characteristic tetrahedron (bold lines) for a quadratic tetrahedron (thin lines) is formed by connecting the four vertex nodes.

Hexahedrons

Hexahedrons can be defined using natural coordinates S ≡ (S1 , S2 , S3 ) formed from the triple Cartesian product of the interval [−1, 1]. In this space, it consists of a cube with vertices at (1, 1, 1), (1, 1, −1), (1, −1, 1), etc. (Figure 3). Expressed in natural coordinates, the shape functions are Ni =

1 8

3 Y

X1 . . . X4 and x1 . . . x4 . Adopting the first node as a reference coordinate and then forming the matrices Dm = X2 − X1 D = x2 − x1

(1 + sgn(Sik )Sk )

(12)

where Sik denotes the value of Sk at node i. From these, we can calculate ∂Ni /∂S, convert this to ∂Ni /∂X using ∂Ni ∂S ∂Ni −1 ∂Ni = = J ∂X ∂S ∂X ∂S and then use (6) to compute B for any value of S. Since the shape functions are nonlinear with respect to X, B is not constant and so the formation of Ke from (7) involves numerical integration as described in Appendix A.1.

Quadratic Tetrahedrons

Quadratic tets can be defined using natural coordinates S ≡ (S1 , S2 , S3 , S4 ) which correspond to the barycentric coordinates of a unit tetrahedron (Figure 4). The shape functions are N1 N2 N3 N4

= S1 (2S1 − 1) = S2 (2S2 − 1) = S3 (2S3 − 1) = S4 (2S4 − 1)

N5 = 4S1 S2 N6 = 4S2 S3 N7 = 4S3 S1

N8 = 4S1 S4 N9 = 4S2 S4 N10 = 4S3 S4

(13)

The natural coordinates are redundant, satisfying the identity S1 + S2 + S2 + S4 = 1. To maintain invertibility, the element Jacobian J (Section 3.2) is augmented to contain this identity:

J(S) =

∂X/∂S 1 1 1 1

(14)

Using this Jacobian, we can evaluate B for any value of S using the same manner as for hexahedrons. Again, the shape functions are nonlinear with respect to X and so Ke must be formed by numerical integration as described in Appendix A.2.

5

x3 − x1

X4 − X1

x4 − x1

it is fairly easy to show ([Irving et al. 2006]) that x = D D−1 m (X − X1 ) + x1

k=1

4.3

X3 − X1

Stiffness Warping Implementation

For tetrahedral elements, F can be computed directly from the node positions. Let the rest and deformed node positions be given by


and hence

∂x = F = DD−1 m , ∂X from which R can be determined as per (10).

(15)

For nonlinear elements, F varies within the element so any value used for determining R must be determined using some representative value. In [Irving et al. 2006], this is interpolated from shape functions, while [Wicke et al. 2007] use the general purpose registration techniques of [Horn et al. 1988]. Alternatively, the approach we present in this paper is to identify a characteristic tetrahedral structure within the element, use it in conjunction with (15) to compute an approximation to F, and then from this extract R. It should be noted that this approximation to F does not affect the computation of Ke , which is precomputed using the element’s shape functions as described in the previous section. For quadratic tets, we form the characteristic tetrahedron from the four vertex nodes, as shown in Figure 5. Given these, we can compute D (along with a precomputed D−1 m ) and hence determine R exactly as for the tetrahedral case. For a hexahedron, we first tessellate it into five tetrahedra using a Freudenthal cut (Figure 6). We then take the innermost tetrahedron to be the characteristic tetrahedron, as this is connected to all the hexahedron’s faces and has the greatest symmetry. Since there are two possible tessellations associated with a Freudenthal cut, there are two possible inner tetrahedron to choose from. In this paper, the selection was made arbitrarily. Alternatively, one could choose both inner tetrahedra and average the two resulting R values, thereby obtaining a value that is symmetric with respect to all vertices.

6 Performance Experiments To evaluate our warping method, we implemented a number of test cases and compared their performance using the tetrahedral, hexahedral, and quadratic elements described above. Those described here include: hollow tubes hanging from a mobile horizontal platform (Figure 7), a thick sheet dropped onto a rigid sphere (teaser

corner nodes. This was repeated for all three element types at different resolutions, with the results shown in Figure 9. Tet A

Tet B

Tet D

Tet C

Tet E

Again, tetrahedral elements converge slowly, approaching a limit around (17280 elements, 4375 nodes). Hexahedral elements converge faster, to a limit near (1024 elements, 1445 nodes), and quadratic elements converge quickly, to a limit near (360 elements, 773 nodes). The limit behavior is also a bit different for the tetrahedral model, suggesting a limitation in modeling capabilities of the tetrahedra. Computation times (per 10 msec step) were 935, 280, and 70 msec, respectively, and so here quadratic elements were more than ten times as efficient as tetrahedral ones.

Figure 6: Freudenthal tessellation of a hexahedron. The innermost tetrahedron, E, is used as the characteristic tetrahedron.

0.5 (#elements, #nodes) 240, 72 800, 240 1400, 420 2160, 648 2400, 720 3080, 924 3840, 1152

0.4

x position

0.3 0.2 0.1 0 -0.1 -0.2

(a) Tetrahedrons

-0.3 0.5

1

1.5

2

2.5

3

3.5

4

time (s) 0.5

Figure 7: Vertically suspended tubes (left), and user interaction with a tapered ellipsoidal object (right).

(#elements, #nodes) 48, 72 160, 240 280, 420 432, 648 616, 924

0.4

x position

0.3

figures), and a tapered ellipsoidal object manipulated interactively (Figure 7). Our tests were implemented using Java 1.5 and the ArtiSynth biomechanical modeling platform (www.artisynth.org). The computer used was a 2.8 GHz Intel Pentium IV processor running Linux. To solve the matrices of (9), we used the Pardiso direct sparse solver [Schenk and Gärtner 2006].

0.2 0.1 0 -0.1 -0.2

(b) Hexahedrons

-0.3 0.5

1

1.5

2

2.5

3

3.5

4

time (s) 0.5

6.1

Dynamic fidelity

0.3

Tetrahedral elements converge the slowest, approaching a limit behavior around (2160 elements, 648 nodes). Hexahedral and quadratic elements both converge quickly, at (160 elements, 240 nodes) and (100 elements, 220 nodes), respectively. Computation times (per 10 msec step) were 51, 17, and 12 msec, respectively, indicating that quadratic elements are four times as efficient as tetrahedral ones. The experiment for the sheet consisted dropping the sheet onto a rigid sphere and monitoring the z coordinate of one of the bottom


x position

For the tube and sheet examples, we analyzed the dynamic behavior at different element resolutions to determine which element choice is the most computationally efficient near the convergence limit (the point at which increased resolution does not affect the simulation and hence indicates ideal behavior). All simulations were performed using (9) with a step size h of 10 msec., and with elasticity parameters of E = 500000 and ν = 0.33. The experiment for the hanging tube consisted of a half-sawtooth motion in which the suspending platform was moved first forward and then backward along the horizontal x axis, while monitoring the x coordinate of an outer node at the tube’s bottom. This was repeated for tetrahedral, hexahedral, and quadratic elements, at a variety of resolutions, with the results shown in Figure 8.

(#elements, #nodes) 30, 66 100, 220 175, 385 270, 594

0.4

0.2 0.1 0 -0.1 -0.2

(c) Quadratic tetrahedrons

-0.3 0

0.5

1

1.5

2

2.5

3

3.5

4

time (s)

Figure 8: Convergence results for the suspended tube subject to a half-sawtooth input.

6.2

Visual Evaluation

Visual demonstrations of low resolution, interactive versions of the tube, sheet, and tapered ellipsoid examples are presented in the accompanying video. Resolutions were adjusted to equalize the perstep simulation times, allowing evaluation of the quality obtained for equal computational effort. In the tube example, the tetrahedral model exhibits a higher frequency (also seen in Figure 8), indicating a stiffer response (a known characteristic of tetrahedral modes). In

10000

-4.5

z position

-5

-5.5

time per step (ms)

(#elements, #nodes) 160, 75 960, 324 5120, 1445 10000, 2646 17280, 4375 23660, 5832

-6

-6.5

(a) Tetrahedrons

1000

100

10

tets

-7 1

2

3

4

5

quads

time (s)

1

-4.5

100 (#elements, #nodes) 32, 75 192, 324 576, 845 1024, 1445 1620, 2166 2000, 2646

-5

z position

hexs

6

-5.5

1000

10000

number of nodes

Figure 10: Computation times vs. the number of nodes, in the sheet example, for tetrahedral, hexahedral, and quadratic elements.

-6

7 Conclusion -6.5

(b) Hexahedrons -7 1

2

3

4

5

6

time (s) -4.5 (#elements, #nodes) 20, 71 80, 227 360, 773 640, 1317 810, 1643

z position

-5

-5.5

We have shown that the popular stiffness-warping technique for tetrahedral deformable body simulation can be extended to nonlinear elements by identifying a characteristic tetrahedron for the element from which the warping rotation can be extracted. Experimental tests and convergence analysis using hexahedrons and quadratic tets indicate that this can yield greater simulation fidelity for equivalent computational cost, or a several-fold gain in computational speed for high-fidelity models. Quadratic tets would seem to hold the largest advantage, since they appear to be roughly as efficient to compute as hexahedra, and are easier to fit to arbitrary objects. The can be created from an tetrahedral mesh by simply adding nodes to the edge midpoints. If a higher resolution surface mesh is available, one can fit the elements to the surface by projecting the outer edge nodes onto the mesh.

-6

-6.5

(c) Quadratic tetrahedrons -7 1

2

3

4

5

6

time (s)

Figure 9: Convergence results for the sheet colliding with the sphere.

Our aim has been to show that applying stiffness-warped smallstrain linear elastic models to nonlinear elements is in fact quite straightforward and can offer noticeable performance advantages.

Acknowledgments all examples, particularly that of the tapered ellipsoid, the nonlinear models exhibit a somewhat less stiff and more fluid motion. The stiffness warping approach appears to perform well, with no noticeable distortions in any of the test examples, indicating the practicality of the characteristic tetrahedron approach. For reference, behavior without warping is illustrated in Figure 2.

6.3

Computation time

Figure 10 shows the per-step computation times for different element types as a function of the number of nodes. The hexahedral elements are the costliest, because their nodes exhibit greater connectivity, leading to a denser solve matrix in (9). However, this is probably not a general result. If we consider a partitioning of space based on a rectilinear grid and Freudenthal cuts within the grid cells, then, ignoring boundaries, the nodal connectivity for tetrahedra, hexahedra, and quadratic tets is, respectively, 12, 26, and 30. This suggests that the per-node computational cost of quadratic tets will sometimes be slightly larger than for hexahedra.


This work was supported by the Natural Sciences and Engineering Research Council of Canada (NSERC), through a Strategic Grant. We would also like to thank Ian Stavness for valuable help.

References BARBI Cˇ , J., AND JAMES , D. L. 2005. Real-time subspace integration for St. Venant-Kirchhoff deformable models. In SIGGRAPH ’05: ACM SIGGRAPH 2005 Papers, ACM Press, New York, NY, USA, 982–990. B ONET, J., AND W OOD , R. D. 1997. Nonlinear Continuum Mechanics for Finite Element Analysis. Cambridge University Press. B RIDSON , R., M ARINO , S., AND F EDKIW, R. 2003. Simulation of clothing with folds and wrinkles. In SCA ’03: Proceedings of the 2003 ACM SIGGRAPH/Eurographics symposium on Computer animation, Eurographics Association, Aire-la-Ville, Switzerland, Switzerland, 28–36. C HOI , K.-J., AND KO , H.-S. 2002. Stable but responsive cloth. In SIGGRAPH ’02: Proceedings of the 29th annual conference

on Computer graphics and interactive techniques, ACM Press, New York, NY, USA, 604–611. C HOI , M. G., AND KO , H.-S. 2005. Modal warping: Realtime simulation of large rotational deformation and manipulation. IEEE Transactions on Visualization and Computer Graphics 11, 1, 91–101. C OTIN , S., D ELINGETTE , H., AND AYACHE , N. 1999. Real-time elastic deformations of soft tissues for surgery simulation. IEEE Transactions on Visualization and Computer Graphics 5, 1, 62– 73. G OKSEL , O., S ALCUDEAN , S. E., D I M AIO , S. P., ROHLING , R., AND M ORRIS , W. J. 2005. 3d needle-tissue interaction simulation for prostate brachytherapy. In MICCAI, 827–834. H ONG , M., J UNG , S., C HOI , M.-H., AND W ELCH , S. W. J. 2006. Fast volume preservation for a mass-spring system. IEEE Comput. Graph. Appl. 26, 5, 83–91.

T ERZOPOULOS , D., P LATT, J., BARR , A., AND F LEISCHER , K. 1987. Elastically deformable models. In SIGGRAPH ’87: Proceedings of the 14th annual conference on Computer graphics and interactive techniques, ACM Press, New York, NY, USA, 205–214. W ICKE , M., B OTSCH , M., AND G ROSS , M. 2007. A finite element method on convex polyhedra. Computer Graphics Forum (Eurographics) 26, 3, 355–364.

A Computing Stiffness Matrices using Gaussian Quadrature To compute the integral in equation (7), we use Gaussian quadrature. This involves replacing the integral with a weighted sum of integrand values evaluated at selected points, with the points and weights depending on the domain and the accuracy desired. For example, a one-dimensional 2-point quadrature on the interval [−1, 1] approximates an integral by

H ORN , B. K. P., H ILDEN , H. M., AND N EGAHDARIPOUR , S. 1988. Closed-form solution of absolute orientation using orthonormal matrices. J. Opt. Soc. Am. A 5, 7, 1127. H UGHES , T. J. R. 2000. The Finite Element Method: Linear Static and Dynamic Finite Element Analysis. Dover Publications. I RVING , G., T ERAN , J., AND F EDKIW, R. 2006. Tetrahedral and hexahedral invertible finite elements. Graph. Models 68, 2, 66– 89. I RVING , G., S CHROEDER , C., AND F EDKIW, R. 2007. Volume conserving finite element simulations of deformable models. In SIGGRAPH ’07: ACM SIGGRAPH 2007 papers, ACM Press, New York, NY, USA, 13.

Z

−1

¨ M ULLER , M., AND G ROSS , M. 2004. Interactive virtual materials. In GI ’04: Proceedings of Graphics Interface 2004, Canadian Human-Computer Communications Society, School of Computer Science, University of Waterloo, Waterloo, Ontario, Canada, 239–246. ¨ M ULLER , M., D ORSEY, J., M C M ILLAN , L., JAGNOW, R., AND C UTLER , B. 2002. Stable real-time deformations. In SCA ’02: Proceedings of the 2002 ACM SIGGRAPH/Eurographics symposium on Computer animation, ACM Press, New York, NY, USA, 49–54. N ESME , M., PAYAN , Y., AND FAURE , F. 2005. Efficient, physically plausible finite elements. In Eurographics (short papers), J. Dingliana and F. Ganovelli, Eds. P ENTLAND , A., AND W ILLIAMS , J. 1989. Good vibrations: model dynamics for graphics and animation. In SIGGRAPH, ACM, J. J. Thomas, Ed., 215–222. ROTH , S. H. M., G ROSS , M. H., T URELLO , S., AND C ARLS , F. R. 1998. A Bernstein-Bezier based approach to soft tissue simulation. Computational Graphics Forum (Eurographics) 17, 3, 285–294. ¨ S CHENK , O., AND G ARTNER , K. 2006. On fast factorization pivoting methods for sparse symmetric indefinite systems. Elec. Trans. Numer. Anal 23, 158–179.


f (x)dx ≈

2 X

wi f (xi ) = w1 f (x1 ) + w2 f (x2 )

(16)

i=1

with w1 = w2 = 1 and x1 , x2 = ±

p

1/3.

Our volume integrals are done using natural coordinates S, which makes it easier to select quadrature points and weights.

A.1

Hexahedrons

Applying a P point quadrature to each natural coordinate interval, we obtain a product quadrature which approximates (7) with

JAMES , D. L., AND PAI , D. K. 2003. Multiresolution Green’s function methods for interactive simulation of large-scale elastostatic objects. ACM Trans. Graph. 22, 1, 47–82. M EZGER , J., AND S TRASSER , W. 2006. Interactive soft object simulation with quadratic finite elements. In AMDO, 434–443.

1

Ke =

P X P P X X

wi wj wk BTijk E Bijk Jijk

i=1 j=1 k=1

where J is the determinant of the element Jacobian J (section 3.2) and relates volumes in natural coordinates S to volumes in material coordinates X. The subscripts ijk indicate values for B and J evaluated at the i-th value of S1 , the j-th value of S2 , and the k-th value of S3 . As is common for hexahedral elements, we used a simple 2-point quadrature (P = 2), as shown in (16), p with weights w1 = w2 = 1 and each coordinate of S sampled at ± 1/3.

A.2

Quadratic Tetrahedrons

The integral (7) is discretized into Ke =

4 1X wk BTk E Bk Jk 6 k=1

where J is the determinant of the augmented element Jacobian J (equation 14), and the subscript k denotes values for B and J evaluated at Sk , where = {a, b, b, b} = {b, a, b, b} = {b, b, a, b} = {b, b, b, a} √ √ for a = (5 + 3 5)/20 and b = (5 − 5)/20. The weights in the sum are all wk = 1/4. S1 S2 S3 S4