FIELD OF THE INVENTION

This disclosure is directed to simulating the motion of thin and long objects (â€œstrandsâ€ such as hair) in the field of computer graphics and more specifically for animation.

BACKGROUND OF THE INVENTION

The technical problem of simulating by computer graphics the dynamics of strands in terms of the forward dynamics problem for rigid multi-body serial chains is well known; see S. Hadap, Hair Simulation, (Ph.D. Thesis, Dept. of Informatics, Univ. of Geneva, Switzerland, 2003). One method uses recursive algorithms for solving the forward dynamics problem of rigid multi-body serial chains. Such algorithms are well known; see R. Featherstone, Robot Dynamics Algorithms, (Kluwer Academic Publishers, 1987). Recursive and residual algorithms for the simulation of multi-body systems are also known; see J. I. Rodriguez, J. M. Jimenez, F. J. Funes, J. G. de Jalon, â€œRecursive and Residual Algorithms for the Efficient Numerical Integration of Multi-Body Systems,â€ Multibody System Dynamics, 11(4), (Springer Netherlands, May 2004): 295-320.

SUMMARY

Examples of objects whose motion can be animated well using the present method and apparatus include long hair, wavy or curly hair, braids, animal ears, tails, foliage, and jewelry. The method is implemented in computer software, of the type that animation artists use in production of computer animated motion pictures or videos.

The present method advantageously simulates accurately the physical effects of inertia, centrifugal and Coriolis forces, gravity, wind or other force fields on strands. The method handles correctly non-zero bend and twist in the rest shape of the strands, as well as changes in bend and twist during the simulation. In many systems for simulating strands the joint constraints between the strand elements (â€œlinksâ€ ) are satisfied only approximately, a procedure which is prone to numerical errors. The present method enforces these constraints exactly. The method uses a formulation based on the minimal set of coordinates needed to specify the strands' configuration. This simplifies the task of specifying the rest shape and ensures optimal use of computing resources. The method employs algorithms for kinematics and dynamics computations that are optimal; the number of operations grows linearly in the number of links. The method correctly computes collisions between strands and scene geometry. The method simulates correctly the effects of constraints and collisions, including static and dynamic friction associated with collisions. The method uses analytical methods for simulating the effects of constraints. Using the present method, it is possible to art-direct intuitively and effectively the physically based simulation of the animated objects.

In the present method, the equations of motion for strands are formulated as a system of differential algebraic equations (DAE). Analytical constraints are computed as external forces, and thus the user is able to impose analytical constraints in the framework of DAE. An approximate representation of static and dynamic friction as unilateral constraints is also provided. For solving the system of DAEs corresponding to strand dynamics, the DAE solver software DASPK well known in the computer field is used in one embodiment.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a flowchart of the present method.

FIG. 2 shows a strand as a multi-body system.

FIG. 3 shows constraints acting on the FIG. 2 system.

FIG. 4 shows collisions with an obstacle.

FIG. 5 shows an apparatus for performing the FIG. 1 method.

DETAILED DESCRIPTION

The present method uses a formulation of strand dynamics in terms of Differential Algebraic Equations (â€œDAEâ€ ). DAE are well known; see K. E. Brenan, S. L. Campbell, L. R. Petzold, Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations, (SIAM, 1996). The motion equations corresponding to strands are highly non-linear and have large numerical stiffness. Non-linear behavior is due to Coriolis forces and the non-linear elastic model used to limit high deformations. Numerical stiffness in the simulation of strand dynamics is caused by small rotational inertia due to thin geometry, large bending and torsional stiffness-to-mass ratio, and non-straight rest shape.

The DAE approach is used here to improve the numerical stability of the strand system in the presence of non-linearity and numerical stiffness. Improved numerical stability allows larger time steps in the simulation, with a corresponding improvement of the system's interactive response time. FIG. 1 thus shows the present method 12 in flowchart form, as would be carried out by a computer program (software) in the context of animation of strands. In the first block 16, a conventional data structure is provided representing the strands. For instance, the strands are represented as instances of a class (â€˜classâ€™ is used here in the sense of an object-oriented programming language construct, see for example B. Stroustroup â€œThe C++ Programing Languageâ€ , Addison-Wesley 2001). The data that is being grouped by the class that represents a strand includes:

- position, velocity, acceleration, orientation, angular velocity and angular acceleration of the strand as a whole.
- position, velocity, acceleration, orientation, angular velocity and angular acceleration of an animated strand that is used as reference for art-direction.
- parameters used to implement the elasticity model of the strand: torsional stiffness and anisotropic bend stiffness coefficients.
- a vector of instances of a class representing the segments that make up the strand.
- a vector of data structures representing joints between segments.
- a vector of data structures representing joint velocities.
  Next at 20, motion paths in space are defined for the base link (see FIG. 2) for each strand of interest. Next at 22, possible dynamic forces are computed (see FIG. 3). Then at 26, collisions with any obstacles are detected (see FIG. 4) Then at 30, the response forces from the detected collisions are computed. Then at 34, the DAE solver is invoked to solve the resulting strand simulation for the current animation frame. At 38, it is determined if any additional frames are to be simulated (animated). Detail of each block of FIG. 1 is provided hereinafter.

DAEs are a well known general form of differential equations, where the equations are expressed in implicit rather than explicit form:

F(y,{dot over (y)},t)=0 â€ƒâ€ƒ(1)

Compare Eq. (1) with a differential equation given in the more familiar explicit form:

{umlaut over (q)}=M^âˆ’1G(q, {dot over (q)}, t) â€ƒâ€ƒ(2)

where the highest order derivatives are separated on the left hand side of the equations. Because of this freedom, DAEs form a much larger class of differential equations than the class of explicit differential equations.

The dynamic state of a strand is completely described for purposes of animation by the state vector:

$\begin{matrix} y = [\begin{matrix} q \\ \dot{q} \end{matrix}] & (3) \end{matrix}$

where q are is a vector of joint variables and {dot over (q)} is a vector of the joint velocities. Both vectors have size n, where n is the total number of degrees of freedom allowed by the joints which couple adjacent strand links. In general, the total number of degrees of freedom n is different from the number of links m, because joints may have more than one degree of freedom.

The DAE solver advances the solution from state [y, {dot over (y)}], to [y, {dot over (y)}]_t+1. The well known DASPK software is used as the DAE solver in one embodiment. The DAE solver requires a definition of the residue function:

rs(t)â‰¡F(y(t), {dot over (y)}(t), t) â€ƒâ€ƒ(4)

There are 2n equations in the DAE system describing the dynamics of a strand: the first n are the force balance equations for the n links, see Equations (8) and (10) below. The last n are consistency conditions between

$y = [\begin{matrix} q \\ \dot{q} \end{matrix}] î¢ž î¢ž and î¢ž î¢ž \dot{y} = [\begin{matrix} \dot{q} \\ \overset{Â¨}{q} \end{matrix}],$

due to {dot over (q)} appearing in both.

A recursive formulation of the kinematic and dynamic equations is used here. This is well known; see R. Cottle, J. S. Pang, R. Stone, The Linear Complementarity Problem, (Academic Press, 1992). The recursive formulation has the advantages that it achieves linear time in n, while using reduced coordinates and thus satisfying exactly all the constraints between the bodies making up the strands.

In the recursive formulation, the equations of motion of the system are expressed only implicitly as recurrence relations between physical quantities. The computation of the residual function rs(t) has two recursive steps: forward kinematics: given the velocity and acceleration of the root link v₀, a₀and the state vector [q, {dot over (q)}] determine the velocities and accelerations of all links; and forward dynamics: given all the above plus the inertia tensors of the links and external forces determine the residual (unbalanced) forces for all links.

To describe the dynamics of the multi-body system, the present method uses spatial algebra and the associated well known notation, see R. Featherstone, Robot Dynamics Algorithms, (Kluwer Academic Publishers, 1987). Consider a serial branched multi-body system (a strand) having m links connected by m joints as shown in FIG. 2. There are no loops in this multi-body system. The base link is denoted by link₀and is typically constrained to a motion path (see FIG. 1). This means that the velocity and acceleration v₀, a₀of the base link are specified by the animation system instead of being a result of the dynamics simulation. The other links are numbered in a first manner, i.e. ensuring that the parent of a link has always a smaller index than the link itself.

The kinematics equations in spatial notation between link_jand its parent link_iare:

{circumflex over (v)}_j={circumflex over (X)}_j{circumflex over (v)}_i+Å _j{dot over (q)}_jâ€ƒâ€ƒ(5)

