Soft bodies early preview

[Nathanael]

Here's a (very) early preview of soft bodies implementation for Bullet.

http://bulletphysics.com/ftp/pub/test/i ... ftBody.zip

Basic keys:
- ']' Next sub demo.
- [space] Reset scene.

All comments/critics are welcome, just don't ask THE question, "when?" ^^.

Enjoy.

[Dirk Gregorius]

Nice demos! Very well done

I assume this is an implementation of a position based solver (e.g. like Jacobsen or Muller)? How do you handle the coupling between the cloth and the rigid bodies? Do you solve both in one loop while correcting position for the cloth and velocities for the rigid body e.g. using the displacement of the particle divided by the timestep. Or is the cloth and rigid body handled on the position level?


Cheers,
-Dirk

[Nathanael]

Nice demos! Very well done

Thanks,

Yes, the solver itself is based on Matthias Muller's work, mainly because its simple to implement and provide velocities which are needed for some future features.

Coupling between rigid and soft bodies is indeed currently done inside the (soft bodies) constraints solver loop, through impulses conversion from position displacements, it work well for a demo, but the constraints solver generate too much 'noise/jumps' so i had to damp those impulses with a arbitrary factor, its not clean, and make things harder to use, so I'm currently investigating alternative approaches will not introducing coupling between soft and rigid solvers.

[DevO]

Great Demos!!!

Are you working alone or with Erwin on this?

As I can see the collision has as usual problems with tunneling.
What are you using for collision between soft-bodies and rigid-bodies?

[Dirk Gregorius]

So you solve the rigid bodies and the particles in two different loops and only in the particle loop you apply impulse derived from the particle displacements? Did you try to solve both in the same loop, e.g.

Code: Select all

for ( n iterations )
  {
  for ( all rigid bod constraints c )
    c->Satisfy();

  for ( all particle constraints c )
    c->Satisfy();

  for ( all coupling constraints c )
     c->Satisfy();
  }

The final sweep would adjust positions for the particles and apply impulses on the rigid bodies derived from the particle correction.


Since you mention the Muller paper. For the soft bodies you simply use tetrahedral meshes with volume preserving constraints and distance constraints between the nodes? Basically:

C(p1, p2) = |p2 p1|- L0 = 0
C(p1, p2, p3, p4) = |(p2 - p1) * (p3 - p1) x (p4 - p1)| / 6 - V0 = 0

[Nathanael]

Great Demos!!!

Are you working alone or with Erwin on this?

As I can see the collision has as usual problems with tunneling.
What are you using for collision between soft-bodies and rigid-bodies?

Thanks a lot,
Yes, tunneling is quite obvious at this point, and the effect is worsen by the fact that the soft solver is one time step late (see Dirk comment).

I am working with Erwin for the integration into Bullet.

For collision i currently use a dynamic signed distance field built on the fly during simulation via gjk/epa.

[Nathanael]

So you solve the rigid bodies and the particles in two different loops and only in the particle loop you apply impulse derived from the particle displacements? Did you try to solve both in the same loop, e.g.
The final sweep would adjust positions for the particles and apply impulses on the rigid bodies derived from the particle correction.

Since you mention the Muller paper. For the soft bodies you simply use tetrahedral meshes with volume preserving constraints and distance constraints between the nodes? Basically:

C(p1, p2) = |p2 p1|- L0 = 0
C(p1, p2, p3, p4) = |(p2 - p1) * (p3 - p1) x (p4 - p1)| / 6 - V0 = 0

Yes i rigid and soft bodies are solved in two different loops, i didn't try to solve them in the same loop because I'm sure it will work much better but that create a coupling between the two that I'd like to dodge. just few of my rationals:
- Iterations number for rigid simulation should, in an ideal world be set to +inf, because we want to solve the system, but in soft bodies, it's more a parameter of the system, that can vary from bodies to bodies for 'aesthetic' reasons.
- Almost by definition, soft bodies create 'soft' constraints on rigid bodies, it should be possible to get nice results using only impulses outside the loop (at least i try to).
- Code/architecture, i don't know or maintain the rigid simulation of Bullet, and peoples who do surely don't want to have include the soft bodies part each time a change occur or new ideas are tested for rigid bodies. That's a trade off between accuracy and simplicity.

For volume, three different methods are demonstrated,
- the demo named 'Pressure' use forces, as describe in 'Pressure soft bodies model' from Maciej Matyka.
- the demo named 'Volume' uses forces proportional to (V-V0)*area per vertex in the direction of the normal, there is no tetrahedrons involved, just global volume.
- the demos named 'TorusMatch' and 'BunnyMatch' use polar decomposition to extract linear transforms (just rotation for the moment) from the deformed model, and 'pulling' vertices's toward the reference shape, see 'shape matching' from Matthias Muller (him again ). Still no tetrahedron.

The reason i don't make use tetrahedrons is twofold, first, authoring, i don't know of any commonly available 3d package that output tetrahedral meshes, second, the three previously mentioned methods work well for volume conservation, and i have ideas to 'emulate' finite elements (for plasticity and fracture, the next step after volume conservation) using the properties of polar decomposition.

Nat.

[Dirk Gregorius]

