1 Introduction
Intricate physical and mathematical models for simulating soft or rigid bodies are the centre of intensive research efforts. Computer animation, laparoscopic haptic surgery simulations, graphical special effects or robotic manipulation are just a few of the fields where such models play a key role. The focus of our research is to address the issue of carrying out stable, physicsbased simulations at interactive update rates. For this purpose, we have conceived an elementary framework for building deformable objects with soft or stiff constraints. Using this framework, we can test the efficiency and impact of several, wellknown numerical integrators. To address our goal, we look find an equilibrium between an integrator’s computational overhead, its precision, symplecticity and, most importantly, its inherent stability.
Following the aforementioned goals, we present a short literature survey on the numerical simulation of soft or semirigid bodies in section 2. The mathematical apparatus for creating and simulating constrained objects is explained in section 3, and our own, deformable linear object model case study is presented in section 4. We present a general update strategy that supports any type of explicit or semiimplicit integration method in section 5. Finally, the efficiency of several popular integrators is analysed and discussed in section 6, and we conclude this research in section 7.
2 Related work
Stability as a central attribute to Physicsbased object simulation within the Computer Graphics world was addressed by directly consider inherently stable integrators. The most popular method employed is the first order, implicit Euler scheme. Servin et al. [SLM06] treated infinitely stiff springs as kinematic constraints and developed a system capable of simulating elastic behaviour at the cost of solving a large sparse linear system for each iteration. The hybrid method of Schroeder et al. [SKZF11] uses explicit updates for the elastic forces and implicit strategies for other components. This complex idea tries to reduce the effect of suppressing material vibrations due to the use of pure implicit integrators. By alternating between implicit, semiimplicit and explicit methods, Volino et al. [VMT05] demonstrated how different cloth materials can be modeled to overcome the inabilities of one integrator to support a certain material property. Some object models allow using for computing an analytical force approximation for an implicit solver as demonstrated by Mesit et al. [MGC07]. Their method of simulating gas filled soft bodies can thus easily support any integrator, but it relies on the topological structure for all analytical derivations. Finally, structured massspring models using pure implicit integrators were popularized by Baraf and Witkin [BW98] in their seminal paper on stable cloth simulation and by Desbrun et al. [DSB99] where inverse dynamics were used to tackle outstretching artefacts. For a discussion of the performance of some popular implicit solvers we direct the user to the paper of Hauth et al. [HE01].
Explicit methods are more popular due to their reduced complexity. Finite element simulations using explicit integrators were performed by Fierz et al. [FSH10]. To stabilize the numerical process, the authors proposed modifying the stiffness matrices of illconditioned tetrahedral elements. This allowed their simulation to use higher time steps in spite of the explicit updates. Although less computationally demanding and stable than their implicit counterparts, explicit integrators can share a fair amount of stability, energy conserving capabilities and time reversibility. In this respect, Tsai et al. [TKL04] present a comprehensive survey on symplectic integrators used in molecular dynamics.
In a different class of their own, positionbased dynamics (see Bender et al. [BMOT13] for a survey) offer a workaround for any force or impulse based system, avoiding the intricacies of using integrators. This family of methods is, however, inaccurate for scenarios where velocities and forces need to be measured, hence we mention it as an alternative for special effects applications.
3 Constrained object animation
To better understand the elements involved in most physicsbased animation scenarios, we will briefly introduce a simple constraint enforcing system. Supporting this kind of simulation mechanism requires a discrete geometrical sampling of the initial shape of the object. For exemplification purposes, we consider a deformable linear object, embedded in a 3D space, whose crosssection is considerably smaller than its length. Its elastic properties that determine how its discrete structure changes are implemented by adding geometric constraints enforced by potentials (see the work of Teschner et al. [THMG04] for a more detailed and generalized application). These potentials provide a direct measure of how the structure of points differs from its initial configuration.
Generally, if are the vertices of a constrained group, then a constraint function is defined as follows:
(1) 
This function incorporates additional information relating the current vertex configuration of the group to the initial geometrical image through specific scalar measures of length, area, angle, volume, etc. For example, a length based constraint incorporates the initial or rest lengths of directly connected vertices as a reference for the measure of deviation. Generally, an energy function produces only positive amounts and can be written as:
(2) 
For a configuration, we now consider the restriction of the energy function at a node . This function is written as:
(3) 
Since the gradient vector coincides with the direction of maximum potential increase, it is natural to consider a penalty function that points in the opposite direction:
(4) 
This last equation can be written in an equivalent form:
(5) 
The simplest example of such a behaviour is the case of a linear spring connecting two points and . The natural constraint function is defined as:
(6) 
where denotes the spring’s rest length. From this constraint function, a corresponding elastic deformation potential energy can be expressed as:
(7) 
We define the following operator:
(8) 
also known in some works as the variational or functional derivative operator.
An elastic spring force at node can be derived using equation 7:
After applying several derivative computation rules (chain rule and the derivative of a product of two functions), we find the general elastic force expression:
(9) 
where is the linear spring’s stiffness coefficient.
A deformable object defined as a closed surface can be discretized by dividing its interior volume into tetrahedral cells. Apart from linear springs, it is desirable to enforce local constraints aimed at preserving the volumes under deformation. Such forces are easy to introduce by deriving them from a volume preserving constraint function involving tetrahedral cells:
(10) 
where is the volume of the tetrahedron in its rest configuration. Subsequently, the volume preserving potential is:
(11) 
Using the same method as for linear springs, we compute the force at the vertex by using the operator:
(12) 
After conveniently arranging the results from the derivation of equation (12), we can write down the expanded formula of this force:
(13) 
where is an addedin stiffness coefficient.
Examining figure 1, we can see how the volumetric force acts to prevent volume changes. For example, if point is shifted and the volume increases as a consequence, the direction of the volumetric force is given by the cross product vector . Inconsistent or ”flipped” cell configurations, as well as degenerate tetrahedra can and are likely to be encountered during a simulation. Compressible bodies prevent themselves from being completely flattened by acting as nonlinear spring elements. When their state is close to a collapse, the elastic forces should increase asymptotically towards infinity. Since the volumetric force component acts like a vertex to face spring in our system, we add a nonlinear spring component, as suggested in [CM97]. Hence, if is the current volume of a tetrahedron, then we can write the expression of the improved volumetric force as:
(14) 
where is the sign function defined on the set of real values.
Another scenario that can prevent the simulation from recovering its initial rest shape in the absence of perturbing forces is the inversion phenomenon. The authors of [ST08] present an improved mechanism, capable of handling the inversion of a tetrahedral structure for finite element simulations. This process relies on finding a rotation that best aligns the deformed cell with its undeformed counterpart and then deriving penalty forces. While we could have used a similar approach, the above modification works for cases where negative volumes are reported. The elastic penalty forces described by equation (14) are capable of acting against the inversion. Additionally, other constraintbased forces can be derived (e.g. area preserving forces for the triangular faces of a tetrahedron or angle preserving forces).
4 Deformable linear object model
Using constrained mass point configuration, we can now describe the steps required to build a tetrahedral cellbased deformable object along the geometric image of a support curve:

A curve discretization: where is a sample 3D point.

A set of frames: define the orientation vectors , and . At each point , a coordinate frame is attached.

Volumetric cells: for each pair of neighbouring vertices, , three tetrahedra are constructed: , , and (as depicted in figure 2 ).
The tetrahedral cells allow retrieving consistent information about local torsion and curvature changes from one link to another adjacent segment. A direct reference for how the object twists around the line is given by the relative orientation of to . Structural resistance is added by substituting the edges of the tetrahedra with linear springs. These constraints tend to act like curvature springs since they oppose the bending of the object around the line.
To account for plausible twisting behaviour, we introduce a quaternion based constraint system that acts like a torsional spring at each node of the object. Considering there are three connected nodes, and , in this order, the torsional spring tends to reposition the point such that the resulting configuration is closer to the initial, rest configuration. Such a process requires finding a suitable axis and computing the relative orientation of the vector with respect to the and vectors. We achieve this behaviour by computing a torsion compensating quaternion:
(15) 
where , . The angle represents the angle between the and vectors with respect to the axis , and the superscript designates values for the initial, undeformed state of the object. The torsion quaternion, can then be used to rotate the vector to minimize the torsion offset. Instead of directly rotating this vector, a force will be applied to the node, thus mimicking the effects of an angular spring. Figure 3 depicts how the torsion quaternions are derived with respect to local geometry.
5 Simulation update logic
In general, the steps required for a complete update of the object’s state can be assembled as follows:

Numerical integration: compute the current acceleration from the current state , and add any force contributions from the collision solver stage. Using an explicit integration method, the new state is computed.

Approximate forces: . These values are to be used in the collision response stage.

