# Position Based Dynamics

## Abstract

The most popular approaches for the simulation of dynamic systems in computer graphics are force based. Internal and external forces are accumulated from which accelerations are computed based on Newton’s second law of motion. A time integration method is then used to update the velocities and finally the positions of the object. A few simulation methods (most rigid body simulators) use impulse based dynamics and directly manipulate velocities. In this paper we present an approach which omit the velocity layer as well and immediately works on the positions. The main advantage of a position based approach is its controllability. Overshooting problems of explicit integration schemes in force based systems can be avoided. In addition, collision constraints can be handled easily and penetrations can be resolved completely by projecting points to valid locations. We have used the approach to build a real time cloth simulator which is part of a physics software library for games. This application demonstrates the strengths and benefits of the method.

## 1. Introduction

Research in the field of physically based animation in computer graphics is concerned with finding new methods for the simulation of physical phenomena such as the dynamics of rigid bodies, deformable objects or fluid flow. In contrast to computational sciences where the main focus is on accuracy, the main issues here are stability, robustness and speed while the results should remain visually plausible. Therefore, existing methods from computational sciences can not be adopted one to one. In fact, the main justification for doing research on physically based simulation in computer graphics is to come up with specialized methods, tailored to the particular needs in the field. The method we present falls into this category.

The traditional approach to simulating dynamic objects has been to work with forces. At the beginning of each time step, internal and external forces are accumulated. Examples of internal forces are elastic forces in deformable objects or viscosity and pressure forces in fluids. Gravity and collision forces are examples of external forces. Newton’s second law of motion relates forces to accelerations via the mass. So using the density or lumped masses of vertices, the forces are transformed into accelerations. Any time integration scheme can then be used to first compute the velocities from the accelerations and then the positions from the velocities. Some approaches use impulses instead of forces to control the animation. Because impulses directly change velocities, one level of integration can be skipped.

In computer graphics and especially in computer games it is often desirable to have direct control over positions of objects or vertices of a mesh. The user might want to attach a vertex to a kinematic object or make sure the vertex always stays outside a colliding object. The method we propose here works directly on positions which makes such manipulations easy. In addition, with the position based approach it is possible to control the integration directly thereby avoiding overshooting and energy gain problems in connection with explicit integration. So the main features and advantages of position based dynamics are

• Position based simulation gives control over explicit integration and removes the typical instability problems.
• Positions of vertices and parts of objects can directly be manipulated during the simulation.
• The formulation we propose allows the handling of general constraints in the position based setting.
• The explicit position based solver is easy to understand and implement.

• 基于位置的仿真给出了显式积分的控制，消除了典型的不稳定性问题。
• 在模拟过程中，顶点和对象部分的位置可以直接被操纵。
• 我们建议的公式允许在基于位置的设置中处理一般的约束。
• 基于明确位置的solver很容易理解和实现。

## 3. Position Based Simulation

In this section we will formulate the general position based approach. With cloth simulation, we will give a particular application of the method in the subsequent and in the results section. We consider a three dimensional world. However, the approach works equally well in two dimensions.

### 3.1. Algorithm Overview

We represent a dynamic object by a set of $$N$$ vertices and $$M$$ constraints. A vertex $$i\in [1,\cdots ,N]$$ has a mass $$m_i$$, a position $$\mathbf x_i$$ and a velocity $$\mathbf v_i$$.A constraint $$j\in [1,\ldots ,M]$$ consists of

• a cardinality $$n_j$$,
基数(和该约束相关的顶点的个数)$$n_j$$
• a function $$C_j:\mathbb R^{3n_j} \to \mathbb R$$,
约束函数$$C_j:\mathbb R^{3n_j} \to \mathbb R$$
• a set of indices $$\{i_1,\ldots i_{n_j}\}, i_k \in [1,\ldots N]$$,
相关顶点索引集合$$\{i_1,\ldots i_{n_j}\}, i_k \in [1,\ldots N]$$
• a stiffness parameter $$k_j \in [0 \ldots 1]$$ and
刚度系数$$k_j \in [0 \ldots 1]$$
• a type of either equality or inequality.
以及约束类型，equality 或者 inequality

