My explicit euler cloth code

[mobeen]

Hi all,
In my quest to understand the cloth dynamics, I have developed a basic cloth (currently explicit integration has been implemented including Explicit Euler, Midpoint Euler and RK4 Methods). I want to share my code (visual studio 2009 sln attached) with the community to get the feedback from the experts. I am also looking for suggestions to improve the performance further. The code is as minimalistic as i possibly could make. For all the demo source codes I could find, all are using C++ objects which hides the details in the implementation. Hopefully, my code is going to make it easier for others (esp. beginners) to roll in their own cloth demos.
You will need freeglut, glew and glm libraries to build it. My objective is to share working code for cloth simulation.

I hope to get comments from all esp. Erwin.
Note that since explicit methods are unstable the code might explode but it depends on the params u set.
In the meantime, I am going through the implicit integration scheme and would like to do a minimalistic demo on that too.

Controls:
Left click and drag to rotate the view.
Left click and drag on any mass to select and then move it.

Snapshot:

Regards,
Mobeen

[mobeen]

50+ views, 4 downloads but no reply

[mobeen]

70+ view, 6 downloads no reply?
Could anyone give me some feedback please?

[DevO]

Hi!

First this is alway great if someone post source-code here, thanks you !

Second I think Explicit Euler and such integrations should not be used today.
Semi-implicit Euler that is used by Bullet is one preferable integration method.
Spring appears to be outdated too.

You could try to look at "Position Based Dynamics".
Here is good starting point for this.

https://github.com/roxlu/pbd

[mobeen]

Hi DevO,
Thanks for your response. Well I was reading these http://www.matthiasmueller.info/realtim ... enotes.pdf coursenotes from realtime physics course at SIGGRAPH 2010. As I was reading through the material, I thought of implementing it into my own applications and so I wrote the demo I shared.

Second I think Explicit Euler and such integrations should not be used today.

Could u tell me a reason why? I know it is unstable. Is there any other reason.
Infact the position based dynamics approach uses explicit Euler integration in the first step to calculate the new position and velocity. It is not assigned to the current position/velocity directly though.

Spring appears to be outdated too.

Who said so? The position based dynamics has given it a new fancy name called "Stretching Constraint" it is infact a spring if u look at the formulation.

The paper and technique u have given is also elaborated in these notes in great detail. Infact, the course was arranged by Matthias Mueller himself who introduced the position based dynamics. I will surely give it a try. Meanwhile, currently I am trying my hands on the implicit integration approach "Large steps in cloth simulation" by Baraff and Witkin but it is not that straight forward implementation. I hope to do it soon. So to sum up, this is how I am going to proceed. I want to create a single source file (main.cpp only) solution for all of these cloth simulation methods.

1) Explicit integration - done
2) Implicit integration - underway
3) Position based dynamics - added to queue
4) Semi implicit integration - added to queue

Do you think there is any newer approach that I must have a look into and which I have missed in the above list?

[DevO]

Well "position based dynamics" seems to use Semi-implicit Euler in the first step.

Semi-implicit Euler lock like this:
vel_new = vel_old + dt * force/mass
pos_new = pos_old + dt * vel_new

but Explicit Euler lock like this:
pos_new = pos_old + dt * vel_old
vel_new = vel_old + dt * force/mass

in PBD you just replaces pos_new with pos_predicted.

Springs modifying velocity based on position and can overshot so the are unstable.
Constraints in PBD are actually Position-Constraints and modifying predicted position based on predicted position.
Of course there are Velocity-Constraints but they work not like Springs.

You can also ready this "Soft Constraints" paper by Erin Catto.
http://code.google.com/p/box2d/downloads/list

Soft Constraints are like spring but should be stable using Semi-implicit Euler.

Of course Implicit integration is stable but you need to use solver linear like CG and this could be slow and of course harder to implement.

[mobeen]

Thanks for the clarification DevO.

[mobeen]

HI DevO,
I have started implementing the position based dynamics approach. I have setup the stretching constraint fine however, I am not sure how to handle the bending constraint. The demo code that u linked is doing the constraints using a different approach. I tried to do the stretching and bending constraints as given in it but it failed. So i went to look into the original paper and this is the stretching constraint I came up with after reading through the paper.

