Real-Time Physics Simulation Forum

Re: Baraff/Witkin's integration matrix not positive: what to do?

2020-09-17T02:22:21+00:00

Finally, it took 11 years but it was worth waiting: a new paper by Theodore Kim addresses the problem, and one of the references points to this forum thread. See http://www.tkim.graphics/FEMBW or https://twitter.com/_TheodoreKim/status ... 0997393409

Statistics: Posted by Erwin Coumans — Thu Sep 17, 2020 2:22 am

Re: Baraff/Witkin's integration matrix not positive: what to do?

2009-12-04T16:22:10+00:00

Well, it is not PD only for large time steps and large stiffness/damping constants, but as their paper is entitled "large steps in cloth simulation", they should have realized this

I realized it already with my second test of my implementation where the simulation just exploded. After some hours of further testing, I discovered the non-PD-ness of the matrix.

Statistics: Posted by mathis — Fri Dec 04, 2009 4:22 pm

Re: Baraff/Witkin's integration matrix not positive: what to do?

2009-12-01T06:44:08+00:00

Oh, I now realize they are the same, so this mean the matrix in baraff 's original paper of stretch enegy is not a PD matrix ..strange

Statistics: Posted by fishboy82 — Tue Dec 01, 2009 6:44 am

Re: Baraff/Witkin's integration matrix not positive: what to do?

2009-11-30T09:19:25+00:00

The formulas are the same for a uv-aligned triangle. Let (x0, x1, x2) be such a triangle with the uv-coordinates (0, 0), (L1, 0), (0, L2). L1 and L2 are the rest lengths of the edges x0-x1 and x0-x2, respectively. This fulfills the requirement of the paper that |dw/du| = |dw/dv| = 1 if w is the rest state w(u,v) = (u,v,0), where w in general maps (u,v) into the current 3d vertex position.

Now for the current state w, the paper approximates

(dw/du, dw/dv) = (x1-x0, x2-x0) * ((u1-u0, u2-u0), (v1-v0, v2-v0))^-1

In our case, the matrix ((u1-u0, u2-u0), (v1-v0, v2-v0)) is a diagonal with the entries L1, L2, so the formula simplifies to

dw/du = (x1 - x0) / L1
dw/dv = (x2 - x0) / L2

Now the paper uses C(x0,x1,x2) = a * (|dw/du| - bu, |dw/dv| - bv) as constraint and k/2 C^T C as energy. This leads to the energy

E = a^2 k / 2 * ( ((x1 - x0)/L1 - bu)^2 + ((x2 - x0)/L2 - bv)^2 )

which is the sum of the energy of two springs with rest lengths bu*L1, bv*L2 and spring constants a^2 k / L1^2, a^2 k / L2^2. And this energy leads to a non-PD system matrix for large time steps and extensions less than the rest lengths.

Statistics: Posted by mathis — Mon Nov 30, 2009 9:19 am

Re: Baraff/Witkin's integration matrix not positive: what to do?

2009-11-30T08:46:20+00:00

Hi mathis :
Have you ever tested to derive the derivation from Baraff's stretch enegry directly but not use your spring form? Because I think there is still a little difference beacasue baraff use the particle's position as varaible rather than the extension of the spring , could this induce the K matrix to be different of yours,just my guess .

Statistics: Posted by fishboy82 — Mon Nov 30, 2009 8:46 am

Re: Baraff/Witkin's integration matrix not positive: what to do?

2009-11-30T08:20:16+00:00

Yes, they only use stretch/bend/shear-terms, but the stretch-term is a spring. They have a more general setting with triangles and uv-coordinates so that the behaviour is not so dependent on the tesselation, but for a triangle with one u- and one v-parallel edge, the formula for the stretch-term is exactly that of a spring (or actually two springs, one for the u-edge and one for the v-edge, but that doesn't matter).

Statistics: Posted by mathis — Mon Nov 30, 2009 8:20 am

Re: Baraff/Witkin's integration matrix not positive: what to do?

2009-11-28T16:28:44+00:00

Is it possible that they don't use the spring-like energy . I remeber they only use stretch ,bend and shear energy. May be these energy will cause a PD matrix but not all energy you derive will be a PD matrix?

Statistics: Posted by fishboy82 — Sat Nov 28, 2009 4:28 pm

Re: Baraff/Witkin's integration matrix not positive: what to do?

2009-11-28T16:20:23+00:00

They state: "The matrices we ultimately hand to our CG method are positive definite." and "Let us define the symmetric positive definite matrix A by A = M - h df/dv - h^2 df/dx."

If you read mathematical papers on the CG method, there is the requirement of positivity. Maybe that you are lucky and CG converges also for some special situations with indefinite matrices, but that's not guaranteed. Theoretically, if A is indefinite for large h, you can also find a h where A is not invertible (the h where A switches from positive definite to indefinite). But with CG we want to solve Ax=b, so we need an invertible A.

Statistics: Posted by mathis — Sat Nov 28, 2009 4:20 pm

Re: Baraff/Witkin's integration matrix not positive: what to do?

2009-11-28T08:52:20+00:00

Hi:
I roughly read the paper, but it seems baraff only mentioned that the K matrix(d2E/dx2 )
is a symmetric Matrix not a PD matrix (cha 4.1).
and in 6.1 he mentioned that CG methods, however, require symmetric matrices and there is also a label "In fact, they work best on positive definite symmetric matrices", I don't knwon CG Method so I suppose that CG should work not only with a PD matrix but symmetric is enough?

Statistics: Posted by fishboy82 — Sat Nov 28, 2009 8:52 am

Re: Baraff/Witkin's integration matrix not positive: what to do?

2009-11-26T10:14:13+00:00

if you read the cloth papers people always state that the matrix is PD.

Yes, this is also my impression. These papers seem to be a bit optimistic.

The case of non-PD is not academic, it occurs already for typical time steps together with stiff cloth. For positive definiteness, you get the condition m + h^2 k (1 - L/|x|) > 0, so for |x| < L:

h^2 < m / (k * (L/|x| - 1))

This is the typical relationship between stiffness and time step which is already known from explicit euler integration, but with an additional division by (L/|x| - 1). Now one could argue that for stiff springs, |x| is always close to L, so the division is by a small number and large time steps are possible. But in a test with damped stiff structural and bending springs (bending springs being based on angles), I had to use time steps of about 1/1000 second.

Statistics: Posted by mathis — Thu Nov 26, 2009 10:14 am