For my CPSC 533 course project I implemented a simple spring and damper dynamics simulation in 3D. This system has several interesting features.
A major problem that can occur when working with systems of springs, dampers and point masses is a stiff system of equations. This typically happens when some of the springs are significantly stronger than the others. A system of equations is considered stiff when the maximimum step size that you can take is constrained by stability conditions and not by the error term.
There are essentially two ways of dealing with stiff systems of equations. The first way is to decrease the step size to meet the stability constraints. This solution is undesirable. Decreasing the step size means that more calculations have to be done to generate the same number of frames of animation. The second solution is to use a method with a larger region of absolute stability. The second solution is to use implicit integration methods. A-stable methods are numerical methods with regions of absolute stability that contain the entire left half-plane of z = h(lambda). Implicit integration methods are all A-stable. This makes them ideal for dealing with this sort of problem.
I implemented a backward Euler solver for my system based on the methods described by Baraff and Witkin[1]. This approach uses a Taylor series expansion of the force vector to achieve a first order approximation. My implementation works directly with equation 6 from [1] and does not use the condition modeldescribed later in that paper.
Because my backward Euler method dealt with equation 6 I was unable to use the conjugate gradient method for solving the resulting system of linear equations. The conjugate gradient method used in [1] requires a symmetric linear system. To solve the system of equations I used the LU decomposition algorithm found in Numerical Recipes[2].
The system allows some user interaction through the use of simple controllers. A number of binary sensors may be added to the system. The binary sensors return a one when the spring they are monitoring is within specified length bounds and returns zero otherwise. A single spring can have multiple sensors monitoring it at the same time.
The control manipulates the system by changing the characteristics of the springs. Every spring has an alternate rest lest, spring constant, and damper constant which it uses when the control takes effect. This method of control can be justified by the fact that there is some spring that could be added to the first spring that would give a combined result equivalent to the new values.
Every spring has a weight associated with each binary sensor. Weights are not constrained to any particular range. If the combined total of the weights exceeds one then the control for that spring is activated. The control must remain activated for a user specified number of frames before the spring will react. Each spring can have a different delay associated with it. Once a spring has been activated sufficiently long for the control to take effect it remains activated until after the weighted sum of the sensor inputs remains below one for the same delay before it can return to its normal state.
A simple collision model was implemented for this project. When a mass crosses through the ground it is forced to lie on the surface of the ground. Any objects at or below the surface of the ground have friction applied to them. Every vertex in the ground has a roughness parameter. Every mass in the system has two roughness parameters, one for moving forward and one for moving backward. When a vertex is not moving straight forwards or backwards the roughness is interpolated. The coefficient of friction is taken as the roughness of the ground multiplied by the roughness of the mass (wrt the masses direction of travel).
Every mass stores pointers to up to two other masses. One of these masses is considered forward of the current mass and the other mass is behind. When friction needs to be calculated the vector between the forward mass position the current mass position is added to the vector defined by the current mass position minus the backward mass position. The vector that results from this is normalised and used to specify the “front” of the mass. Taking the cross product of this vector with the normalised force vector gives cos(theta) where theta is the angle between the front of the current mass and the direction it is accelerating. To perform the interpolation of the mass’ roughness add one to cos(theta) and divide the result by two. This gives a function with the following properties:
These are all desirable properties for an interpolation function. The interpolated roughness is given by
forward_roughness * [cos(theta) + 1] / 2 + backward_roughness * {1 - [cos(theta) + 1] / 2}
Once the coefficient of friction has been determined a Coulomb friction model is used.
Other implicit methods such as the implicit midpoint or the trapezoidal method could be implemented. Because both of these methods are 2nd order, a closer Taylor series approximation to the force vector would be required to see a significant improvement. Unlike backward Euler, neither of these methods has stiff decay. This only becomes important if the equations are very stiff so it may not be a problem.
More sensor types would be a useful addition. Other binary sensors could detect ground contact or follow a target.
Linking the size of the time steps to the render time. Currently the simulation looks a little jumpy due to the differences between real time and simulated time.
Finally, a better user interface would be very useful. The current interface is pretty kludgy.
[1] Baraff, David and Andrew Witkin. “Large Steps in Cloth Simulation,” Proceedings of SIGGRAPH 1998, ACM SIGGRAPH, pp. 43-54.
[2] Press W., Teukolsky S., Vetterling W., Flannery B. (1992), “Numerical Recipes in C”, 2nd ed. Cambridge University Press, Cambridge.