Constraint $$j$$ with type equality is satisfied if $$C_j(\mathbf x_{i_1},\ldots ,x_{i_{n_j}}) = 0$$. If its type is inequality then it is satisfied if $$C_j(\mathbf x_{i_1},\ldots ,x_{i_{n_j}}) \ge 0$$. The stiffness parameter $$k_j$$ defines the strength of the constraint in a range from zero to one.

Based on this data and a time step $$\Delta t$$, the dynamic object is simulated as follows:

(01) forall vertices $$i$$

(02) 　initialize $$\mathbf x_i = \mathbf x^0_i,\mathbf v_i = \mathbf v^0_i,w_i = 1/m_i$$

(03) endfor

(04) loop

(05) 　forall vertices $$i$$ do $$\mathbf v_i \gets \mathbf v_i +\Delta t w_i \mathbf f_{ext}(\mathbf x_i)$$

(06) 　dampVelocities($$\mathbf v_1,\ldots ,\mathbf v_N$$)

(07) 　forall vertices $$i$$ do $$\mathbf p_i \gets \mathbf x_i + \Delta t \mathbf v_i$$

(08) 　forall vertices $$i$$ do generateCollisionConstraints($$\mathbf x_i \to \mathbf p_i$$)

(09) 　loop solverIterations times

(10) 　　projectConstraints($$C_1,\ldots ,C_{M+M_{coll}}, \mathbf p_1,\ldots ,\mathbf p_N$$)

(11) 　endloop

(12) 　forall vertices $$i$$

(13) 　　$$\mathbf v_i \gets (\mathbf p_i - \mathbf x_i)/ \Delta t$$

(14) 　　$$\mathbf x_i \gets \mathbf p_i$$

(15) 　endfor

(16) 　velocityUpdate($$\mathbf v_1,\ldots,\mathbf v_N$$)

(17) endloop

Lines (1)-(3) just initialize the state variables. The core idea of position based dynamics is shown in lines (7), (9)-(11) and (13)-(14). In line (7), estimates $$\mathbf p_i$$ for new locations of the vertices are computed using an explicit Euler integration step. The iterative solver (9)-(11) manipulates these position estimates such that they satisfy the constraints. It does this by repeatedly project each constraint in a GaussSeidel type fashion (see Section 3.2). In steps (13) and (14), the positions of the vertices are moved to the optimized estimates and the velocities are updated accordingly. This is in exact correspondence with a Verlet integration step and a modification of the current position [Jak01], because the Verlet method stores the velocity implicitly as the difference between the current and the last position. However, working with velocities allows for a more intuitive way of manipulating them.

The velocities are manipulated in line (5), (6) and (16). Line (5) allows to hook up external forces to the system if some of the forces cannot be converted to positional constraints. We only use it to add gravity to the system in which case the line becomes $$\mathbf v_i \gets \mathbf v_i + \Delta t\mathbf g$$, where $$\mathbf g$$ is the gravitational acceleration. In line (6), the velocities can be damped if this is necessary. In Section 3.5 we show how to add global damping without influencing the rigid body modes of the object. Finally, in line (16), the velocities of colliding vertices are modified according to friction and restitution coefficients.

The given constraints $$C_1,\ldots ,C_M$$ are fixed throughout the simulation. In addition to these constraints, line (8) generates the $$M_{coll}$$ collision constraints which change from time step to time step. The projection step in line (10) considers both, the fixed and the collision constraints.

The scheme is unconditionally stable. This is because the integration steps (13) and (14) do not extrapolate blindly into the future as traditional explicit schemes do but move the vertices to a physically valid configuration $$\mathbf p_i$$ computed by the constraint solver. The only possible source for instabilities is the solver itself which uses the Newton-Raphson method to solve for valid positions (see Section 3.3). However, its stability does not depend on the time step size but on the shape of the constraint functions.

The integration does not fall clearly into the category of implicit or explicit schemes. If only one solver iteration is performed per time step, it looks more like an explicit scheme. By increasing the number of iterations, however, a constrained system can be made arbitrarily stiff and the algorithm behaves more like an implicit scheme. Increasing the number of iterations shifts the bottleneck from collision detection to the solver.

### 3.2. The Solver