Code: Select all

void ApplyStretchConstaint(int i) {
   DistanceConstraint c = d_constraints[i];
   glm::vec3 dir = X[c.p1]-X[c.p2];
   float len = glm::length(dir);
   float invMass = w_i+w_i;
   glm::vec3 dP = (w_i/invMass) * (len-c.rest_length ) * (dir/len) * c.k;
   if(c.p1!=0 && c.p1 != numX)
      tmp_X[c.p1] -= dP;
   if(c.p2 != 0 && c.p2 != numX )
      tmp_X[c.p2] += dP;
}

This is working fine.

Now for the bend constraint the paper calculates the current dihedral angle and subtracts from the initial dihedral angle. I can obtain this angle by simple acos the dot product btw the two adjacent triangle normals and then subtract from the initial dihedral angle to get the change in the two dihedral angles. Now what should I do with this value? Should I rotate the two triangles by this amout?
From my thoughts, I think since hte bend constraint tries to approximate the stretch constraint across the two tris, I would need to rotate the tris across the common edge of the two tris. I may be wrong though. Any ideas? [EDIT: I did not read the full paper the appendix contains the definition of the bending and stretching constraints I will do through it and this http://bulletphysics.org/Bullet/phpBB3/ ... 3&start=15 thread also forced me to read the paper again]

Just updating on the status of my different cloth simulation demos.

(Green means finished, Blue means in progress and Red means finished but with some issues)
1) Explicit Integration

• a. Explicit Euler (this thread )
  b. Mid point Euler
  c. RK 4

2) Implicit Integration

• a. Implicit Euler (http://www.bulletphysics.org/Bullet/php ... f=4&t=7025)
  b. Baraff and Witkin’s Approach (http://www.bulletphysics.org/Bullet/php ... f=4&t=7051)

3) Semi Implicit Integration

• a. Semplectic Euler

4) Verlet Integration
5) Position based Dynamics

[DevO]

The code that I linked uses method from this paper.
A Triangle Bending Constraint Model for Position-Based Dynamics

This one is faster and easier to implement but has one problem, it does not work if two triangles are folded together so the angle is zero.

To project the bending constrain for the PBD paper one could try to use ideas from this paper.
Simulation of clothing with folds and wrinkles

An Accurate Model for Bending describes how to calculate the projection direction for all 4 (u1,u2,u3,u4) points of the bending constraint.

I see you have created opencloth project but why are you create one executable for all the simulation methods?
Would it be not better for comparison if only one executable witch all the method could be used?
So one could change the method dynamical (at runtime) and see the difference.

[Dirk Gregorius]

For bending constraints you just add another distance constraint over the shared edge of two triangles.

[mobeen]

Hi DevO and Dirk,
Thanks for the information.

I see you have created opencloth project but why are you create one executable for all the simulation methods?

Well currently, I am just trying to implement all of the algorithms individually to make sure they work fine alone. Once this is done, the next step will be to make a simple application combining all of the techniques so you may switch among them and see the results.

The code that I linked uses method from this paper.
A Triangle Bending Constraint Model for Position-Based Dynamics

Yeah, this I have already implemented nice and simple. Its the bending constraint in the original paper that is not working. I have noted the other paper that you have linked and will surely read that in a day or two.

For bending constraints you just add another distance constraint over the shared edge of two triangles.

Yeah this is what I thought too and i think this will work fine. However, I wanted to implement the technique exactly as given in the position based dynamics paper.

Thanks,
Mobeen

[mobeen]

Ok,
I am getting some results however after a while the cloth mesh disappears. I tried to debug it and i find that after some iterations, the values in q1,q2,q3 and q4 (eq. 25-28) in the original paper become 0. This is how I have setup the bend constraint may be someone here might know whats wrong.

Code: Select all

