Projective Dynamics

Mass-spring System And Implicit Method #

For node $i$ with mass $m_i$ experiencing force $\boldsymbol{f_i}$

which goes to

Gathering all nodes we have

where

and

Considering $\boldsymbol{x^{t+1}}$ as unkown variable to be found, we have

Construct a function

where $\partial E/\partial \boldsymbol{x}=-\boldsymbol{f}$, we will see that

so we are finding $\boldsymbol{x}$ such that $\frac{\partial F(\boldsymbol{x})}{\partial \boldsymbol{x}}=0$, which is also a necessary condition for local minimum of $F$. That is, if we find a $\boldsymbol{x}$ at which $F$ is locally minimal, we find a $\boldsymbol{x}$ such that $\frac{\partial F(\boldsymbol{x})}{\partial \boldsymbol{x}}=0$. Finding local minimum can be done by Newton-Raphson’s method.

We will need to calculate the hessian of $F$.

Position Based Dynamics #

High stiffness can cause problem when updating position using newton’s law (time and force). If we have another rule to directly manuplate position while conserving neccesary physical quantity (e.g. momentum), the problem can be solved. That is PBD. PBD is non-physical.

We will define the constraint as

where $d$ is the rest length between node $i$ and $j$. In general, we want all nodes satisfying constraints

But sadly in practice, we usually cannot satisfy all constraints at once, and must iterate over all constraints. Take one constraint $C_i$ for example. If the system has experienced an external force and moved to $\boldsymbol{x^{new}}$, the constraint is still expected to be satisfied

Here PBD makes an assumption: the direction of $\Delta\boldsymbol{x}$ is same with gradient of $\boldsymbol{C}_i$

Substitute it into previous equation and we will get

and

For the correction of an individual point $\boldsymbol{x_{k}}$, we have

where the scaling factor

If the points have individual masses, we weight the corrections $\Delta\boldsymbol{x_k}$ by the inverse masses $w_k=1/m_k$

Finally, we have

and

Following is a concrete example in one-spring case.

What about multiple springs? The Gauss-Seidel approach projects each spring sequentially in a certain order. We cannot ensure the satisfaction of every constraint. But the more iterations we use, the better those constraints are satisfied.

Instead of representing all the elastic potential as a penalty-like term, PBD treats all the constraints as infinitely hard constraints. For instance, in a mass-spring system simulation, PBD also treats all springs as inextensible springs.

Projective Dynamics #

For simplification, we will denote $\boldsymbol{y}:=\boldsymbol{x^{t}}+\Delta t\boldsymbol{v^{t}}$, so

The optimization problem offers an interesting insight into Position Based Dynamics. If we define an energy potential EP BD as the sum of squares of the PBD constraints, we notice that the PBD constraint projection solver attempts to minimize F using a Gauss-Seidel-like method. Unfortunately, the PBD solver does not take the inertial term explicitly into account, and further the projection of each individual constraint is not guaranteed to decrease energy because the effect of the remaining constraints is ignored. Nevertheless, PBD is typically quite effective at reducing F and can be therefore understood as a heuristic variant of the implicit Euler method.

Mass-spring System #

We wll start introducing PD from simple mass-spring system containing $n$ vertices and $m$ springs.

In the first step, we wll reformulate elastic energy

where $\boldsymbol{x_1,x_2}$ are spring endpoints, $r$ is the rest length, and $k$ is the spring stiffness. The key to reformulation is the following fact showing that the spring potential is a solution to a specially designed constrained minimization problem.

That is, from

we goes to

where $i$ indicates the i-th spring and $i_{1},i_{2}$ indicates two end points of the i-th spring.

It might sound like a bad idea to reformulate our original unconstrained minimization problem into a constrained optimization problem with an extra auxilary variable $\boldsymbol{d_i}$ which almost doubles the degrees of the freedom of the system. However, the new problem can be really easily solved by alternating optimization (also known as local-global optimization). In the local step, we fix $\boldsymbol{x}$ and find $\boldsymbol{d}$. In the global step, we fix $\boldsymbol{d}$ and find $\boldsymbol{x}$ by solving a large convex quadratic optimization problem iteratively.

where

where $G_i\in R^{n\times1}$ is the incidence vector of the i-th spring, i.e. $G_i(i_1)=1$ and $G_i(i_2)=-1$ and zero otherwise. Similarly, $S_i\in R^{m\times1}$ is a indicator vector for the i-th spring, $S_i=\delta_i$. Based on this, we have

Notice here we only discuss one dimension (x/y/z). If considering 3 dimensions, we will need to do $L^3=L\otimes I_{33},G^3=G\otimes I_{33}$. $\otimes$ denotes Kronecker product. It is easy to see both $M/\Delta t^2$ and $L$ are symmetric and positive semi-definite.

ARAP Material #

We start from how to calculate deformation gradient $F$ for a tetrahedron with 4 vertices $\boldsymbol{x_1,x_2,x_3,x_4}$ in deformed state and $\boldsymbol{X_1,X_2,X_3,X_4}$ in reference state . We know that

So $F=D_sD_m^{-1}$ where $D_s=[\boldsymbol{x_1-x_4};\boldsymbol{x_2-x_4};\boldsymbol{x_3-x_4}]$ and $D_m=[\boldsymbol{X_1-X_4};\boldsymbol{X_2-X_4};\boldsymbol{X_3-X_4}]$.

For an object containing $m$ tetrahedrons and $n$ vertices, we denote $F_i^T=D_{mi}^{-T}D_{si}^T$ for i-th tetrahedron where $D_{si}^T\in R^{3\times3}$ can be represented as a linear combination of the position vector $\boldsymbol{x}$ (note here $\boldsymbol{x} \in R^{n\times 3}$ is different from above in $R^{3n\times 1}$) as $D_{si}^T=A_i\boldsymbol{x}$ where

where the non-zero columns are corresponding to the indices of the first second third and fourth vertex of the i-th element. If we denote $G_i^T=D_{mi}^{-T}A_i$ then $F_i^T=G_i^T\boldsymbol{x}$.

Consider we now define potential density as

where $F=RS$ is the polar decomposition of $F$. So the potential is equivalent with ARAP. Then

$p_i$ can also be represented with selector matrix $p_i=S_i^T\boldsymbol{p}$ where $\boldsymbol{p}\in R^{3m\times3}$ is the stacked version of $p_i$. And $S_i=\delta_i\otimes I_{33}\in R^{3m\times 3}$.

So the energy can be rewritten in

Gradient Descent #