The input to the solver are the $$M+M_{coll}$$ constraints and the estimates $$\mathbf p_1,\ldots,\mathbf p_N$$ for the new locations of the points. The solver tries to modify the estimates such that they satisfy all the constraints. The resulting system of equations is non-linear. Even a simple distance constraint $$C(\mathbf p_1,\mathbf p_2) = |\mathbf p_1 - \mathbf p_2| - d$$ yields a non-linear equation. In addition, the constraints of type inequality yield inequalities. To solve such a general set of equations and inequalities, we use a Gauss-Seidel-type iteration. The original Gauss-Seidel algorithm (GS) can only handle linear system. The part we borrow from GS is the idea of solving each constraint independently one after the other. However, in contrast to GS,solving a constraint is a non linear operation. We repeatedly iterate through all the constraints and project the particles to valid locations with respect to the given constraint alone. In contrast to a Jacobi-type iteration, modifications to point locations immediately get visible to the process. This speeds up convergence significantly because pressure waves can propagate through the material in a single solver step, an effect which is dependent on the order in which constraints are solved. In over-constrained situations, the process can lead to oscillations if the order is not kept constant.
solver的输入是$$M+M_{coll}$$个约束，和作为顶点新位移的一组预估位置$$\mathbf p_1,\ldots,\mathbf p_N$$。solver会尝试修改这些预估位置，使其满足所有约束条件。由此产生的方程组是非线性的。即使是一个简单的距离约束$$C(\mathbf p_1,\mathbf p_2) = |\mathbf p_1 - \mathbf p_2| - d$$也会生成一个非线性方程。此外，类型为inequality的约束也会产生不等式。为了解决这样一组方程和不等式，我们使用了Gauss-Seidel类型的迭代法。原始的Gauss-Seidel算法(GS)只能解决线性系统问题。我们从GS中借用的思想，是相继地独立解决每个约束(GS算法是将前面求解的结果应用到未求解的等式中，应该是这个意思吧)。然而，与GS相比，求解约束是一个非线性的操作。我们反复遍历所有的约束，并将粒子投影到有效的位置，仅考虑给定的约束条件。与雅可比迭代法相比，对点位置的修改可以立即对过程可见。这大大加快了收敛速度，因为压力波可以在一个单独的求解过程中通过材料传播，这种效应依赖于约束求解的顺序。在过度约束的情况下，如果顺序不保持不变，则过程可能导致振荡。

### 3.3. Constraint Projection

Projecting a set of points according to a constraint means moving the points such that they satisfy the constraint. The most important issue in connection with moving points directly inside a simulation loop is the conservation of linear and angular momentum. Let $$\Delta \mathbf p_i$$ be the displacement of vertex $$i$$ by the projection. Linear momentum is conserved if

$\sum_i m_i \Delta \mathbf p_i = \mathbf 0 \tag{1}$

which amounts to conserving the center of mass. Angular momentum is conserved if

$\sum_i \mathbf r_i \times m_i \Delta \mathbf p_i = \mathbf 0 \tag{2}$

where the $$\mathbf r_i$$ are the distances of the $$\mathbf p_i$$ to an arbitrary common rotation center. If a projection violates one of these constraints, it introduces so called ghost forces which act like external forces dragging and rotation the object. However, only internal constraints need to conserve the momenta. Collision or attachment constraints are allowed to have global effects on the object.

The method we propose for constraint projection conserves both momenta for internal constraints. Again, the point based approach is more direct in that we can directly use the constraint function while force based methods derive forces via an energy term (see [BW98, THMG04]). Let us look at a constraint with cardinality $$n$$ on the points $$\mathbf p_1, \ldots ,\mathbf p_n$$ with constraint function $$C$$ and stiffness $$k$$. We let $$\mathbf p$$ be the concatenation $$[\mathbf p^T_1, \ldots , \mathbf p^T_n]^T$$. For internal constraints, $$C$$ is independent of rigid body modes, i.e. translation and rotation. This means that rotating or translating the points does not change the value of the constraint function. Therefore, the gradient $$\nabla_{\mathbf p}C$$ is perpendicular to rigid body modes because it is the direction of maximal change. If the correction $$\Delta \mathbf p$$ is chosen to be along $$\nabla_{\mathbf p}C$$ both momenta are automatically conserved if all masses are equal (we handle different masses later). Given $$\mathbf p$$ we want to find a correction $${\Delta \mathbf p}$$ such that $$C(\mathbf p + \Delta \mathbf p) = 0$$. This equation can be approximated by