Ã¢_j={circumflex over (X)}_jÃ¢_i+Å _j{umlaut over (q)}_j+{circumflex over (v)}_j{circumflex over (Ã—)}Å _j{dot over (q)}_jâ€ƒâ€ƒ(6)

The spatial velocity {circumflex over (v)}_j, spatial acceleration Ã¢_jand spatial inertia matrix ÃŽ_jare all defined in the local frame {circumflex over (F)}_jthat is rigidly attached to the principal inertia axis of link_j. The joint variables q_jdefine the spatial transformation {circumflex over (X)}_jthat may be used to transform spatial quantities defined in the parent's frame {circumflex over (F)}_ito the frame {circumflex over (F)}_j. The spatial matrix of the joint axis connecting link_jto its parent link is Å _j. The computation of the forward kinematics of the system amounts to applying Equations (5-6) for j=1 . . . m while starting from the given velocity {circumflex over (v)}₀and acceleration Ã¢₀of the root link.

In the DAE framework, computing forward dynamics means computing a scalar residue value for each degree of freedom of a joint. These residue values are associated with the links' accelerations. First, the force balance condition for each link are used to obtain residual forces. This computation is again recursive and, due to data dependencies, starts from the outermost links and proceeds towards the base link.

Consider one of the outermost links and denote it by link_m. The force balance equation for link_mis:

{circumflex over (r)}Å _m=ÃŽ_mÃ¢_m+{circumflex over (v)}_m{circumflex over (Ã—)}ÃŽ_m{circumflex over (v)}_mâˆ’{circumflex over (f)}_mâ€ƒâ€ƒ(7)

where {circumflex over (f)}_mis the external spatial force acting on link_mand {circumflex over (r)}Å _mis the spatial force residue.

The residue for link_mis obtained by projecting the force residue {circumflex over (r)}Å _monto the associated joint's motion subspace defined by the spatial joint axis Å _m. The joint actuation force associated with the strand deformation model Q_mis accounted for, which also acts in the joint's motion subspace (the deformation model is explained below):

rs_m=Å _m^S{circumflex over (r)}Å _mâˆ’Q_mâ€ƒâ€ƒ(8)

Here rs_mhas a number of components equal to the dimension of the joint's motion subspace, i.e. to the number of degrees of freedom of the joint connecting link_mto the rest of the multi-body system.

The force balance equation for link_i, which is the parent of link_j, is:

$\begin{matrix} {\hat{rs}}_{i} = {\hat{I}}_{i} î¢ž {\hat{a}}_{i} + {\hat{v}}_{i} î¢ž \hat{Ã—} î¢ž {\hat{I}}_{i} î¢ž {\hat{v}}_{i} - {\hat{f}}_{i} + \underset{j}{âˆ‘} î¢ž {\hat{X}}_{j}^{- 1} î¢ž {\hat{rs}}_{j} & (9) \end{matrix}$

where {circumflex over (F)}_iis the external spatial force acting on link_i. Note that the balance equation accounts for the effect of the spatial force residual of link_j. To bring this force in link_i's frame the inverse of the transform {circumflex over (X)}_jis used. The sum is over the children of link_i. As above, the DAE residual is the projection of the spatial force residual in the joint's motion space:

rs_i=Å _i^S{circumflex over (r)}Å _iâˆ’Q_iâ€ƒâ€ƒ(10)

This process is repeated for each link, progressing towards the base link of the multi-body system.

The residue associated with joint velocities is simply given by consistency conditions in the representation of the dynamic state of the strand. There is {dot over (q)} appearing in both y and {dot over (y)}, which leads to:

rs_i+n={dot over (q)}_i*âˆ’{dot over (q)}_iâ€ƒâ€ƒ(11)

where {dot over (q)}_iand {dot over (q)}_iare the joint velocities in y and {dot over (y)}, respectively.

The base joint in the multi-body system in FIG. 2 is typically constrained to some prescribed motion path. In practice, it is often necessary to impose some additional constraints on the system, e.g. constraining a point on some other link to a motion path, or constraining a joint to an animated value in time, see FIG. 3. These constraints (â€œconstraint pathâ€ in FIG. 3) are transient in nature and often introduce cyclic dependency in the system. The effects of such additional constraints are computed by augmenting the DAE based multi-body formulation. A method for computing the effects of such constraints is known, see E. Kokkevis, â€œPractical Physics for Articulated Characters,â€ Proceedings of Game Developers Conference 2004 (2004) and belongs to a class of techniques called analytical methods. The main problem with this known methodology is that it is not compatible with a DAE solver. A DAE solver integrates the system without explicitly calculating the state derivatives, whereas the known method directly alters the explicitly calculated derivatives to satisfy analytical constraints.

