Impulse-based dynamic simulation library (IBDS)

[mewert]

Ah, you are right. How on earth can they call it a "Fast" method when they are solving the linear system exactly?

Well, fast is a relative term and it's in an academic paper. So I would say "fast" is use in the same way I've seen it used before. I hate being a snob like this, I've got enough of that attitude from those with higher education than myself

Still, they are using a standard method and putting their own label on it. It's just a Newton solver for constrained minimization.

Well, it seems that's what most graphics and physics papers do, apply known methods in a new way to problems in those domains. I do kind of find it annoying that I have to read through lots of convoluted math lingo just to arrive at the conclusion that it could have been summarized quite simply by saying "we are applying known method X to problem Y with the modifications Z". And flesh out Z a bit more.

I think this paper makes some valuable contributions and has some interesting and novel ideas for simulating cloth. I hope it works well in a real-time, extreme strain situation.

[mewert]

a) They might use more sophisticated integration methods to find candidate positions

Nope, doesn't look like it. I like the idea of treating the diagonals of the quadrilateral mesh as springs, though.

b) The use a Newton method as opposed to NGS where they solve the "inner" LS with maybe a CG method

Nope, they use a direct solver: PARDISO. I think they would benifit from using GS with an accumulation vector, though. They also use a non-linear mass matrix ( using triangle area to calcluate the mass ).

It don't quite understand where they handle collision then? Is this another phase after integration and constraint force correction?

They state they treat their solver as a velocity filter. Which, would mean you could do collision resolution as a seperate phase. Potentially having another outer-loop to the solver. If I implemented this I would use a GS solver with simple collision-plane constraints mixed in with the strain-limiting constraints, maybe even throw in some friction for good measure

And the final question to bring this thread again back to the original discussion.

Yeah, we've kind of coopted this thread ( or I have, rather ).

Using position-based constraints in this case is a method of strain-limiting ,or enforcing inextensiblility. Momentum is conserved ( energy doesn't have to be conserved ) by applying the projected position change as a velocity update. Which is a bit different than a pure position based dynamics method, where the velocity is implicitly encoded in the positions. That difference is what I like about this approach.

PBD: the solver outputs X_i+i directly
Inextensible cloth: output dV that was arrived at by finding X_i+1, which is further refined to a specified tolerance.

Almost back on topic!

[Dirk Gregorius]

Nope, doesn't look like it. I like the idea of treating the diagonals of the quadrilateral mesh as springs, though.

As Erin pointed out correctly bending and shear forces are summarized in the potential forces. Since those are not handled as constraints they are modeled as spring and dealt with in the integrator. The term potential reminds me on the Baraff paper. The idea of treating diagonals as springs came from Erin. It makes sense to model those as soft springs and here is actually a reason for symplectic Euler since it handles this stuff better than the parabolic integrator.

Nope, they use a direct solver: PARDISO. I think they would benifit from using GS with an accumulation vector, though. They also use a non-linear mass matrix ( using triangle area to calcluate the mass ).

That's what I meant. They solve with all constraints simulatineously instead individually as we do. I would have chosed a CG solver for the sparse system as "direct" solver first, but maybe other methods are better suited. I am not an expert on this topic. What do you mean with non-linear mass matrix? I understand that you can assign e.g. 1/3 area times "density "of each adjacent triangle to a vertex. Do you mean this or do they change the mass matrix as the system deforms?

Which is a bit different than a pure position based dynamics method, where the velocity is implicitly encoded in the positions. That difference is what I like about this approach.

Isn't this exactly what the "Position Based Dynamics" paper does?


So here is how I understand the paper:

1) Integrate system forward considering only gravity, bending and shear forces (later are springs) to build potential velocities and positions
2) Apply the change in velocities found by their method and correct the positions.
3) Handle collisions somehow

So what phases am I missing?


I go now and read the paper again, since I couldn't find all this stuff in it. For me this was more "much ado about nothing" :-/

[mewert]

it's more like this ( I added the collision phases ):

1) project positions forward to the end of the time-step as a candidate new position
2) collection collision constraints based on projected positions
3) enforce constraints using lagrange multiplier solution to manifold projection. This the variables you are solving for are position, but you output an impulse. Do this wrapped in non-linear iteration loop.
4) apply constraint impulse
5) integrate forward to new positions


How does this differ from PBD?

In PBD, step 4 and 5 don't really exist. velocity is implicitly Xi+1 - Xi.

In step 3, PBD actually modifies the jacobian as you are solving for the free variables. With inexensible cloth your jacobian stays fixed when solving the constraint system ( it can change after you have a full solution, though ). This is the biggest difference.

I see a number of known types of ways of solving cloth ( with the author of, what for me; is the reference paper I would reach for ). I'm sure my list is going to irk someone