$C(\mathbf p + \Delta \mathbf p) \approx C(\mathbf p) + \nabla_{\mathbf p}C(\mathbf p) \cdot \Delta \mathbf p = 0 \tag{3}$

Restricting $$\Delta \mathbf p$$ to be in the direction of $$\nabla_{\mathbf p}C$$ means choosing a scalar $$\lambda$$ such that

$\Delta \mathbf p = \lambda \nabla_{\mathbf p}C(\mathbf p) \tag{4}$

Substituting Eq.(4) into Eq.(3), solving for $$\lambda$$ and substituting it back into Eq.(4) yields the final formula for $$\Delta \mathbf p$$

$\Delta \mathbf p = - \frac{C(\mathbf p)}{|\nabla_{\mathbf p}C(\mathbf p)|^2} \nabla_{\mathbf p}C(\mathbf p) \tag{5}$

which is a regular Newton-Raphson step for the iterative solution of the non-linear equation given by a single constraint. For the correction of an individual point $$\mathbf p_i$$ we have

$\Delta \mathbf p_i = -s\nabla_{\mathbf p_i}C(\mathbf p_1,\ldots,\mathbf p_n), \tag{6}$

where the scaling factor

$s = \frac{C(\mathbf p_1,\ldots,\mathbf p_n)}{\sum_j| \nabla_{\mathbf p_j}C(\mathbf p_1,\ldots,\mathbf p_n)|^2} \tag{7}$

is the same for all points. If the points have individual masses, we weight the corrections $$\Delta \mathbf p_i$$ by the inverse masses $$w_i=1/m_i$$. In this case a point with infinite mass, i.e. wi = 0, does not move for example as expected. Now Eq. (4) is replaced by

$$\Delta \mathbf p_i = \lambda w_i \nabla_{\mathbf p_i} C(\mathbf p)$$ yielding

$s = \frac{C(\mathbf p_1,\ldots,\mathbf p_n)}{\sum_j w_j| \nabla_{\mathbf p_j}C(\mathbf p_1,\ldots,\mathbf p_n)|^2} \tag{8}$

for the scaling factor and for the final correction

$\Delta \mathbf p_i = -sw_i\nabla_{\mathbf p_i}C(\mathbf p_1,\ldots,\mathbf p_n). \tag{9}$

To give an example, let us consider the distance constraint function $$C(\mathbf p_1, \mathbf p_2)=|\mathbf p_1 - \mathbf p_2|-d$$. The derivative with respect to the points are $$\nabla_{\mathbf p_1}C(\mathbf p_1, \mathbf p_2)=\mathbf n$$ and $$\nabla_{\mathbf p_2}C(\mathbf p_1, \mathbf p_2)=\mathbf -\mathbf n$$ with $$\mathbf n=\frac{\mathbf p_1-\mathbf p_2}{|\mathbf p_1-\mathbf p_2|}$$. The scaling factor $$s$$ is, thus, $$s = \frac{|\mathbf p_1 - \mathbf p_2|-d}{w_1+w_2}$$ and the final corrections

$\Delta \mathbf p_1 = -\frac{w_1}{w_1+w_2}(|\mathbf p_1 - \mathbf p_2|-d)\frac{\mathbf p_1-\mathbf p_2}{|\mathbf p_1-\mathbf p_2|} \tag{10}$

$\Delta \mathbf p_2 = +\frac{w_1}{w_1+w_2}(|\mathbf p_1 - \mathbf p_2|-d)\frac{\mathbf p_1-\mathbf p_2}{|\mathbf p_1-\mathbf p_2|} \tag{11}$

which are the formulas proposed in [Jak01] for the projection of distance constraints (see Figure 2). They pop up as a
special case of the general constraint projection method.