Thanks for the links regarding the soft bodies. I need to look them up so we can maybe discuss them a little more. Still I would be interested in the volume constraints since it will be fast and simple (if it works). There is a simple open source tool that can create the tetra mesh for you. I need to look it up tomorrow in the company and post it here...

Again - great work!
-Dirk

[Nathanael]

Still I would be interested in the volume constraints since it will be fast and simple (if it works). There is a simple open source tool that can create the tetra mesh for you.

I just finished deriving the tetrahedral volume constraint for positions, as you said, its simple enough, so I'll give it a try. For speed, we need to get rid of the cubic root, but I guess the approximation (power series expansion) that Thomas Jakobsen used for distance constraints should work there too.

It may take a while before i get back to you a that, because right now the priority is Bullet integration. Concerning that tool to generate tetra's, does it work with triangles soup (or a least slightly non manifold/open) meshes?

Thanks again, Nat.

[Dirk Gregorius]

I am sure you did this, but you can simplify the derivation of the volume constraint by making p1 the origin. Basically you have p1, p2, p3 and p4. The volume constraint is:

C = |(p2 - p1) * (p3 - p1) x (p4 - p1)| / 6 - V0

Now we move p1 into the origin:

p1' = p1 - p1 = ( 0,0,0 )
p2' = p2 - p1
p3' = p3 - p1
p4' = p4 - p1

This way the constraint becomes:

C' = |p2' * p3' x p4'| / 6 - V0

The gradients of C and C' are the same, but it is a little bit easier to find the gradients of C' in my opinion. The following holds then:

Define V = p2' * p3' x p4'

dC'/dp4 = dC/dp4 = sign( V ) * (p2' x p3') / 6
dC'/dp3 = dC/dp3 = sign( V ) * (p4' x p2') / 6
dC'/dp2 = dC/dp2 = sign( V ) * (p3' x p4') / 6

dC/dp1 = -( dC/dp4 + dC/dp3 + dC/dp2) = -( dC'/dp4 + dC'/dp3 + dC'/dp2)


I also derived the area constraint for a triangle A = |(p2 - p1) x (p3 - p1)| / 2. Basically I always apply the above technique for all these constraints. Also for the simple distance constraint.


The tetra generator can be found here. You need to look up its limitations since I don't know from the back of my had if it works for your requirements: http://tetgen.berlios.de/


I now look at the references you mention,
-Dirk

[Antonio Martini]

Yes, the solver itself is based on Matthias Muller's work, mainly because its simple to implement and provide velocities which are needed for some future features.

How is your solver related to Matthias Muller's work? is it using the same method?

cheers,
Antonio

[Dirk Gregorius]

I edited my post above and added the results for the gradients. How do you get a cubic root?

[Nathanael]

How is your solver related to Matthias Muller's work? is it using the same method?

Hi Antonio,
In fact, a tiny bit of it, but useful indeed, at the end of the step, Muller explicitly evaluate the velocity where Jacobsen implicitly evaluate it. That's not much, but i didn't thought about it before, that's what we call thinking too much inside the box i guess. , the rest of it is a classic Jacobsen iterative solver.
For bending constraints i use the same Jacobsen distance constraints as for structural constraints, generated from the adjacency graph of the soft body mesh topology.
Some of the demos use polar decomposition the extract a rigid frame from the deformed model, as explained in great details in 'Matrix Animation and Polar Decomposition, Ken Shoemake.', the decomposition implementation method itself come from 'The Polar Decomposition Properties, Applications And Algorithms, Pawel Zielinski, Krystyna Zietak'. Muller make use of it in one paper, but goes beyond linear deformation, and far beyond my understanding...

Nat.

[Nathanael]

I edited my post above and added the results for the gradients. How do you get a cubic root?

First, I'd like you to understand that I am probably better at coding than at solving equations, and probably better at cooking than at coding...
I did drop p1 too, more zeros, better i feel.
Then i am solving (p2'*k) * (p3'*k) x (p4'*k) = V0 for k, doing so, i am looking the scaling factor that bring my tetrahedron to it rest volume.
Let V be p2' * p3' x p4'
I got k = 1 / (V/V0)^1/3 . Applying that factor from the center of mass do the job.
The method work for distance, triangle, tetrahedron, etc... but can't figure why we got different results.

[Dirk Gregorius]

So you do this the following?

c = ( p1 + p2 + p3 + p4 ) / 4

p1 += k * (p1 - c) / |p1 - c|
p2 += k * (p2 - c) / |p2 - c|
p3 += k * (p3 - c) / |p3 - c|
p4 += k * (p4 - c) / |p4 - c|

This would be an interesting approach. I am not sure if this is physical correct though. The displacements should be in the direction of the gradients to be kind of "workless"". I need to think about this. The gradients I posted above are the derivatives of C w.r.t p_i where i = 1..4.


Or do you just do this:

c = ( p1 + p2 + p3 + p4 ) / 4

p1 -= c;
p2 -= c;
p3 -= c;
p4 -= c;

p1 *= k;
p2 *= k;
p3 *= k;
p4 *= k;

p1 += c;
p2 += c;
p3 += c;
p4 += c;



Cheers,
-Dirk