The present approach computes analytical constraints as external forces. This allows imposing analytical constraints in the framework of DAE. Consider a point constraint p_jas depicted in FIG. 3. The trajectory (path) of p_j, starting with the initial conditions, can be uniquely defined by the time varying acceleration a_j^dn_j. As discussed above, here one does not directly evaluate the state derivative vectors {dot over (y)}_t+iin order to integrate the system. Therefore, one cannot simply enforce the acceleration constraint, by directly altering the state derivatives {dot over (y)}_t+ias known. The present method enforces the constraint by applying an external force instead.

Before every DAE integration step, the precise nature of the force f_jn_j, is determined analytically using the above known method. The unit constraint direction n_jis treated as constant and is updated for every integration step. There is a linear relationship between the magnitude of the applied force f_jand the resulting desired acceleration a_j^d:

$\begin{matrix} a_{j}^{d} = \frac{âˆ‚ a_{j}}{âˆ‚ f_{j}} î¢ž f_{j} + a_{j}^{0} & (12) \end{matrix}$

where, a_j⁰is the acceleration in the direction n_jbefore the force is applied.

If there is another constraint point p_iwith force having magnitude f_iin the direction n_i, the resulting accelerations a_i^dand a_j^dwill be given by the following linear system:

$\begin{matrix} [\begin{matrix} a_{i}^{d} \\ a_{j}^{d} \end{matrix}] = [\begin{matrix} \frac{âˆ‚ a_{i}}{âˆ‚ f_{i}} & \frac{âˆ‚ a_{i}}{âˆ‚ f_{j}} \\ \frac{âˆ‚ a_{j}}{âˆ‚ f_{i}} & \frac{âˆ‚ a_{j}}{âˆ‚ f_{j}} \end{matrix}] + [\begin{matrix} a_{i}^{0} \\ a_{j}^{0} \end{matrix}] & (13) \end{matrix}$

Generalizing, for m such constraints one needs to determine the vector of f Îµ R^mof unknown force magnitudes by solving the following linear system:

Kf+a⁰=0. â€ƒâ€ƒ(14)

The Jacobian matrix K Îµ R^mÃ—m(well known in the field) can be evaluated by applying unit test forces at each constraint and evaluating the changes in accelerations. An efficient procedure to evaluate the Jacobian using the framework of Featherstone's Articulated Body Method is in Kokkevis, referenced above. The constraint forces thus determined are applied to the multibody system over the next integration step via the external force {circumflex over (f)}, as disclosed above.

One replaces the constraint direction n by a spatial vector {circumflex over (n)} to generalize the type of the constraints that can be represented. For example, {circumflex over (n)} can represent also joint acceleration constraints. One further extends the constraint system to include unilateral constraints such as collisions, friction and joint limits by posing it as a Linear Complementarity Problem (LCP) which is well known in the field. An LCP solver takes a linear system written in the form:

Mz+q=w â€ƒâ€ƒ(15)

and solves for z such that for every i:

z_iâ‰§0, w_iâ‰§0 and z_iw_i=0. â€ƒâ€ƒ(16)

The following abbreviates these conditions as:

z_iâ‰§0 complementary to w_iâ‰§0. â€ƒâ€ƒ(17)

The LCP of the present constraint system can then be written as:

Kf+a⁰âˆ’a^dâ‰§0 complementary to fâ‰§0 â€ƒâ€ƒ(18)

The LCP states that the forces f, applied only in positive constraint direction, would strive to make the resulting constraint accelerations a equal to desired acceleration a^d. However, force f_iwill be zero if the resulting constraint acceleration a_iis greater than desired acceleration a_i^d. Iterative methods are used to solve the LCP.

Using forces for enforcing the constraints has an advantage here. If the forces are indeterminate, they are projected into the multi-body's motion null-space, thus always giving valid configuration, without any jumps or jittering in the simulation. Further, as the forces are determined analytically, as compared to, e.g., penalty-based formulation, they are small for most types of the constraints. Thus they are well within the stability zone of the integrator when taking 4 to 8 time-steps per animation frame. The only exception to this is the velocity impulse constraint, discussed below. As the constraint may not be satisfied accurately, one augments the constraint acceleration by a proportional-derivative form. To exemplify, for a positional constraint, the desired acceleration and the constraint direction are:

