My explicit euler cloth code
My explicit euler cloth code
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
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
 Attachments

 ClothCPU.zip
 (6.04 KiB) Downloaded 849 times
Re: My explicit euler cloth code
50+ views, 4 downloads but no reply
Re: My explicit euler cloth code
70+ view, 6 downloads no reply?
Could anyone give me some feedback please?
Could anyone give me some feedback please?
Re: My explicit euler cloth code
Hi!
First this is alway great if someone post sourcecode here, thanks you !
Second I think Explicit Euler and such integrations should not be used today.
Semiimplicit 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
First this is alway great if someone post sourcecode here, thanks you !
Second I think Explicit Euler and such integrations should not be used today.
Semiimplicit 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
Re: My explicit euler cloth code
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.
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.
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?
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.
Could u tell me a reason why? I know it is unstable. Is there any other reason.Second I think Explicit Euler and such integrations should not be used today.
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.
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.Spring appears to be outdated too.
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?
Re: My explicit euler cloth code
Well "position based dynamics" seems to use Semiimplicit Euler in the first step.
Semiimplicit 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 PositionConstraints and modifying predicted position based on predicted position.
Of course there are VelocityConstraints 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 Semiimplicit 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.
Semiimplicit 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 PositionConstraints and modifying predicted position based on predicted position.
Of course there are VelocityConstraints 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 Semiimplicit 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.
Re: My explicit euler cloth code
Thanks for the clarification DevO.
Re: My explicit euler cloth code
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.
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
4) Verlet Integration
5) Position based Dynamics
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) * (lenc.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;
}
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
 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)
 a. Semplectic Euler
4) Verlet Integration
5) Position based Dynamics
Re: My explicit euler cloth code
The code that I linked uses method from this paper.
A Triangle Bending Constraint Model for PositionBased 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.
A Triangle Bending Constraint Model for PositionBased 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.

 Posts: 874
 Joined: Sun Jul 03, 2005 4:06 pm
 Location: Kirkland, WA
Re: My explicit euler cloth code
For bending constraints you just add another distance constraint over the shared edge of two triangles.
Re: My explicit euler cloth code
Hi DevO and Dirk,
Thanks for the information.
Thanks,
Mobeen
Thanks for the information.
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.DevO wrote: I see you have created opencloth project but why are you create one executable for all the simulation methods?
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.DevO wrote: The code that I linked uses method from this paper.
A Triangle Bending Constraint Model for PositionBased Dynamics
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.Dirk wrote: For bending constraints you just add another distance constraint over the shared edge of two triangles.
Thanks,
Mobeen
Re: My explicit euler cloth code
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. 2528) in the original paper become 0. This is how I have setup the bend constraint may be someone here might know whats wrong.
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. 2528) 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(1d*d)*(phiphi0[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 = q2q3q4;
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;
}
}
Re: My explicit euler cloth code
You should test for degenerate and special cases, may be this is why you cloth mesh disappears.
Why do you use this ?
and not this
Why do you use this ?
Code: Select all
float q1_len = glm::length(q1); //Sqrt
(q1_len*q1_len)
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.0d*d) == 0.0 !
return; //nothing to do
}
phi = acos(d) ;
if((phiphi0[index]) == 0.0) return; //nothing to do
i_d = sqrt(1d*d)*(phiphi0[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 = q2q3q4;
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;
}
}
Re: My explicit euler cloth code
Hi DevO,
Thanks for the response. is equivalent to
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?
Thanks for the response.
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 correctDevO wrote: Why do you use this ?and not thisCode: Select all
code removed
Code: Select all
Code removed
Code: Select all
float q1_len2 = glm::length(q1)*glm::length(q1);
Code: Select all
float q1_len2 = glm::dot(q1,q1);
Apart from this, do u think the rest of the code is fine?
Re: My explicit euler cloth code
Yes this is correct.For the special case, d==1, this means that the triangles are facing in the opposite direction?
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(1d*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);
Code: Select all
glm::dot(q4,q4);
Well all other parts are looking ok to me, but I am not sure without have tested it.