1) position-based verlet/relaxation ( jacobsen )
2) implicit cloth ( witkin-baraff )
3) explicit euler ( trival )
4) lagrange multiplier projected-manifold constraints for strain-limiting ( Goldenthal, et al ) Maybe there is an earlier paper on this? Or is it an obvious solution to everyone, or am I wrong and it is just an example of one of the other types?

I know that 4 is kind of obvious is some ways, I've had informal chats about implementing cloth using a rigid-body solver, and have seen examples of it too. But I haven't actually seen it presented in a formal fleshed-out way before and I appreciate all the details presented in the inextensible cloth paper.

[raigan2]

My understanding is that in "implicit velocity" methods like Jakobsen/PBD, steps 4,5 _do_ exist -- they're executed immediately after solving each row of the matrix, Gauss-Seidel style.

I might be wrong, but isn't this analogous to what happens in iterative impulse-based methods? Rather than calculating the constraint forces and then integrating forward to get the new velocities, you immediately change the velocities row-by-row.

So you could say that just as PBD implicitly encodes velocity in positions, iterative impulse solvers implicitly encode acceleration/force in velocities. Of course I might be confused..

[mewert]

Right, you update the state ( position, velocity ) for each row, actually not even when you have finished with one row. You aren't just changing velocities, you are changing the jacobian as well. Since the jacobian is also implicit. c = (x1 - x2)/|x1 - x2| You change x1 and you've changed the jacobian. Which is quite different. I don't know if that is a bad thing but it neccesitates updating your collision constraints frequently.

Non-stationary, is that the mathese for that kind of behavior? Anyone with more math skillz care to comment?

I should expand 2) and 3)

2)...
3)...
3a) do 2 & 3 until you are within error tolerance.

[bone]

Of course the parabolic integrator can't be used for rotation integration. I use a Runge-Kutta of fourth order for this.

Just FYI, there's probably a better method out there if you're looking to accurately integrate rotations including the gyroscopic (coriolis) effects:

http://math.ucsd.edu/~sbuss/ResearchWeb ... /paper.pdf

In a nutshell, Mr. Buss provides both an "augmented" 2nd order method and a 3rd order method, both of which give better results than RK4 at a fraction of the cost.

[ronygold]

A few comments about our paper "Efficient Simulation of Inextensible Cloth" discussed here before.

The main difference of this paper, compared to strain-limiting , PBD, and other methods that use Non-linear Gauss-Seidel (NGS) is that our method solves for all constraints simultaneously. This has been nicely explained above. It may seem like a slower approach to some of you, but if you would try to enforce low strain (1% or less) for a large mesh (several thousands of vertices), and you would iterate NGS until it convergence - it would be much slower than solving all constraints at once. You can check figure 5 in the paper to see the difference in performance (compare fast projection to SL-Gauss-Seidel).

When we called our method 'fast' it is because we compared it to other methods under the same configuration (with the same permissible strain). For these settings, to our best knowledge, our method is indeed the fastest; not because we're academic snubs

Now, if you have linear constraints (for example for collision reaction), or a small number of constraints (low resolution mesh), or you don't care about having high strain - then you can perform only a small number of NGS iterations and you would get faster performance than using our method. These assumptions may be valid for real-time cloth simulation, or other applications, but these are not the assumptions we used in our paper.