$\begin{matrix} a_{i}^{d} = {\overset{Â¨}{p}}_{i}^{d} + K_{p} î¢ (p_{i}^{d} - p_{i}) + K_{d} î¢ ({\dot{p}}_{i}^{d} - v_{i}) & (19) \\ a_{i}^{d} = ï˜\dots a_{i}^{d} ï˜† & (20) \\ n_{i}^{d} = \frac{a_{i}^{d}}{a_{i}^{d}} & (21) \end{matrix}$

where {umlaut over (p)}_i^d, {dot over (p)}_i^d, p_i^dare respectively the acceleration, velocity and position of the constraint path and p_i, v_iare the current position and velocity of the constraint.

The unilateral position constraints disclosed above are used to develop a collision response model for the multi-body system. Collision detection is well known in computer graphics. For brevity, the following describes only the requirements from the collision detection system as needed for the present method. Between the current configuration given by the state vector y_tand the extrapolated configuration using the derivative vector {dot over (y)}_tand integration time step h, one finds all the points on the multi-body system that would collide with the obstacles. FIG. 4 shows two such collision positionsâ€”point p_ihas already penetrated the obstacle and point p_jis about to collide. There may be more than one colliding point for a link. Let n_ibe the collision normal, with direction pointing away from the obstacle, and a_iand v_ibe the collision accelerations and velocities respectively, relative to the obstacle.

One applies collision response in two stepsâ€”contacts and impacts. First, the unilateral constraints are computed that would prevent collision points from accelerating towards the obstacle. This is followed by computation of velocity impulses that would prevent collision points from moving towards the obstacle.

One can decompose the collision acceleration and the collision velocity into the normal components a_ni, v_niand tangential components a_ti, v_ti. To prevent any acceleration towards the obstacle, a unilateral constraint is inserted along the collision normal direction n_i. The unilateral constraint corresponding to the friction acts in the tangent plane defined by the collision normal. Instead of sampling the tangent plane into discrete set of tangents to formulate a complex and numerically expensive friction model based of the discrete frictional cone, the present method uses a technique as follows. The unilateral constraint direction corresponding to friction is set as:

$\begin{matrix} t_{i} = unit î¢ (a_{ti} + \frac{v_{ti}}{h}) & (22) \end{matrix}$

If both a_tiand v_tiare zero, the previously determined tangent vector is used for the contact.

Finally, the LCP formulation corresponding to the two unilateral constraints in the direction n_iand t_iat collision position p_iis:

a_niâˆ’a_ni^dâ‰§0 complementary to f_niâ‰§0 â€ƒâ€ƒ(23)

Î¼_if_niâˆ’f_tiâ‰§0 complementary to Î»_iâ‰§0 â€ƒâ€ƒ(24)

(a_tiâˆ’a_ti^d)+Î»_iâ‰§0 complementary to f_tiâ‰§0 â€ƒâ€ƒ(25)

The desired normal acceleration a_ni^dis set to be proportional to the penetration depth d_iif the point is penetrating, see Equations 19-21, or zero if the collision point is outside the obstacle. The desired tangential acceleration a_ti^dis set to (âˆ’v_ti/h). The LCP formulation will compute the required amount of normal force f_nito remove the normal acceleration a_ni. The tangential force f_tiwill be at most Î¼_if_ni, and it will try to remove any tangential non-zero velocity component (the dynamic friction case) or if the tangential velocity is zero it will try to remove any tangential acceleration (the static friction case).

The present method uses impulses to arrest the collision normal velocity v_ni. Only those contacts are considered that have the normal velocity component v_ni<0. For the impulse computations, one can use the same acceleration constraints discussed above by setting:

a_ni^d=âˆ’(1+v)v_ni, â€ƒâ€ƒ(26)

where v is the coefficient of restitution associated with the collision. Instead of applying potentially large forces, one alters the joint velocities {dot over (q)}_t. Although theoretically this would invalidate the consistent solution set y_t, {dot over (y)}_t, in practice the solver is tolerant to the error.