void UpdateBendingConstraint(int index) {
   size_t i=0;
   BendingConstraint c = b_constraints[index]; 
   float d = 0, phi=0,i_d=0;
   glm::vec3 n1=glm::vec3(0), n2=glm::vec3(0);
   
   glm::vec3 p1 = tmp_X[c.p1];
   glm::vec3 p2 = tmp_X[c.p2]-p1;
   glm::vec3 p3 = tmp_X[c.p3]-p1;
   glm::vec3 p4 = tmp_X[c.p4]-p1;

   glm::vec3 p2p3 = glm::cross(p2,p3);		n1 = glm::normalize(p2p3);
   glm::vec3 p2p4 = glm::cross(p2,p4);		n2 = glm::normalize(p2p4); 
  
   d	= glm::dot(n1,n2);
   phi = acos(d) ;
   i_d = sqrt(1-d*d)*(phi-phi0[index]) ;

   glm::vec3 p2n1 = glm::cross(p2,n1);
   glm::vec3 p2n2 = glm::cross(p2,n2);
   glm::vec3 p3n2 = glm::cross(p3,n2);
   glm::vec3 p4n1 = glm::cross(p4,n1);
   glm::vec3 n1p2 = -p2n1;
   glm::vec3 n2p2 = -p2n2;
   glm::vec3 n1p3 = glm::cross(n1,p3);
   glm::vec3 n2p4 = glm::cross(n2,p4);
   
   float lenp2p3 = glm::length(p2p3);
   float lenp2p4 = glm::length(p2p4);
  
   glm::vec3 q3 =   (p2n2 + n1p2*d)/ lenp2p3;
   glm::vec3 q4 =   (p2n1 + n2p2*d)/ lenp2p4;
   glm::vec3 q2 =  (-(p3n2 + n1p3*d)/ lenp2p3) - ((p4n1 + n2p4*d)/lenp2p4);
   glm::vec3 q1 = -q2-q3-q4;

   float q1_len = glm::length(q1);
   float q2_len = glm::length(q2);
   float q3_len = glm::length(q3);
   float q4_len = glm::length(q4); 
  
   float sum = M[c.p1]*(q1_len*q1_len) + 
                   M[c.p2]*(q2_len*q2_len) + 
                   M[c.p3]*(q3_len*q3_len) + 
                   M[c.p4]*(q4_len*q4_len);	
				 
   glm::vec3 dP1 = -( (M[c.p1] * i_d) /sum)*q1;
   glm::vec3 dP2 = -( (M[c.p2] * i_d) /sum)*q2;
   glm::vec3 dP3 = -( (M[c.p3] * i_d) /sum)*q3;
   glm::vec3 dP4 = -( (M[c.p4] * i_d) /sum)*q4;
  
   if(c.p1!=0 && c.p1 != numX) {
      tmp_X[c.p1] += dP1*c.k;
   }
   if(c.p2!=0 && c.p2 != numX) {
      tmp_X[c.p2] += dP2*c.k;
   }
   if(c.p3!=0 && c.p3 != numX) {
      tmp_X[c.p3] += dP3*c.k;
   }	
   if(c.p4!=0 && c.p4 != numX) {
      tmp_X[c.p4] += dP4*c.k;
   } 
}

[DevO]

You should test for degenerate and special cases, may be this is why you cloth mesh disappears.

Why do you use this ?

Code: Select all

float q1_len = glm::length(q1); //Sqrt
(q1_len*q1_len)

and not this

Code: Select all

float q1_slen = glm::length2(q1); //no Sqrt

Code: Select all

