Stan's cloth simulation demo/discussion

[Erwin Coumans]

(moved reply to general discussion, see original topic here: http://continuousphysics.com/Bullet/php ... .php?t=401)

1) baraff & witkin large steps in cloth simulation SIGGRAPH 98
2) SIGGRAPH 2005 physics presentation.
3) Jeff Laner's demo and gamedev/gamasutra article (good starting point for anybody)
4) I just recently found this webpage for a UW phd student who posted an implemention (with source) of a backward integrator as a maya plugin

Thanks for the great looking cloth demo and explanation!
Those 4 links are missing, can you provide them, especially the last one? The first one is http://www.cs.cmu.edu/~baraff/papers/index.html

- Sequential impulses (SI) can be used to simulate cloth. SI can be used to simulate oscillating pendulums.

Would you drop the angular terms and just add particles (point rigid bodies) with point to point constraints (with a distance) using the commonly use velocity based 'relaxation' scheme like PGS / SI?

Also, Stan, please add this link to a nice cloth demo including source code(Verlet based if I remember correctly):
http://www.maxgarber.com/projects/cloth/

[asmithers]

I tried supporting both baraff and verlet simulatenously

I started off using a baraff based approach to my cloth sim, then moved to a modified verlet. Some of the effects I made for baraff, namely rigidity came across quite nicely to the verlet model, although the data wasn't valid for both models. My baraff model always ended up exploding after a few minutes. The energy would never settle.

I do like stan's implementation, there is something pleasing with this implementation. Its extremely similar to my own in that it shows bows when it collides to the floor, this being a bi-product of my attempt to implement choi's algorithm for buckling. (choi02).

when I get time ( dont know when at moment ), I'll go through this algorthim in more detail, it looks promising based on the results. Although seems to suffer issues with rapid change of positions... something that was tough for me to solve when attaching to characters.

Have you tried working on moving cloth at large speed intervals, Ie a flag on presidential car for example ?

andi

[melax]

Hey, thanks for checking it out. Yes I'll definitely be improving the writeup and adding the suggested links. I' should probably add a link back to this bullet thread as well.

About the car on the president's flag - Im imagining that the frequency of the flag flapping would be very high. So it might be better to use a lower cost simulation but use more timesteps.

Agreed, rapid changes in position and velocity are a challenge for simulating cloth. It simply takes iterations for forces to propogate across the mesh. Both Verlet/constraint-iteration and backward integration handle these much better than forward euler since forward euler cant really do stiff systems. If an iteration scheme can rearrange the iteration order and make some assumptions as to which end of the cloth is "attached" and which is free, then it could have an additional "responsiveness" advantage here.

I'm on vacation for a couple weeks. During this time, I might not be able to reply to any questions in a timely manner.

thanks

[mewert]

Impressive stability, I couldn't make it explode.

We use a verlet/relaxation-based cloth simulator, which has been working pretty well for a reasonably high poly mesh on a fast moving character. What Erin is calling Sequential Position something is what I like to call a relaxation solver. The reason it is often coupled with verlet integrator is so that momentum is indeed changed. Since in verlet integration, the velocity is represented by the differnece between the last position and the current position. So as you update the postion of a particle, you change it's velocity at the same time ( assumes all particles have the same mass ). This also is the reason it's semi-implicit. Since the velocity is always in-sync with the results of your solver wrt time. Forward integrators are unstable because the results of the solver are a step out of sync with the application of those results ( one way to look at it ).

We had a back-euler cloth solver that worked for square patches. Once we started applying the cloth to a shirt in continuous colliding contact with a moving character, things started going awry. Stability was a big issue and often the CG solver would refuse to converge. I think we had some problems with our math, but by that point the verlet/relaxation cloth was in a good state and we dropped the back-euler cloth.

My main problem with the verlet/relaxation cloth is that it is very difficult to tune. Artists could not do it unsupervised and it has a very limited range as far as material types it can simulate. Also, collision detection quickly becomes the dominant factor for CPU time, since you must do CD more than once a frame per vertex. Back-euler is O(vertex^2) for the CG solver, but the cost of collision detection is much cheaper so maybe it could be faster compared to the O( vertex ) relaxation solver.

