My explicit euler cloth code

mobeen · **Joined:** Thu May 05, 2011 11:47 am **Posts:** 61

HI DevO,
Thanks for the reply once again. Are these special cases discussed in the pos based dynamics paper? Atleast, I could not find it. I will read through the other paper to see if it has something to offer. What do u suggest how to handle these cases?

DevO · **Joined:** Fri Mar 31, 2006 7:13 pm **Posts:** 90

mobeen wrote:

HI DevO,
Thanks for the reply once again. Are these special cases discussed in the pos based dynamics paper? Atleast, I could not find it. I will read through the other paper to see if it has something to offer. What do u suggest how to handle these cases?

No the special cases are usually ignored in the parers.
First test if it work if you just ignore (return;) this cases.
The case d == 1.0 is the common case for you simulation because you start with the plane and all angles are 180° at the beginning.
So probably you phi0 is 0.0 (180°) too, in this case this code will be correct.

Code:

   if(d == 1.0){ //180° case, the triangles are planar
        phi = 0.0;  //acos(1.0) == 0.0
        if(phi == phi0[index]) return; //nothing to do 
   }

Code:

   if(d == -1.0){ //0° case, the triangles are facing in the opposite direction, folded together.
       phi = PI;  //acos(-1.0) == PI
       if(phi == phi0[index]) return; //nothing to do 
      //in this case one just need to push vertices 1 and 2 in n1 and n2 directions, so the constrain will do the work in second iterations.
   }

Dirk Gregorius · **Posted:** Thu Jul 21, 2011 3:14 pm

The bending constraint using the dihedral angle is not really useful for games. It is much more expensive without improving the simulation quality significantly. Even PhysX is either not using it anymore or at least recommending to use the method I suggested earlier. You can check their manual.

Here are some other links related to position based dynamics:
http://www.matthiasmueller.info/publications/hpbd.pdf
http://www.matthiasmueller.info/publica ... Meshes.pdf
http://www.matthiasmueller.info/publica ... ticles.pdf
http://web.archive.org/web/200706102238 ... dc2001.htm

DevO · **Joined:** Fri Mar 31, 2006 7:13 pm **Posts:** 90

The method presented in this paper is probably only a bit slower as a simple distance constraint.
A Triangle Bending Constraint Model for Position-Based Dynamics
Probably one could use the same optimisation to remove Sqrt() as in distance constraint too.
One problem is if the angle is 0° then it does not work.

For simple distance constraint 180° case will not work, because the vertices will be moved in the plane.

Dirk Gregorius · **Posted:** Thu Jul 21, 2011 4:35 pm

This is all very theoretical. I will be convinced if I see a *significant* improvement on a game character using a special bending constraint which justifies the extra cycles.

mobeen · **Joined:** Thu May 05, 2011 11:47 am **Posts:** 61

HI DevO and Dirk,
Thanks for the responses.

Since I am doing the code exactly as in the original paper, the triangle order is given in a form where the first triangle is CCW and the second is CW winded. This will turn the direction of normals to 180 degrees. So my initial phi are all PI. The first and second vertex are forming the shared edge. For d==-1 special case in the condition when I need to move vertex 1 and 2, are these the shared vertices that I move or the other non-shared vertices? My intuition says it has to be the shared vertices.

Doing these changes, now i see some vertices popping up and down. I am currently only using the bending constraints, I have slowed the simulation speed and the resulting mesh looks like a crumpled piece of paper. Is this how it should look with only the bending constraints? Note that the alternating normals are flipped.

Dirk I will go through those links thanks for letting me know.

EDIT: I just added in the distance constraints however, the cloth seems very stiff and rigid as if it is inside a viscous fluid. It converges pretty quickly though. On rapid movements, it explodes? Any ideas? I could attach the code if u would want to have a look?

DevO · **Joined:** Fri Mar 31, 2006 7:13 pm **Posts:** 90

Hi mobeen,

I think that you triangle order is the source of you problems.
This bending constraint will only work if all triangles have the same orientation, ether all CW or all CCW.
The phi0 should be 0.0 for the plane.

Quote:

For d==-1 special case in the condition when I need to move vertex 1 and 2, are these the shared vertices that I move or the other non-shared vertices?

Shared vertices are vertices that form the edge ? If so then you need to move another 2 vertices in opposite direction so the triangles will be unfolded again.

mobeen · **Joined:** Thu May 05, 2011 11:47 am **Posts:** 61

Hi DevO,
Changing both triangle order to same, makes the cloth to crumpled even more. I have double checked the phi0 values and they are all 0.

My current code makes cloth triangles from a quad. Lets say using v0 v1 v5 v6 now the common edge that is shared by the triangles is v0v6 and the non shared (apex vertices are v5 and v1). This is how they are arranged exactly as in the position based dynamics paper. (I m using ASCII fig here so bare with me)

Code:

   v0                 v1              v2... 
      +-----------+-----------+ ...
     /  \             /               /
    /       \        /               /
   /            \   /               /
 +-----------+-----------+ ...
v5              v6

I am attaching the code if u think it may help. See if you can spot the problem/s.

Thanks for your help DevO.

DevO · **Joined:** Fri Mar 31, 2006 7:13 pm **Posts:** 90

Hi mobeen,

I have mode some changes to you code, so it locks more like PBD.
Unfortunately it is still not not fully correct.
For example you "center of mass's damping" seems to be broken.

New array W is to store inverse particle mass if this is zero then particle is fixed and will not move but the constrains need to know this!