Correct positions: applying a position based dynamics constraint projection, the linear springs are additionally relaxed. Such a process is equivalent to an iterative GaussSeidel solver and it requires translating the object vertices by small displacements. This step helps satisfying stiff constraints and avoids outstretching artifacts.

Collision resolution: pairwise link collisions are detected and response forces and impulses are derived. The positions and velocities are corrected, storing response forces in accumulators, . These force residues are fed back to the numerical integration scheme and are included in the acceleration component the next iteration will use.
These steps must then be performed in this order, leading to a plausible update effect at interactive framerates.
6 Results
We have used our model to simulate a simple laparoscopic suturing task. Given a physical substrate, a thread was driven through a tissuelike structure (see figure 4). The wire model we used is ideal for testing the behaviour of an integrator where both stiff and soft constraints drive the simulation. A simple methodology was employed to compute a score sheet for each integration scheme. We tracked the speed and stability of several methods while varying the timestep. The results were quantized by considering the explicit Euler method as a reference. The order of accuracy is not particularly important as it is orthogonal to our stability goals. Due to their energyconserving properties, symplectic integrators are favoured when competing with other methods that achieve similar scores.
Criteria  

Method  Max  Time  Symplectic  Order  Score 
Explicit Euler  0.0098s  8ms  NO  1.225  
Symplectic Euler  0.0294s  8ms  YES  3.675  
Midpoint  0.0153s  8ms  NO  1.9125  
HalfStep  0.0168s  10ms  NO  1.68  
Verlet  0.023s  10ms  YES  2.3  
ForestRuth  0.021s  12ms  YES  1.75  
Symplectic Midpoint  0.028s  8ms  YES  3.5  
RungeKutta 4  0.027s  14ms  NO  1.928  
Moified HalfStep  0.0216s  10ms  NO  2.16 
Analyzing table 1, the symplectic Euler method (discussed by Cromer [Cro81]) is the most stable for the total iteration time it needs (the score is obtained as the ratio between the maximum time step and the time needed to execute one simulation iteration). We modified the Midpoint method to achieve symplecticity, obtaining the second highest score. Given the fact that the Midpoint family of methods is accurate to the second order for the position terms (while being a first order method for the velocity terms), we recommend it for applications where accuracy is of some importance. The Verlet method ( [Ver67] ), popular for molecular dynamics simulations, also benefits from its relatively high stability, symplecticity and second order accuracy. A third order method, the ForestRuth integration scheme [FR90]
, represents the best alternative for applications where accuracy is a key element. The RungeKutta fourth order method, although supporting relatively high time steps, is the most time consuming and probably not a good choice for real time applications.
Although symplectic integrators excel in scenarios where their energy preserving features are central (e.g. where only conservative forces are involved), we have found this family of integrators to outperform their explicit integrator counterparts. As a last remark, we have also modified the Half Step method ([HNW93]
) to support a semiimplicit update for the middle estimation
. This modification significantly increases the method’s stability, as seen in table 1.As an additional benchmark study, we have used the pendulum equation, . The explicit Euler is clearly the most unstable, introducing ghost energies as seen in the phase space diagram comparison with respect to the midpoint method in figure 4(a). On the other hand, in figure 4(b), the symplectic Euler and Midpoint methods have a much higher stability range, with the latter yielding slightly lower energy variations. As second order explicit integrators, the modified half step method is also more stable than the original version (as depicted in figure 5(a)). Nevertheless, both methods tend to add ghost energies, but are much more stable than the first order explicit Euler. As the order of accuracy increases, larger time steps can be used (as it is the case with RungeKutta methods), but the performance impact does not justify such a tradeoff. Finally, for models where no damping or nonconservative forces are involved, the ForestRuth and Verlet methods are the natural choices (comparative phase space plots in figure 5(b)). Nevertheless, in case accuracy can be sacrificed, we still recommend using the first order symplectic Euler method as it has the highest stability versus complexity score.
7 Conclusion
In this work we have discussed the implications of employing an explicit integration method for updating a soft or semirigid body simulation. For real time applications where accuracy is not a goal, we recommend using a symplectic integrator as it has the best performance and stability score. Even when nonconservative forces are involved, this class of integrators is able to cope with stiff constraints. However, for applications where accuracy cannot be sacrificed, both the RungeKutta or the ForestRuth integrators can be used. The latter choice is accurate up to the third order and conserves energy, while the RungeKutta progressively loses small amounts of energy, counting for a slight increase in stability.
On a final note, implicit integrators, mentioned in the related work section 2, are not the usual choice for real time applications, thus motivating our investigation towards an explicit integrator alternative.
Appendix A Explicit Integration Methods
In the following appendix, we present the minor modifications of the explicit integration methods that were tested in our simulation application. These numerical methods make use of an acceleration function, , a fixed timestep, , and compute new positions and velocities, given the previous state .
Symplectic Midpoint  Modified Half Step 

