Our main goal for today is to extend our discussion about particles from last class and model systems in which the particles may be connected. This will allow us to animate cloth (see the demo at the end of the notes), which is what you'll implement in this week's lab.

Spring-mass systems

For some systems, such as an elastic string or cloth, we can model particles as being connected to each other through a mini spring. For example, we might model an elastic string as a bunch of springs attached to each other (see the figure above). The external force acting on each particle (a vertex in the mesh) would then depend on the position of the neighboring springs and is described using Hooke's law. Hooke's law states that the force $\vec f$ needed to compress a spring by a displacement of $\Delta \vec x$ is proportional to $\Delta \vec x$. The proportionality coefficient $k$ is a constant that is a property of the spring. $\vec f = k \Delta \vec x$. The direction of this external force is in the direction of the displacement, hence the spring exerts a force opposite the displacement. Applying Newton's second law to the $i$-th particle in a general spring-mass system (similar to the one in the image above) gives:

$$ m \vec a_i = m\vec g + k \left( \ell_{i,i+1} - \lVert\vec p_{i+1} - \vec p_{i}\rVert \right)\cdot\frac{(\vec p_{i+1} - \vec p_{i})}{\lVert\vec p_{i+1} - \vec p_{i}\rVert} + k \left( \ell_{i,i-1} - \lVert\vec p_{i-1} - \vec p_{i}\rVert\right)\cdot\frac{(\vec p_{i-1} - \vec p_{i})}{\lVert\vec p_{i-1} - \vec p_{i}\rVert} $$

where $k$ is the stiffness of the spring (here, assumed constant for each spring) and $\ell_{ij}$ is the length of the spring at rest. Note that the direction of the spring force is directly opposite to the direction of the spring displacement with respect to the original spring length. Careful attention is needed at the boundaries of the domain. Depending on the problem setup, one of these will likely be a wall, so $\vec p_{\mathrm{wall}} = 0$. The other end might also be fixed (think of holding a piece of string at both ends) or it might be free, in which case we could approximate this boundary condition as $\vec p_{\mathrm{end}} = \vec p_{\mathrm{end}-1}$.

In the case of a cloth simulation, we have a grid of particles connected by springs, i.e. a 2d mesh instead of a 1d mesh we had above. The only extra complication is that instead of only having structural springs between neighboring vertices, we also have shear and bending springs:

• structural spring: creates a force between a vertex at $(i,j)$ and any neighboring vertices $(i \pm 1,j)$, $(i,j \pm 1)$.
• shear spring: creates force on a vertex at $(i,j)$ caused by a neighboring vertex along the diagonals surrounding vertex $(i,j)$.
• bending spring: creates a force on a vertex at $(i,j)$ due to bending using vertices at $(i \pm 2, j)$ and $(i,j\pm 2)$.

We also need to be careful in the way we model the forces arising from springs near the cloth boundaries. Of course, you can add external forces such as gravity, wind, etc. to your equations of motion.

(from Interactive Computer Graphics, Angel)

There are other ways to do cloth simulations with triangle meshes (instead of a grid). Many of these techniques also make use of adaptive methods (see this article), whereby anisotropic triangles are placed in regions in which the cloth currently exhibits a lot of curvature. Adaptation has the benefit of reducing computation time since you can potentially use fewer triangles to resolve the cloth motion with the same accuracy. The difficulty arises in how to automatically determine where these triangles should be placed (and how to orient them) - this is related to my research so please feel free to ask me about it!

A simpler way to model cloth with Verlet integration