Anyways, would be great to see Stan's cloth working on an animated character. Man, are we ever a demanding bunch, eh?

[/img]

[Oscar Civit Flores]

Hi all,

There was a nice presentation on "verlet clothes on characters" at GDC 2005 from the guy who added cloth physics to the Source engine for Vampire: The Masquerade:

https://www.cmpevents.com/sessions/GD/P ... ePants.pps

It's a very good presentation and the ideas there actually *work*.

I've been implementing a similar system lately for an upcoming PC game with some modifications to handle fast moving characters with attached clothes (it's a soccer game, so players move quite fast) and it works very well. Changing relaxation order to handle attached points first and biasing it to "move more free particles than attached ones" effectively speeds up convergence.

The main problem I've found are cloth-character collisions, which interact a lot with the length constraints. I have to process them several times per frame to obtain good results, and, also important, perform a final collision treatment step to ensure that, even though the length/curvature constraints may be slightly violated, the cloth patch is not intersecting the character.

Also, as somebody else pointed out, establishing the material properties in a relaxation based system is a bit tricky. I've tried to obtain parameters independent of the timestep and the relaxation iteration count and map them into the range [0..1] to make them intuitive, but currently this is still an experiment.

Oscar

[Genscher]

Hello!

That is really a nice implementation and cloth demo. I also liked the very good structured code which is really good for understanding.
I had also the luck to get the old "freecloth" implementation of Baraf & Witkins paper, but is was really confusing and not helping.