We have not considered the type and the stiffness $$k$$ of the constraint so far. Type handling is straight forward. If the type is equality we always perform a projection. If the type is inequality, the projection is only performed if $$C(\mathbf p_1,\ldots,\mathbf p_n)<0$$. There are several ways of incorporating the stiffness parameter. The simplest variant is to multiply the corrections $$\Delta \mathbf p$$ by $$k\in[0 \ldots 1]$$. However, for multiple iteration loops of the solver, the effect of $$k$$ is non-linear. The remaining error for a single distance constraint after $$n_s$$ solver iterations is $$\Delta \mathbf p(1-k)^{n_s}$$. To get a linear relationship we multiply the corrections not by $$k$$ directly but by $$k'=1-(1-k)^{1/n_s}$$. With this transformation the error becomes $$\Delta \mathbf p(1-k')^{n_s}=\Delta \mathbf p(1-k)$$ and, thus, becomes linearly dependent on $$k$$ and independent of $$n_s$$ as desired. However, the resulting material stiffness is still dependent on the time step of the simulation. Real time environments typically use fixed time steps in which case this dependency is not problematic.

### 3.4. Collision Detection and Response

One advantage of the position based approach is how simply collision response can be realized. In line (8) of the simulation algorithm the $$M_{coll}$$ collision constraints are generated. While the first $$M$$ constraints given by the object representation are fixed throughout the simulation, the additional $$M_{coll}$$ constraints are generated from scratch at each time step. The number of collision constraints $$M_{coll}$$ varies and depends on the number of colliding vertices. Both, continuous and static collisions can be handled. For continuous collision handling, we test for each vertex $$i$$ the ray $$\mathbf x_i \to \mathbf p_i$$. If this ray enters an object, we compute the entry point $$\mathbf q_c$$ and the surface normal $$\mathbf n_c$$ at this position. An inequality constraint with constraint function $$C(\mathbf p) = (\mathbf p - \mathbf q_c) \cdot \mathbf n_c$$ and stiffness $$k = 1$$ is added to the list of constraints. If the ray $$\mathbf x_i \to \mathbf p_i$$ lies completely inside an object, continuous collision detection has failed at some point. In this case we fall back to static collision handling. We compute the surface point $$\mathbf q_s$$ which is closest to $$\mathbf p_i$$ and the surface normal $$\mathbf n_s$$ at this position. An inequality constraint with constraint function $$C(\mathbf p) = (\mathbf p - \mathbf q_s) \cdot \mathbf n_s$$ and stiffness $$k=1$$ is added to the list of constraints. Collision constraint generation is done outside of the solver loop. This makes the simulation much faster. There are certain scenarios, however, where collisions can be missed if the solver works with a fixed collision constraint set. Fortunately, according to our experience, the artifacts are negligible.

Friction and restitution can be handled by manipulating the velocities of colliding vertices in step (16) of the algorithm. The velocity of each vertex for which a collision constraint has been generated is dampened perpendicular to the collision normal and reflected in the direction of the collision normal.

The collision handling discussed above is only correct for collisions with static objects because no impulse is transferred to the collision partners. Correct response for two dynamic colliding objects can be achieved by simulating both objects with our simulator, i.e. the $$N$$ vertices and $$M$$ constraints which are the input to our algorithm simply represent two or more independent objects. Then, if a point $$\mathbf q$$ of one objects moves through a triangle $$\mathbf p_1, \mathbf p_2, \mathbf p_3$$ of another object, we insert an inequality constraint with constraint function $$C(\mathbf q, \mathbf p_1, \mathbf p_2, \mathbf p_3) = \pm (\mathbf q - \mathbf p_1) \cdot [(\mathbf p_2 - \mathbf p_1) \times (\mathbf p_3 - \mathbf p_1)]$$ which keeps the point $$\mathbf q$$ on the correct side of the triangle. Since this constraint function is independent of rigid body modes, it will correctly conserve linear and angular momentum. Collision detection gets slightly more involved because the four vertices are represented by rays $$\mathbf x_i \to \mathbf p_i$$. Therefore the collision of a moving point against a moving triangle needs to be detected (see section about cloth self collision).

### 3.5. Damping

In line (6) of the simulation algorithm the velocities are dampened before they are used for the prediction of the new positions. Any form of damping can be used and many methods for damping have been proposed in the literature (see [NMK05]). Here we propose a new method with some interesting properties:

(1) $$\mathbf x_{cm} = (\sum_i \mathbf x_i m_i)/(\sum_i m_i)$$

(2) $$\mathbf v_{cm} = (\sum_i \mathbf v_i m_i)/(\sum_i m_i)$$

(3) $$\mathbf L = \sum_i \mathbf r_i \times (m_i \mathbf v_i)$$

(4) $$\mathbf I = \sum_i \mathbf{\tilde r}_i \mathbf{\tilde r}^T_i m_i$$