The dynamic strand is composed from a hierarchy of input curves defined in a local frame, that defines the initial rest shape of the corresponding multi-body system. One embodiment of the present method provides the user (animator) with high-level control over the direct numerical simulation by means of relevant physical parameters of the dynamic strand, such as mass per unit length, strand radius, bend stiffness and damping, torsional stiffness and damping, gravity, and air drag.

The user can animate all the parameters and specify their length-wise anisotropic variation. The collision parameters coefficient of restitution, coefficient of static friction, and coefficient of dynamic friction are defined per obstacle surface. The strand may have additional anisotropic weights over the collision parameters, along with their length-wise variation.

Three-degrees-of-freedom spherical joints are used in one implementation of the present strand system. The joint variable of the i-th joint is expressed as a quaternion q_iand the joint velocity as a vector w_i. One decomposes the quaternion defining the relative transformation between two links into a twist component Î¸_taround the local y-axis and a pure bend component Î¸_baround a bend axis b. A nonlinear stiffness model is defined as follows:

$\begin{matrix} Q_{b} = K_{b} î¢ž b î¢ž î¢ž \tan î¢ (\frac{({Î¸}_{b} - {Î¸}_{b}^{0})}{2}) & (27) \\ Q_{t} = K_{t} î¢ž y î¢ž î¢ž \tan î¢ (\frac{({Î¸}_{t} - {Î¸}_{t}^{0})}{2}) & (28) \end{matrix}$

where Î¸_b⁰and Î¸_t⁰describe the rest configuration, K_tis the torsional stiffness coefficient and K_bis the anisotropic bend stiffness coefficient. The response model is almost linear for small deformations. However, the non-linear response prevents excessive deformations.

General external force fields are supported using the plug-in aspect of the well known Maya animation software and that of other animation software. The user can attach any complex combination of time-varying fields such as wake, turbulence, fluid simulations and general event driven scripted force fields. The user can further specify length-wise anisotropic weights for the external force fields. The user can optionally include these fields in computing the Jacobians numerically, as discussed above.

Often the dynamical simulation is driven by very extreme and brisk animated motions. Although the present method is a robust solution able to cope with such constraints, animators sometimes want the resulting motion to be selectively less violent (abrupt). A time scale Î² is provided here to control the amount of energy pumped into the strand system. Î² is a factor by which the velocity and acceleration of the base joint can be scaled and it is typically between zero and one; however the user can set it to values larger than one to accentuate the motion.

The local dynamics Î³ is another similar parameter which blends out velocity and acceleration of some local dynamics reference frame.

Ã¢₀=Î²(Ã¢₀âˆ’Î³Ã¢_ref) â€ƒâ€ƒ(29)

{circumflex over (v)}₀=Î²({circumflex over (v)}₀âˆ’Î³{circumflex over (v)}_ref) â€ƒâ€ƒ(30)

One scenario that is frequent in animation is a braid or long hairs of a character that fly away when the character starts running or rides a horse. In this case the local dynamics reference frame can be simply set to the character's hip joint, and with an appropriate value of the local dynamics parameter one can control the amount of flyaway in the result of the simulation.

Often the animator wants to reach a certain pose for an animated character at a certain time to make the character expressive. In order to support this posable dynamics requirement the present method inserts â€œspringsâ€ between a dynamic strand and the desired animated pose at the tip of each segment. This results in introducing penalty based constraints. The user will typically animate the stiffness and damping of the springs, to change gradually from a purely dynamics simulation to key-frame animation and back.

The present method is typically implemented in the form of computer software (code) executed on a computer. Coding the software, in light of this disclosure, would be routine; a suitable computer language is C+. FIG. 5 shows the structure of such a program 68 in terms of a group of modules each corresponding to activity in one of the steps of the FIG. 1 method. Thus in FIG. 5, memory 70 stores the strand data structures. A motion path component 74 (module) defines the motion path for the base link of each strand. A dynamics force module 78 evaluates the dynamic forces for each strand. A collision detector module 82 performs the collision detection for the strands, while collision response module 86 evaluates the collision response forces DAE solver 88 solves the differential algebraic equations for all strands in the current frame. The software itself would typically be stored on a computer readable medium (memory), such as a tape, disk, hard disk, flash memory, etc., for execution by the computer processor.

This disclosure is illustrative and not limiting; further modification will be apparent to those skilled in the art in light of this disclosure and are intended to fall within the scope of the appended claims.