I encountered some problems when implementing the cloth simulation algorithm from Baraff & Witkin 98's Large Steps in Cloth Simulation.

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)

  1. The partial derivative equation

Rewrite the aforementioned equality of motion into $\mathbf A \Delta \mathbf v = \mathbf b$ form, ( 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$, which 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 x})\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$ (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}$, $\mathbf p$ and $\mathbf q$ are two oithogonal constraint direction, $\mathbf S$ is constraint matrix (Baraff98 chapter 5.1)

  • $\frac{\partial \mathbf f}{\partial \mathbf v}$ and $\frac{\partial \mathbf f}{\partial \mathbf x}$ are force derivatives, which are symmetrical matrices


  1. 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 x})$.

  1. symmentrical

In block view: treat $\mathbf A$ as an $N\times N$ block matrix, $\mathbf M^{-1}$ and $\mathbf S$ are both block diagonal matrices, makes $\mathbf A$ a block symmetrical matrix. While in normal view: $\mathbf A$ as $3N\times 3N$ matrix, $\mathbf S$ is symmetrical, but product of $\mathbf M^{-1} \mathbf S$ and $\frac{\partial \mathbf f}{\partial \mathbf v}$ (or $\frac{\partial \mathbf f}{\partial \mathbf x}$) 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$?

  1. 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. To apply PCG method, 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?

Some docs I also refered:


1 Answer 1


Baraff and Witkin propose to incorporate constraints by modifying the linear system by the constraints matrices $\mathbf S_i$. As they state in the beginning of Section 5.2, the resulting system is not symmetric anymore. Therefore, the modified linear system must be computed with a linear solver that can treat non-symmetric systems (as - if I got you correctly - you proposed with the normal matrix $\mathbf A^T \mathbf A$). However, Baraff and Witkin introduce the modified conjugate gradient method, which is a variant of the original conjugate gradient method, that computes a solution of the initial linear system with considering the constraints at the same time, without the need to take care of the non-symmetry. In their paper they do not solve the normal equations ($\mathbf A^T \mathbf A$), but rather the original system $\mathbf A \Delta \mathbf v = \mathbf b$, so your questions ask something which is beyond of the paper.

  • 1
    $\begingroup$ Follow your answer, I found out that the Modified Conjugate Gradient method is proposed to deal with the non-symmetric semi-positive-definite linear equation(section 5.2). Thanks a lot. $\endgroup$ Commented Feb 16, 2016 at 1:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.