(5) $$\omega = \mathbf I^{-1}\mathbf L$$

(6) forall vertices $$i$$

(7) 　$$\Delta \mathbf v_i = \mathbf v_{cm} + \omega \times \mathbf r_i - \mathbf v_i$$

(8) 　$$\mathbf v_i \gets \mathbf v_i + k_{\text{damping}}\Delta \mathbf v_i$$

(9) endfor

Here $$\mathbf r_i = \mathbf x_i - \mathbf x_{cm}$$, $$\mathbf {\tilde r}_i$$ is the 3 by 3 matrix with the property $$\mathbf {\tilde r}_i \mathbf v = \mathbf r_i \times \mathbf v$$, and $$k_{\text{damping}} \in [0 \ldots 1]$$ is the damping coefficient. In lines (1)-(5) we compute the global linear velocity $$\mathbf x_{cm}$$ and angular velocity $$\omega$$ of the system. Lines (6)-(9) then only damp the individual deviations $$\Delta v_i$$ of the velocities $$\mathbf v_i$$ from the global motion $$\mathbf v_{cm} + \omega \times \mathbf r_i$$. Thus, in the extreme case $$k_{\text{damping}} = 1$$ , only the global motion survives and the set of vertices behaves like a rigid body. For arbitrary values of $$k_{\text{damping}}$$, the velocities are globally dampened but without influencing the global motion of the vertices.

### 3.6. Attachments

With the position based approach, attaching vertices to static or kinematic objects is quite simple. The position of the vertex is simply set to the static target position or updated at every time step to coincide with the position of the kinematic object. To make sure other constraints containing this vertex do not move it, its inverse mass $$w_i$$ is set to zero.

## 4. Cloth Simulation

We have used the point based dynamics framework to implement a real time cloth simulator for games. In this section we will discuss cloth specific issues thereby giving concrete examples of the general concepts introduced in the previous section.

### 4.1. Representation of Cloth

Our cloth simulator accepts as input arbitrary triangle meshes. The only restriction we impose on the input mesh is that it represents a manifold, i.e. each edge is shared by at most two triangles. Each node of the mesh becomes a simulated vertex. The user provides a density $$\rho$$ given in mass per area [$$kg=m^2$$]. The mass of a vertex is set to the sum of one third of the mass of each adjacent triangle. For each edge, we generate a stretching constraint with constraint function

$C_{stretch}(\mathbf p_1, \mathbf p_2)=|\mathbf p_1 - \mathbf p_2|-l_0,$

stiffness $$k_{stretch}$$ and type equality. The scalar $$l_0$$ is the initial length of the edge and $$k_{stretch}$$ is a global parameter provided by the user. It defines the stretching stiffness of the cloth. For each pair of adjacent triangles $$(\mathbf p_1, \mathbf p_3, \mathbf p_2)$$ and $$(\mathbf p_1, \mathbf p_2, \mathbf p_4)$$ we generate a bending constraint with constraint function

$C_{bend}(\mathbf p_1, \mathbf p_2, \mathbf p_3, \mathbf p_4) = acos\left(\frac{(\mathbf p_2 - \mathbf p_1) \times (\mathbf p_3 - \mathbf p_1)}{|(\mathbf p_2 - \mathbf p_1) \times (\mathbf p_3 - \mathbf p_1)|} \cdot \frac{(\mathbf p_2 - \mathbf p_1) \times (\mathbf p_4 - \mathbf p_1)}{|(\mathbf p_2 - \mathbf p_1) \times (\mathbf p_4 - \mathbf p_1)|}\right)-\varphi_0,$

stiffness $$k_{bend}$$ and type equality. The scalar $$\varphi_0$$ is the initial dihedral angle between the two triangles and $$k_{bend}$$ is a global user parameter defining the bending stiffness of the cloth (see Figure 4). The advantage of this bending term over adding a distance constraint between points $$\mathbf p_3$$ and $$\mathbf p_4$$ or over the bending term proposed by [GHDS03] is that it is independent of stretching. This is because the term is independent of edge lengths. This way, the user can specify cloth with low stretching stiffness but high bending resistance for instance (see Figure 3).

Eqns. (10) and (11) define the projection for the stretching constraints. In the appendix A we derive the formulas to project the bending constraints.