As for a different topic:
Fast projection method is NOT equivalent to a Newton solve, as fast projection changes its target function for every iteration (which does not happen during Newton's iterations). So we're not solving the same equation as trying to solve a constrained minimization problem using multiple Newton's iterations.

For bending we used forces -- either quad based forces ( "Stable But Responsive Cloth") or hinge based forces ("Discrete Shells").

For collision treatment we used Bridson's "Robust treatment of collisions, contact and friction for cloth animation".

I hope this answers some of the questions regarding our paper.

Rony

[Dirk Gregorius]

Thanks for pointing this out. Actually once I read through the references (e.g. Hairer) it became much clearer. Anyway, regarding collision handling, this is then just another velocity filter applied after your strain limiting pass?

And in the context of real-time cloth. Your quad mesh is basically as you say warp- or weft aligned. Did you consider solving those "chains" with a direct solver in an iterative GS loop? This would be a block GS then. The chains would have a nice matrix pattern that can be solved very efficient. Using so called "stiff subsolver" is something we use for rigid bodies.

[ronygold]

Yes, collision reaction is another velocity filter pass.

As for solving for chains, we did not consider this approach. Solving for each chain will be very cheap, but still we may need many chains and many iterations. It should reduce the number of GS iterations by an order of magnitude, but each iteration will be more expensive. It may be a very good idea for real-time cloth simulation as the algorithm will not be iterated until convergence and actually solving for a chain in real-time is possible.
For complex topologies there is quite large number of triangles, especially along the seams between different patterns, which makes the decomposition to threads a little more complicated.

BTW, how many vertices it is possible to simulate using current real-time cloth simulators?

[Erin Catto]

Hi Rony,

Thanks for visiting the forum and answering some of our questions.

Our system can simulate roughly 2 vertices per microsecond (3.5 on a dual-core). For a game running at 60Hz, we have 16.6 milliseconds per frame. For cloth we are willing to give about 1ms per frame. So that's roughly 2000-3500 total vertices on screen.

Fast projection method is NOT equivalent to a Newton solve, as fast projection changes its target function for every iteration (which does not happen during Newton's iterations).

I don't understand your assertion. According to this:

http://en.wikipedia.org/wiki/Newton's_method

the most basic Newton method uses a new function value f(x_n) and function derivative f'(x_n) each iteration. Are you changing to a new function each iteration: f1(x_1), f2(x_2), ...?

[Dirk Gregorius]

I agree with Erin's numbers. In Lair we did only simple flags and banners which could be solved very fast and also could be solved in parallel on several SPUs. This was never a bootleneck so I didn't made any real measurements. I think the upper limit of simulataineous solved banners was 200 during large battles where the vertex count varies maybe from 100 - 600. Also banners with more than 500 vertices were the minority.

Yes, collision reaction is another velocity filter pass.

Why do you treat them separately and not in the same relaxation loop? E.g. if I would use a chain subsolver I would relax the chains with the contact constraints together. If you look at this from this perspective I would argue that your stretch solver is just an even "stiffer" subsolver? This is just out of curiousity, your cloth looks really amazing!

Cheers,
-Dirk

[ronygold]

Thanks for the feedback and the input regarding real time cloth simulation!

Why do you treat them separately and not in the same relaxation loop? E.g. if I would use a chain subsolver I would relax the chains with the contact constraints together.

Ideally one would like to solve for both inextensibility and collision reaction together, we tried to do it but found that it is something hard to get (with reasonable performance) and that the velocity filter worked well enough for our needs.

We tried to do this in two ways:
1. to treat collisions as constraints and to add them to the linear system and solve for inextensibility and constraints at once. This did not work very well for several reasons: 1. it was slow because the number of constraints changes for every iteration and we had to reallocate my matrix for every time step/iteration. 2. We encountered some robustness and convergence problems, especially with edge-edge collisions for moving objects. 3. Using this approach, one can get into over constrained or contradicting constraints situations which are hard to deal with.

2. Apply collision reaction after every linear solve. This approach worked, but it slowed the convergence considerably which made the overall performance slower.

the most basic Newton method uses a new function value f(x_n) and function derivative f'(x_n) each iteration. Are you changing to a new function each iteration: f1(x_1), f2(x_2), ...?

You are correct, each 'fast projection' iteration considers a different 'f'. For every such 'f' a single newton iteration is performed and then a new 'f' is computed. The process is repeated until the constraints are met up to our threshold. This is not the same as performing several Newton iterations on the same 'f'. See section 4.3 in our paper for more details about why and how this 'f' changes.

[Dirk Gregorius]

Maybe you allow for one more question regarding cloth setup. You write that you cloth mesh was imported as quat-dominant mesh. I like the idea and espaecially your argumentation about the benefits of this approach. My question is how you find the bending springs? So far I would create one constraint for each edge and IIRC two shear springs for each diagonal in a quad? I should maybe also ask how you handle the seldom triangle faces for completeness.

For our projects I simple used triangle meshes since this is easy for an artist to model. All meshes could be flaged as cloth and would then be handled by the cloth simulation. Creating the physical model is simple from this data. I create one constraint for each edge and a "softer" constraint for each diagonal across a shared edge between two adjacent triangles (like Bridson).

Maybe you can share some of you experiences with cloth authoring? You mention models from OptiTex - are those models available? It might be an interesting experiment to run the models from your paper through a much simpler solver and compare the results...


Cheers,
-Dirk

[ronygold]

My question is how you find the bending springs?

Once inextensibility constraints are applied, the mesh can be triangulated. In fact, the mesh can start as triangular and the constraints can be applied to a subset of the edges. Then, you can use any triangle based bending force.

We used hinge based forces such as "Discrete Shells" (which is the same as the bending force described by Bridson in 2003). Since the mesh deforms almost isometrically, you can also use simpler (but more efficient) hinge based models such as "cubic shells" or "A Quadratic Bending Model for Inextensible Surfaces".

There exist quad based bending forces such as the one used in "Stable but Responsive Cloth" which we also used. For our sparse set of triangles we did not apply bending resistance, we could have done better job here, like using triangular based force for them, but since we saw no visual artifact from not doing so, we did not bother.

Regarding our models, I doubt if we can share them. I will have to ask OptiTex -- if they will agree, I will put them on-line and post a message in this forum.