I encountered some problems when implementing the cloth simulation algorithm from Baraff & Witkin 98's [Large Steps in Cloth Simulation](http://www.cs.cmu.edu/~baraff/papers/sig98.pdf).


----------


**Baraff & Witkin 98**

Consider the cloth as a particle system in 3D, which consists of N particles.

1. Equation of motion

In each time step, the implicit integration of particle system is:

$$
\begin{pmatrix} \Delta \mathbf x \\ \Delta \mathbf v \\ \end{pmatrix} = h\begin{pmatrix} \mathbf v_0 + \Delta \mathbf v \\ \mathbf M^{-1} \mathbf f(\mathbf x_0 + \Delta \mathbf x, \mathbf v_0 + \Delta \mathbf v) \\ \end{pmatrix}
$$

where $h$ is the time step (scalar)


2. The partial derivative equation

Take the 1st order Taylor expansion of $\mathbf f(\mathbf x_0 + \Delta \mathbf x, \mathbf v_0 + \Delta \mathbf v)$, and introduce constraint matrix $\mathbf S$ (make mass matrix from $\mathbf M^{-1}$ to $\mathbf M^{-1}\mathbf S$), we have: 

$$
(\mathbf I - h\mathbf M^{-1} \mathbf S\frac{\partial \mathbf f}{\partial \mathbf v} - h^2\mathbf M^{-1} \mathbf S\frac{\partial \mathbf f}{\partial \mathbf v})\Delta \mathbf v = h\mathbf M^{-1} \mathbf S(\mathbf f_0 + h\frac{\partial \mathbf f}{\partial \mathbf x} \mathbf v_0) + \mathbf z
$$ 

Where all of the matrices are treated as $N\times N $ block matrices, eack block sizes $3\times 3 $; vectors as N block vector, each sized $3\times 3 $ (with N number of 3D forces/positions/velocities).

 - $\mathbf I$ is identical matrix

 - $ \mathbf M = diag(\mathbf M_1, \mathbf M_2, ..., \mathbf M_N) $, $\mathbf M_i = diag(m_i, m_i, m_i)$, mass of each particle

 - $ \mathbf S = diag(\mathbf S_1, \mathbf S_2, ..., \mathbf S_N) $, 
$  \mathbf S_i = \begin{cases}  \mathbf I,  & \text{$ndof(i) = 3$} \\ \mathbf I - \mathbf p_i \mathbf p_i^T,  & \text{$ndof(i) = 2$} \\ \mathbf I - \mathbf p_i \mathbf p_i^T - \mathbf q_i \mathbf q_i^T,  & \text{$ndof(i) = 1$} \\ \mathbf 0,  & \text{$ndof(i) = 0$} \end{cases}$, constraint matrix (Baraff98 chapter 5.1)

 - $\frac{\partial \mathbf f}{\partial \mathbf v}$ and $\frac{\partial \mathbf f}{\partial \mathbf x}$ are symmetrical matrices

.

  
3. solve the PDE

The paper solves the aforementioned PDE by **Modified Preconditioning Conjugate Gradient Method** (Baraff98 chapter 5.2) for $\Delta \mathbf v$, then update the velocity $\mathbf v$ and position $\mathbf x$ of each particle.


----------


**My Questions**

Solve $\mathbf A\mathbf x=\mathbf b$ with PCG method requires the matrix $\mathbf A$ to be symmetrical and positive-definite, where according to the aforementioned PDE, $\mathbf A = (\mathbf I - h\mathbf M^{-1} \mathbf S\frac{\partial \mathbf f}{\partial \mathbf v} - h^2\mathbf M^{-1} \mathbf S\frac{\partial \mathbf f}{\partial \mathbf v})$.

1. symmentrical

Treat $\mathbf A$ as an $N\times N$ block matrix, $\mathbf M^{-1}$ and $\mathbf S$ are both block identity matrices, makes $\mathbf A$ a block symmetrical matrix. While $\mathbf A$ as $3N\times 3N$ matrix, $\mathbf S$ is not symmetrical, makes $\mathbf A$ not symmetrical. 

**Question 1**: should I treat $\mathbf A\mathbf x=\mathbf b$ as a block matrix rather than normal matrix? If not, how to make $\mathbf A$ symmetrical?

I then modified the PDE as $\mathbf A^T\mathbf A\mathbf x=\mathbf A^T\mathbf b$, $\mathbf A^T\mathbf A$ is symmetrical.

**Question 2**: How to make the `fliter` procedure (Baraff98 chapter 5.3) compatible with $\mathbf A^T\mathbf A\mathbf x=\mathbf A^T\mathbf b$?

2. positive-definite

The blocks $\mathbf S_i$ in $\mathbf S$ may become zero block when $ndof(i) = 0$, which makes $\mathbf A^T\mathbf A$ a positive-semidefinite matrix, according to Baraff98 chapter 5.3, 

> The CG method (technically, the preconditioned CG method) takes a symmetric **positive semi-definite** matrix $\mathbf A$, a **symmetric positive definite** preconditioning matrix $\mathbf P$ of the same dimension as $\mathbf A$ ,a vector $\mathbf b$ and iteratively solves $\mathbf A \Delta \mathbf v = \mathbf b$.

**Question 3**: How to find the symmetrical and positive-definite preconditioning matrix $\mathbf P$, with $\mathbf A^T\mathbf A$ symmetrical and semi PD?