### 4.2. Collision with Rigid Bodies

For collision handling with rigid bodies we proceed as described in Section 3.4. To get two-way interactions, we apply an impulse $$m_i \Delta \mathbf{p}_i / \Delta t$$ to the rigid body at the contact point, each time vertex $$i$$ is projected due to collision with that body. Testing only cloth vertices for collisions is not enough because small rigid bodies can fall through large cloth triangles. Therefore, collisions of the convex corners of rigid bodies against the cloth triangles are also tested.

### 4.3. Self Collision

Assuming that the triangles all have about the same size, we use spatial hashing to find vertex triangle collisions [THM03]. If a vertex $$\mathbf q$$ moves through a triangle $$\mathbf p_1$$, $$\mathbf p_2$$, $$\mathbf p_3$$, we use the constraint function

$C(\mathbf q,\mathbf p_1, \mathbf p_2, \mathbf p_3) = (\mathbf q-\mathbf p_1) \cdot \frac{(\mathbf p_2 - \mathbf p_1)\times (\mathbf p_3 - \mathbf p_1)}{|(\mathbf p_2 - \mathbf p_1)\times (\mathbf p_3 - \mathbf p_1)|} -h \tag{12}$

where $$h$$ is the cloth thickness (see Figure 5). If the vertex enters from below with respect to the triangle normal, the constraint function has to be
$$h$$是布的厚度(见图5)。如果顶点从三角形背面沿着法线方向穿过，那么约束函数应该是

$C(\mathbf q,\mathbf p_1, \mathbf p_2, \mathbf p_3) = (\mathbf q-\mathbf p_1) \cdot \frac{(\mathbf p_3 - \mathbf p_1)\times (\mathbf p_2 - \mathbf p_1)}{|(\mathbf p_3 - \mathbf p_1)\times (\mathbf p_2 - \mathbf p_1)|} -h \tag{13}$

to keep the vertex on the original side. Projecting these constraints conserves linear and angular momentum which is essential for cloth self collision since it is an internal process. Figure 6 shows a rest state of a piece of cloth with self collisions. Testing continuous collisions is insufficient if cloth gets into a tangled state, so methods like the ones proposed by [BWK03] have to be applied.

### 4.4. Cloth Balloons

For closed triangle meshes, overpressure inside the mesh can easily be modeled (see Figure 7). We add an equality constraint concerning all $$N$$ vertices of the mesh with constraint function

$C(\mathbf p_1, \ldots , \mathbf p_N) = \left(\sum^{n_{\text{triangles}}}_{i=1}(\mathbf p_{t^i_1} \times \mathbf p_{t^i_2}) \cdot \mathbf p_{t^i_3}\right) - k_{\text{pressure}}V_0 \tag{14}$

and stiffness $$k=1$$ to the set of constraints. Here $$t^i_1$$,$$t^i_2$$ and $$t^i_3$$ are the three indices of the vertices belonging to triangle $$i$$. The sum computes the actual volume of the closed mesh. It is compared against the original volume $$V_0$$ times the overpressure factor $$k_{\text{pressure}}$$. This constraint function yields the gradients

$\nabla_{\mathbf p_i}C = \sum_{j:t^i_1=i}(\mathbf p_{t^j_2} \times \mathbf p_{t^j_3}) + \sum_{j:t^i_2=i}(\mathbf p_{t^j_3} \times \mathbf p_{t^j_1}) + \sum_{j:t^i_3=i}(\mathbf p_{t^j_1} \times \mathbf p_{t^j_2}) \tag{15}$

These gradients have to be scaled by the scaling factor given in Eq. (7) and weighted by the masses according to Eq. (9) to get the final projection offsets $$\Delta \mathbf p_i$$.

## References

[THM03] TESCHNER M., HEIDELBERGER B., MüLLER M., POMERANERTS D., GROSS M.: Optimized spatial hashing for collision detection of deformable objects. Proc. Vision, Modeling, Visualization VMV 2003 (2003), 47–54

[BWK03] BARAFF D., WITKIN A., KASS M.: Untangling cloth. In Proceedings of the ACM SIGGRAPH (2003), pp. 862–870.

## Appendix A:

### Bending Constraint Projection

Position Based Dynamics【译】

(1)
(2)

0条