References
 [BMOT13] Jan Bender, Matthias Müller, Miguel A. Otaduy, and Matthias Teschner. Positionbased methods for the simulation of solid objects in computer graphics. In EUROGRAPHICS 2013 State of the Art Reports. Eurographics Association, 2013.
 [BW98] David Baraff and Andrew Witkin. Large steps in cloth simulation. In Proceedings of the 25th annual conference on Computer graphics and interactive techniques, SIGGRAPH ’98, pages 43–54, New York, NY, USA, 1998. ACM.
 [CM97] Lee Cooper Cooper and Steve Maddock. Preventing collapse within massspringdamper models of deformable objects. In In The Fifth Internation Conference in Central Europe on Computer Graphics and Visualization, pages 70–78, 1997.
 [Cro81] Alan Cromer. Stable solutions using the euler approximation. American Journal of Physics, 49(5):455–459, 1981.
 [DSB99] Mathieu Desbrun, Peter Schröder, and Alan Barr. Interactive animation of structured deformable objects. In Proceedings of the 1999 conference on Graphics interface ’99, pages 1–8, San Francisco, CA, USA, 1999. Morgan Kaufmann Publishers Inc.
 [FR90] E. Forest and R. D. Ruth. Fourthorder symplectic integration. Phys. D, 43(1):105–117, April 1990.
 [FSH10] Basil Fierz, Jonas Spillmann, and Matthias Harders. Stable explicit integration of deformable objects by filtering high modal frequencies. Journal of WSCG, 18(13):81–88, February 2010.
 [HE01] Michael Hauth and Olaf Etzmuß. A high performance solver for the animation of deformable objects using advanced numerical methods. Comput. Graph. Forum, 20(3):319–328, 2001.

[HNW93]
E. Hairer, S. P. Nørsett, and G. Wanner.
Solving ordinary differential equations I (2nd revised. ed.): nonstiff problems
. SpringerVerlag New York, Inc., New York, NY, USA, 1993.  [Jim12] P. Jiménez. Survey on modelbased manipulation planning of deformable objects. Robot. Comput.Integr. Manuf., 28(2):154–163, April 2012.
 [MGC07] Jaruwan Mesit, Ratan K. Guha, and Shafaq Chaudhry. 3d soft body simulation using massspring system with internal pressure force and simplified implicit integration. JCP, 2(8):34–43, 2007.
 [NMK06] Andrew Nealen, Matthias Müller, Richard Keiser, Eddy Boxerman, and Mark Carlson. Physically based deformable models in computer graphics. Computer Graphics Forum, 25(4):809–836, 2006.
 [SKZF11] Craig A. Schroeder, Nipun Kwatra, Wen Zheng, and Ronald Fedkiw. Asynchronous evolution for fullyimplicit and semiimplicit time integration. Comput. Graph. Forum, 30(7):1983–1992, 2011.
 [SLM06] Martin Servin, C. Lacoursiére, and N. Melin. Interactive simulation of elastic deformable materials. In SIGRAD 2006 Conference Proceedings, pages 22–32, 2006.
 [ST08] Ruediger Schmedding and Matthias Teschner. Inversion handling for stable deformable modeling. Vis. Comput., 24(7):625–633, July 2008.
 [THMG04] Matthias Teschner, Bruno Heidelberger, Matthias Muller, and Markus Gross. A versatile and robust model for geometrically complex deformable solids. In Proceedings of the Computer Graphics International, CGI ’04, pages 312–319, Washington, DC, USA, 2004. IEEE Computer Society.
 [TKL04] ShanHo Tsai, M. Krech, and D.P. Landau. Symplectic integration methods in molecular and spin dynamics simulations. Brazilian Journal of Physics, 34:384 – 391, 06 2004.
 [Ver67] Loup Verlet. Computer ”experiments” on classical fluids. i. thermodynamical properties of lennardjones molecules. Phys. Rev., 159:98–103, Jul 1967.
 [VMT05] P. Volino and N. MagnenatThalmann. Accurate garment prototyping and simulation. ComputerAided Design & Applications, 2(5):645–654, 2005.
Comments
There are no comments yet.