Next: 6.4 Simplifications Up: 6 Numerical Implementation Previous: 6.2 Calculation of Element

## 6.3 Discrete Dynamics Equations

In order to integrate the D-NURBS ordinary differential equations of motion (17) in an interactive modeling environment, it is important to provide the modeler or designer with visual feedback about the evolving state of the dynamic model. Rather than using costly time integration methods that take the largest possible time steps, it is more crucial to provide a smooth animation by maintaining the continuity of the dynamics from one step to the next. Hence, less costly yet stable time integration methods that take modest time steps are desirable.

The matrices , , and (and , , and ) are symmetric, sparse, and banded. Several algorithms are available for the numerical integration of the D-NURBS ordinary differential equations of motion. The suitability of implicit or explicit integration algorithms is dependent on the bandwidth of the matrices, as determined by the dimensionality of the parametric space and the order of the NURBS basis functions. The matrices for a D-NURBS curve have a single band which has a half-bandwidth of 4k, where k is the order of the NURBS basis. For D-NURBS surfaces, the matrices become block banded, with each block containing n bands similar to those of dynamic curves, where n depends on the order of the NURBS basis in the opposite parametric direction.

We integrate the differential equations (17) through time by discretizing the derivative of over time-steps . The state of the D-NURBS at time is integrated using prior states at time t and . Depending on the choice of physical parameters, (17) may be a stiff system. We use an implicit time integration method in order to maintain the stability of the integration scheme. The implicit method employs discrete derivatives of using backward differences

Making use of the fact that , we obtain the time integration formula

where the superscripts denote evaluation of the quantities at the indicated times, and where the remaining quantities are evaluated at time . For example, we can extrapolate the mass matrix using the formula

and likewise for the other matrices and vectors in (29). The simpler, constant extrapolations , etc., ([15] Section 8.6) also work satisfactorily.

In the interest of efficiency, we do not factorize the matrix expression on the left hand side of (29) in order to solve for . Instead, we employ the conjugate gradient method to obtain an iterative solution [27, 32]. To achieve interactive simulation rates, we limit the number of conjugate gradient iterations per time step to 10. We have observed that 2 iterations typically suffice to converge to a residual of less than . More than 2 iterations tend to be necessary when the physical parameters (mass, damping, tension, stiffness, applied forces) are changed dramatically during interactive sculpting.

Note that when physical parameter values are chosen such that the equations (17) are not stiff, it is much cheaper to employ an explicit time integration method using forward differences. Appendix D discusses the forward difference approach. Note that the explicit method requires values for the matrices only at time t, hence (30) is not needed.

For the D-NURBS curve, we simply replace with in (29) and everything proceeds as in the case of surfaces.

In the case of D-NURBS with linear constraints, we discretize the derivatives of (rather than ). Analogous to (29), the discrete version of (22) is

Since there are fewer degrees of freedom in than in , faster numerical implementation of constrained D-NURBS is possible, provided the constraint matrix is sparse. Note that since the conjugate gradient algorithm requires only gradient vectors, we need not compute , and explicitly. The only extra cost is the computation of and the multiplication of with several vectors in (31).

For nonlinear constraints, at each time step we can apply the conjugate gradient algorithm to solve (27) for the Lagrange multipliers and the constrained generalized accelerations (given known and ). We then integrate and from t to to obtain the constrained generalized velocities and coordinates (e.g., using the simple Euler method ; ).

Next: 6.4 Simplifications Up: 6 Numerical Implementation Previous: 6.2 Calculation of Element

Demetri Terzopoulos | Source Reference