regards,
DevO

mobeen · **Joined:** Thu May 05, 2011 11:47 am **Posts:** 61

Hi DevO,
Thanks for the corrections DevO. Yeah the handling of mass weights like this should have been done from the start.

Quote:

Unfortunately it is still not not fully correct. For example you "center of mass's damping" seems to be broken.

Thanks for the fast feedback DevO u r a life saver

I will update you with another version soon and once u signal me that it is fine, I will commit it to SVN.

EDIT: I tried to isolate the update of distance constraints only but on setting the ks=1 value, the simulation explodes. I cannot figure out why this is so when I am pretty sure the equations are correct. I am trying to see if I could solve this issue. I think I will do the COM damping after this Ks=1 issue is solved.

mobeen · **Joined:** Thu May 05, 2011 11:47 am **Posts:** 61

OK DevO I think I know what was causing the explosion. IN the distance update function, I was checking the original positions each iteration when instead I should check the predicted positions (as mentioned in the pseudocode in the original paper). Specifically, I changed this

Code:

glm::vec3 dir = X[c.p1] - X[c.p2];

to this

Code:

glm::vec3 dir = tmp_X[c.p1] - tmp_X[c.p2];

in the UpdateDistanceConstraint func. and now the explosion is gone (atleast in the distance constraint case).

Now I am working on the center of mass damping.

EDIT: Ok I think I have got something working now.
I think the original paper has some typo in the center of mass damping calculation. The algorithm line 7 in section 3.5 page 5

Code:

delVi = Vcm + cross(w,Ri[i])-V[i]

should be

Code:

delVi = Vcm + cross(Ri[i],w)-V[i]

because the second term is probably calculating the angular torque which should be the cross product of vector ( ri=Xi-Xcm) to the angular velocity (omega). My reference for this is http://en.wikipedia.org/wiki/Torque . what's ur say DevO am I correct?

Now for Kdamp value 0.5, I get the cloth falling down but then as soon as it is about to bend, it undergoes excessive damping. I tried Kdamp value of 1 but then the simulation explodes. I am still working on this to make it stable.

DevO · **Joined:** Fri Mar 31, 2006 7:13 pm **Posts:** 90

Hi,

the parer equation should be correct!
I was using it and it was working for me.

This should work.

Quote:

delVi = Vcm + cross(w,Ri[i])-V[i]

This line appears to be wrong.

Quote:

I += glm::outerProduct(ri, ri)*mass;

Cross product Matrix

There could also be numerical problems because of accumulation.

Kdamp =1.0 actually should work, and the plane should move then like rigid body.

regards,
Devid

mobeen · **Joined:** Thu May 05, 2011 11:47 am **Posts:** 61

Hi DevO,
Thanks for the update. Well I tried to run the code but the cloth rotates and then explodes. Anyways,
why do u recalculate the mass each iteration when it is a constant? and why do u say this

Code:

const glm::vec3 delVi = Vcm + glm::cross(ri,w) - V[i];

is wrong? I checked the definition of torque it is the cross product of ri with w.

Quote:

This line appears to be wrong.

Code:

I += glm::outerProduct(ri, ri)*mass;

Is this correct?

Code:

glm::mat3 tmp = glm::mat3(0,-Ri[i].z,  Ri[i].y, 
              Ri[i].z,       0,-Ri[i].x,
           -Ri[i].y,Ri[i].x,      0);
I += (tmp*Ri[i])*mass;

Could u recheck your old pbd implementation to see where I am wrong?

mobeen · **Joined:** Thu May 05, 2011 11:47 am **Posts:** 61

Ok DevO,
Thanks I think I have got it. The actual calculation of this

Code:

 
glm::mat3 tmp = glm::mat3(0,-Ri[i].z,  Ri[i].y, 
              Ri[i].z,       0,-Ri[i].x,
           -Ri[i].y,Ri[i].x,      0);
I += (tmp*Ri[i])*mass;

should be this

Code:

I += (tmp*glm::transpose(tmp))*mass;

And the rest is all fine. I got this from http://www.sccg.sk/~onderik/phd/ca2010/ ... sson11.pdf

DevO could u check this final version for me. I have corrected the triangle based bending constraint as well. You can undefine USE_TRIANGLE_BENDING_CONSTRAINT if u want to use it.

One more thing, what is your real name and affiliation, I will add in your name and affiliation in the code.

Thanks for the help man. You rock.

Cheers,
Mobeen

DevO · **Joined:** Fri Mar 31, 2006 7:13 pm **Posts:** 90

Quote:

Anyways,
why do u recalculate the mass each iteration when it is a constant?

Because this calculation is really cheap, to store and the query need more time and memory as to recalculate it.

Yes now witch transpose() it should be correct.

In you code you still use constant mass but this is not correct because fixed point should actually have infinite mass, that can be replaced in real implementation by the big number. (1000?)

Quote:

Meanwhile, I am working on the other triangle based bending constraint u told me about.

This should be 3 point constrain only.
The only tricky part is actually properly creating this constrains.
The should be just like distance constrains that are used to replace bending constrains but are using middle point too.

Code:

0-----1-----2
|     |     |
3-----4-----5
|     |     |
6-----7-----8

In this simple case you should have 6 bend constrains.
C(0,1,2), C(3,4,5), C(6,7,8), C(0,3,6), C(1,4,7), C(2,5,8).

Physics Simulation Forum

My explicit euler cloth code

Who is online