void UpdateBendingConstraint(int index) {
   size_t i=0;
   BendingConstraint c = b_constraints[index];
   float d = 0, phi=0,i_d=0;
   glm::vec3 n1=glm::vec3(0), n2=glm::vec3(0);
   
   glm::vec3 p1 = tmp_X[c.p1];
   glm::vec3 p2 = tmp_X[c.p2]-p1;
   glm::vec3 p3 = tmp_X[c.p3]-p1;
   glm::vec3 p4 = tmp_X[c.p4]-p1;

   glm::vec3 p2p3 = glm::cross(p2,p3);
   glm::vec3 p2p4 = glm::cross(p2,p4);

   float lenp2p3 = glm::length(p2p3); 
   if(lenp2p3 == 0.0) { return; } //need to handle this case.
   float lenp2p4 = glm::length(p2p4);
   if(lenp2p4 == 0.0) { return; } //need to handle this case.

   glm::vec3 n1 = p2p3 / lenp2p3; //normalize, lenp2p3 != 0.0
   glm::vec3 n2 = p2p4 / lenp2p4; //normalize, lenp2p4 != 0.0

   d   = glm::dot(n1,n2);
   if(d == -1.0){ // 180° special case, because (1.0-d*d) == 0.0 !
      return; //nothing to do 
   }

   phi = acos(d) ;
   if((phi-phi0[index]) == 0.0) return; //nothing to do 
   i_d = sqrt(1-d*d)*(phi-phi0[index]) ;

   glm::vec3 p2n1 = glm::cross(p2,n1);
   glm::vec3 p2n2 = glm::cross(p2,n2);
   glm::vec3 p3n2 = glm::cross(p3,n2);
   glm::vec3 p4n1 = glm::cross(p4,n1);
   glm::vec3 n1p2 = -p2n1;
   glm::vec3 n2p2 = -p2n2;
   glm::vec3 n1p3 = glm::cross(n1,p3);
   glm::vec3 n2p4 = glm::cross(n2,p4);
 
   glm::vec3 q3 =   (p2n2 + n1p2*d)/ lenp2p3;
   glm::vec3 q4 =   (p2n1 + n2p2*d)/ lenp2p4;
   glm::vec3 q2 =  (-(p3n2 + n1p3*d)/ lenp2p3) - ((p4n1 + n2p4*d)/lenp2p4);
   glm::vec3 q1 = -q2-q3-q4;

   float q1_len2 = glm::length2(q1);
   float q2_len2 = glm::length2(q2);
   float q3_len2 = glm::length2(q3);
   float q4_len2 = glm::length2(q4);
 
   float sum = M[c.p1]*(q1_len2) + M[c.p2]*(q2_len2) +
   M[c.p3]*(q3_len2) + M[c.p4]*(q4_len2);   
            
   glm::vec3 dP1 = -( (M[c.p1] * i_d) /sum)*q1;
   glm::vec3 dP2 = -( (M[c.p2] * i_d) /sum)*q2;
   glm::vec3 dP3 = -( (M[c.p3] * i_d) /sum)*q3;
   glm::vec3 dP4 = -( (M[c.p4] * i_d) /sum)*q4;
 
   if(c.p1!=0 && c.p1 != numX) {
      tmp_X[c.p1] += dP1*c.k;
   }
   if(c.p2!=0 && c.p2 != numX) {
      tmp_X[c.p2] += dP2*c.k;
   }
   if(c.p3!=0 && c.p3 != numX) {
      tmp_X[c.p3] += dP3*c.k;
   }   
   if(c.p4!=0 && c.p4 != numX) {
      tmp_X[c.p4] += dP4*c.k;
   }
}

[mobeen]

Hi DevO,
Thanks for the response.

Why do you use this ?

Code: Select all

code removed

and not this

Code: Select all

Code removed

Well i think I was following the paper a bit too much I guess. Now how could i miss that? Thanks for spotting this big mistake. Infact, the squared length could even be replaced by a dot product if i am correct

Code: Select all

float q1_len2 = glm::length(q1)*glm::length(q1);

is equivalent to

Code: Select all

float q1_len2 = glm::dot(q1,q1);

For the special case, d==-1, this means that the triangles are facing in the opposite direction? I think we would need to adjust the positions so that the dot product is positive again? I am not sure whether this is correct but any ideas?

Apart from this, do u think the rest of the code is fine?

[DevO]

For the special case, d==-1, this means that the triangles are facing in the opposite direction?

Yes this is correct.
Probably the code should better look like this.

Code: Select all

   d   = glm::dot(n1,n2);
   //try to catch invalid values that will return NaN.
   // sqrt(1 - (1.0001*1.0001)) = NaN 
   // sqrt(1 - (-1.0001*-1.0001)) = NaN 
   if(d<-1.0) d = -1.0; else if(d>1.0) d=1.0; //d = clamp(d,-1.0,1.0);

   //in both case sqrt(1-d*d) will be zero and nothing will be done.
   if(d == -1.0){ //0° case, the triangles are facing in the opposite direction, folded together.
      //TODO:
   }
   if(d == 1.0){ //180° case, the triangles are planar
      //TODO:
   }

Code: Select all

glm::length2(q4);

should be the same as

Code: Select all

glm::dot(q4,q4);

but first one is easier to read.

Well all other parts are looking ok to me, but I am not sure without have tested it.