Since I am one of the developers of the (hopefully) comming cloth implementation in blender
[www.blender.org - our cloth patch is based on Provot: http://mediawiki.blender.org/index.php/ ... Simulation ]
, I am really interested in more stable solvers. explicit / middle point does fine for simple simulations and looks even more convinceable (our cloth patch is actually used for some TV show about pirate plags in the USA because they said "your cloth looks better than Maya ones" ) but it has defineatly it's downsides when used on rigged boddies. Then it explodes. Another implicit solver could change this behaviour

For our collision detection, which you mentioned as further work in your cloth demo also, I started to implement oriented inflated bounding volume hierarchies ( http://www.gris.uni-tuebingen.de/public ... proved.pdf ).
This works very nice!
Our C-code implementation (the oriented inflated BVH) is also free available under GPL on the side i mentioned above.
We doesn't tested the cones for self-collision yet.

What is missing from the paper is the part about collision response where your cloth demo comes in. They suggest to use Baraf & Witkins constraint method and they also improve it in some ways. Some newer siggraph presentation combined all the papers to one nice cloth simulation.

If you don't mind I would really like to use your code as reference for a cloth solver implementation in blender.

Actually my dream is to combine bullet's rigid bodies, blender armature and cloth simulation into one thing, but that's a long way to go

Greetings and great work!
Genscher

[melax]

...
If you don't mind I would really like to use your code as reference for a cloth solver implementation in blender. ...

Hi, and thanks for the kind words. You are certainly most welcome to use or reference the code as you see fit.

[Genscher]

Good morning!

Thank you for your fast answer!

My first try now is to use our spring network with your solver. That's easy because we have the same spring-definition and building function

In the moment I'm guessing the type/class definitions of "Array" and "float3x3" which are used in vec3n.cpp.

I would appreciate if there would be a possibility to take a look at the header files ("vecmath.h" and "array.h") so I know how to code the missing classes/typedefs

But thank you again for this great reference!


-----------
eMail is genscher22 at gmx dot net

[Genscher]

Hello!

I worked on integrate melax solver concept and finished it's integration
I have cloth flying in the wind based on melax solver in blender

The problem is, that i have to convert between different storage (C++ to C interface) so the speed of the simulation is going down a bit.
That's why I recently took a look into "Speeding up cloth simulation"/"Decomposing cloth" from Boxman/Asher where they introduce a modified PCG (preconditioned conjugate gradient) algorithm. I also took a look at the "high performance solver" introduced by Hauth et. al.

I implemented a function to create P and P^-1.
Now I am missing the vector z.
Anybody here who want to help building a precondiitoned CG algorithm

Greetings
Genscher

[melax]

Congrats on the blender implementation, I'm glad it all came together.

So far, I have just implemented a modified conjugate gradient solver. The preconditioned cg solver is explained really well in: SHEWCHUK, J. 1994. "An introduction to the conjugate gradient method without the agonizing pain"

Wrapping your head around the modified algorithm with preconditioning is a bit trickier. Baraff/Witkin's paper does provide the pseudocode.

------------
This isn't a big issue for me yet, but since people asked about it...
On the subject of cloth responsiveness...
A straightforward cloth implementation is going to require numerous iterations for a tug on one end to propogate across a spring network with a lot of points. Rather than spending cycles updating parts of the cloth that dont need updating, it only makes sense to try and focus the computation where it is needed.
On of the mentioned power point presentations describes a verlet/relaxation implementation that would order the points and bias the iteration based on where the cloth was attached to a moving body. It may be possible to do something similar for implicit integration. Recall that the conjugate gradient starts with an "initial guess" and then iterates from there to the result. Usually the zero vector is used as a starting point. This could be a good place to insert the extra information that we know about the simulation. i.e. make a better "initial guess" for the delta velocity vector (dV). Then after that you just use the same old conjugate gradient method. This strategy should give you the best of both worlds. You take advantage of "additional knowledge" about cloth construction and attached/manipulated vertices, meanwhile your iteration solver still uses the same algorithm (regardless of the cloth configuration) that can be run in PARALLEL - especially important for any non-single-cpu implementations.

So does this work? I did some (limited) experiments initializing based on how the cloth was tugged by the mouse. I did find a downside. Normally, since I only care about aesthetic behavior, I use a rather high epsilon value for the CG solver to keep the iteration count really low. When it always uses zero as the starting vector, this doesn't matter - but with the starting point being dramatically different from frame to frame, there were small but noticable discrepancies in the behavior at points on the cloth.
Note that my initialization of dV was poor - everything was initialized to the "tug". Ideally you would only want to apply that to things at full stretch starting from the attached point. This would be similar to doing a single one of those verlet/relaxation (loop-dependant) iteration type of steps to initialize dV. In essence, for my tests, I was choosing what would have been a worse "starting point" than zero for certain parts of the mesh. In conclusion, I cant really make any conclusions at this point.
-------------------

In response to the request for a demo with cloth on an animated character... Yes, I agree that is something that I need to do. Unless all that "pretty math" can lead to usable applications with good performance, then it doesn't do much good. Trouble is finding time to work on this


Thanks for the comments all

[Genscher]

[edited because one error was found, but the problem is not solved]
[update in the post below]

[Genscher]

Hello!


Man, sometimes, it's good to have a little break and then going back and watch over the problem again.

I have 2nd-order BDF working with melax source
It's really just as plain as in the paper

Greetings!

Actually, I was a bit too fast with posting this

My BDF implementation seems to behave a little bit weird.
I have some simple cloth hanging down (-9.6 gravity on z-axis) so it should just hang down along the z-axis.
But when using my BDF formula the cloth stabilize in a state with arround ~15 degree away from the z-axis.

See this compare:


[But even in this state you can also see the nice wrinkling behaviour of the BDF in compare to the implicit euler (BE).]

Here's my current BDF-implementation:

Code: Select all

A -= (dFdV * (2*dt/3) +  dFdX * ((4*dt*dt)/9)); // A = I - (dFdV * (2*dt/3) +  dFdX * (4*dt*dt/9))
		
B = (X-Xold) * (1/3); // (X - Xold) * 1/3 
		
temp = (V*8 - Vold*2) * (dt/9); // + (8*V - 2*Vold)*(dt/9);
B = B + temp; 

temp = dFdV * V;  // + ((4*dt*dt)/9)*(F - dFdV*V)
temp = (F - temp) * ((4*dt*dt)/9);
B = B + temp; 

temp = X - Xold;  // - (2*dt/9)*dFdV*(X-Xold);
temp = dFdV*temp;
B = B - temp*(2*dt/9); 	

and the update code:

Code: Select all

Xnew = X + dV;
		
Vnew = Xnew*(3/2) - X*2;
Vnew = Vnew*1 + Xold*(1/2);
Vnew = Vnew * dt;

Maybe someone knows where the error is

[Genscher]

Hello!

In the meantime i added the AIMEX ("adaptive implicit-explicit") scheme to melax source which gives arround 5%-15% speedup.

Just add this code to PreSolveSpring() just before you compute the jacobians:

Code: Select all

/////////////////////////////////
// adaptive implicit explicit scheme 
// from "Speeding up cloth simulation", Boxman & Ascher 2001
////////////////////////////////
	
// handle every bending springs explicit
if(s.type==SPRING_BEND)
	return;

// check the stability equation 
if((dt * (k * dt + 2 *spring_damp))<=0.2)
	return;
////////////////////////////////	

For all interested why this makes the simulation faster and why that works:
The basic of the idea is to handle special types of springs explicit instead of implicit. That's also the reason why we need a stability criterium because our simulation could explode again and we would have lost one of the reasons to go implicit.
Speed:
Explicit solving saves computing time since we don't have to compute the forces for the jacobians.
Another nice effect is, that the matrix A gets more sparse and the CG algorithm is faster in result.

Greetings
Genscher

[melax]

Hey, I hadn't heard of the implicit-explicit idea to just use forward for the weaker forces. Good idea.

I'm sure you noticed this, but just in case you didn't:
Note that when the cloth's springs are initially constructed it creates an entry in the matrices. So you might want to avoid this creation step for any explicit springs. Otherwise, it might still be working through those extra (now zero) matrices each time it does a matrix multiply.


--------
For further optimization, you can take advantage of the symmetry. Specifically there are two useful things to do here. Try only storing one of the 3x3 blocks for anything off the main diagonal. Instead of using a general float3x3 class for each of the blocks, use a symmetric3x3 that stores 6 numbers (float3 diagonal, and float3 offdiagonal). BTW the correct choice is for symmetricmatrix.offdiag[0] to refer to myfloat3x3[1][2] and myfloat3x3[2][1].

Its a straightforward optimization.

One of my objectives was 'Education' - to make all this backwards integration understandable by the typical developer. That becomes more difficult if the reference implementation includes everything and every optimization.

--------

BTW gravity is 9.81 not 9.6 i'm sure that just a typo, but couldn't resist. I usually just use 10.


-------------------
I didn't use a template metaprogramming approach to optimize the expressions. Although its a fun mental exercise to implement, the approach is very obfuscated and only yields a performance gain in release mode. The Optimized Expressions approach, lets you optimize large equations without making the higher level code look ugly. However, i only went as far as what I needed for the formulas I was using. If you tried to write a longer expression then it wouldn't compile for you. (I noticed your use of temporaries). If you want to put the full expression a a single line of code, then the big vector header file just needs to "grow a bit" by hand (instead of using template compiler tricks) to accomodate what you want.

I certainly hope there isn't some sort of mistake in the big vector implementation that is causing your problem.

--------------

I wish I had more time right now to be working on this.

I certainly appreciate your feedback.

[Genscher]

Hello!

First of all: Thank you for your advises! I will take a look at them.

I tried to debug the BDF formula today and was totally confused after I found the error:

Code: Select all

float a = 1 / 3; // result: a = 0!

So I only had to convert all integers to float and the whole BDF formula is working nicely! What a stupid error

I also have to say that all my "improvements" are not to show any faults/slowness in your code, it's just my attempt to understand the physics/maths which lays behind your code

I'm also happy that you didn't used templates so much because there is mostly a problem when it's going to "crossplattform" (especially GCC versus MSVC-compiler problems).


Greetings
Genscher

update: I just tested your advice to leave the AddBlocks() call for bending springs. Almost 50% speedup in compare to the fully implicit method