It can be tricky to model cloth with the methods above because there's a lot of bookkeeping needed to keep track of the forces between neighboring particles, and we need to be careful with how the equations are applied at boundaries. Also, separately updating position and velocity like we did last class can be error prone. Remember, we're not trying to model the physics perfectly, we just want something to look okay (the previous methods aren't exactly physically correct anyway). Let's continue with the spring-mass model, but this time, we'll model the springs as a bunch of constraints and also use Verlet integration to update the particle positions. Instead of updating the position and velocity by computing the right-hand-side of our system of ODEs, Verlet integration updates the positions directly using:

$$ \vec p^{k+1} = 2\vec p^k - \vec p^{k-1} + \frac{\sum\vec{f}_{\mathrm{ext}}\Delta t^2}{m}. $$

where $\vec{x}^{k+1}$ is the position update for the next time step $(k+1)$, $\vec p^k$ is the current position and $\vec p^{k-1}$ is the previous position of the particle (mesh vertex). You could also add more damping by using 1.99 (instead of 2) in the factor on the current position, and subtracting 0.99 times the previous position (instead of 1). We will update the position of each vertex by first moving it according to the external forces (such as gravity or wind), and then trying to satisfy the constraints. We will try to satisfy these constraints by computing how much each spring is stretched, and then moving both end point vertices half of this stretched distance to compute their new location. The way we compute this displacement is important, and we'll use a method that is not only fast, but also adds enough damping to create a realistic-looking cloth. For a spring (constraint) with endpoint vertices $\vec p$ and $\vec q$, we can compute the normalized deflection $\delta$ as

$$ \delta = \frac{\ell - \ell_0}{\ell}, \quad \ell = \lVert\vec q - \vec p\rVert, $$

where $\ell_0$ is the original length of the spring at rest. We can then update each endpoint vertex of the spring according to

$$ \vec p^{\prime} = \vec p + \alpha_p ( \vec q - \vec p), \quad \vec q^{\prime} = \vec q - \alpha_q(\vec q - \vec p). $$

where $\alpha_p$ and $\alpha_q$ represent the displacement applied to each endpoint particle. A simple choice would be to use $\alpha_p = \alpha_q = \frac{1}{2}\delta$, however, heavier particles should be moved less than lighter particles. Instead, we can calculate the fraction of the displacement using the mass of each particle:

$$ \alpha_p = \frac{m_q}{(m_p + m_q)}\delta, \quad\quad \alpha_q = \frac{m_p}{(m_p + m_q)}\delta. $$

We would then immediately set $\vec p \gets \vec p^{\prime}$ and $\vec q \gets \vec q^{\prime}$. Although this sounds inconsistent because we are setting new positions of spring vertices before computing the update for other vertices, this actually works really well and is known as relaxation. By doing only 1 or 2 iterations of the (1) update based on external forces and then (2) update to satisfy constraints, we can actually get a pretty realistic looking cloth! The overall algorithm then looks like:

// cloth update algorithm for one time step (frame)
for numIter iterations: // numIter is usually about 1 or 2
  for each particle:
    particle.forceUpdate(); // compute new position using external forces and Verlet integration
  for each constraint:
    p = constraint.p; // first endpoint vertex (vec3)
    q = constraint.q; // second endpoint vertex (vec3)
    delta = (||q - p|| - L0)/||q - p||; // scalar
    dx = delta * (q - p); // vec3
    mt = mp * mq / (mp + mq); // mp and mq are the masses of particles p and q
    constraint.p += dx * mt / mp;
    constraint.q -= dx * mt / mq;

With this scheme, you can also impose certain vertices to be fixed with this weighting: just set $1/m$ to zero (i.e. $m \to \infty$, a very heavy particle).

Here is an example of a cloth animation using the scheme described above, which you will implement in this week's lab. Try clicking one of the particles (dots) and pull it around. The animation will then re-run and it should look like a swaying cloth!

Other systems to model.

We can also model hair using the techniques we saw today in which each strand of hair is a collection of springs. Based on the type of hair, certain types of springs (structural, shear or bending) may have more of an influence on the motion of the strand. For a recent research paper about modeling tightly-coiled hair, please see this project.

© Philip Claude Caplan, 2023 (Last updated: 2023-12-07)