Regarding strain based dynamics [DONE]
-
- Posts: 122
- Joined: Thu May 05, 2011 11:47 am
Regarding strain based dynamics [DONE]
hi all.
I am trying to implement strain based dynamics and i have a slight problem with the implementation. I assume I am doing something wrong. This is for a triangular mesh so i am not considering the tetrahedral constraint.
There are two formulas given for calculating Strain for tri mesh. One in Eq. 11 which is
S = Q^-1 T P^T P Q^-1
Q is a 2x2 matrix, P is a 3x2 matrix so in the end S is a 2x2 matrix.
Now from explicit formulas given in Eq. 28, strain ij element can be calculated using
Sij = fi.fj = (Pci).(Pcj)
P is a 3x2 matrix
Q is a 2x2 matrix and c[i,j] is one column of it i.e. c[i,j] is 2x1 matrix. So Pci has dimensions 3x2*2x1 = 3x1 hence fi and fj are [3x1] vectors and fi.fj (dot product) will give us a scalar which we will place for ith row and jth col of S so S is a 2x2 matrix.
Now here is where the problem comes in. According to eq. 29 the correction of strain del Sij = fj ci^T+fi cj^T
fj and fi are 3x1 column vectors.
ci is a 2x1 vector and ci^T is 1x2 thus del Sij is going to give me 3x1 * 1x2 = 3x2 matrix
My questions are:
1) How do I add this del Sij to triangle vertex position? IIRC dp was a 3x1 column vector so it was directly added to positions of each vertex?
1) Why are there two formulas for calculating Strain? What will be the difference in output if i take the first or the second one?
I am trying to implement strain based dynamics and i have a slight problem with the implementation. I assume I am doing something wrong. This is for a triangular mesh so i am not considering the tetrahedral constraint.
There are two formulas given for calculating Strain for tri mesh. One in Eq. 11 which is
S = Q^-1 T P^T P Q^-1
Q is a 2x2 matrix, P is a 3x2 matrix so in the end S is a 2x2 matrix.
Now from explicit formulas given in Eq. 28, strain ij element can be calculated using
Sij = fi.fj = (Pci).(Pcj)
P is a 3x2 matrix
Q is a 2x2 matrix and c[i,j] is one column of it i.e. c[i,j] is 2x1 matrix. So Pci has dimensions 3x2*2x1 = 3x1 hence fi and fj are [3x1] vectors and fi.fj (dot product) will give us a scalar which we will place for ith row and jth col of S so S is a 2x2 matrix.
Now here is where the problem comes in. According to eq. 29 the correction of strain del Sij = fj ci^T+fi cj^T
fj and fi are 3x1 column vectors.
ci is a 2x1 vector and ci^T is 1x2 thus del Sij is going to give me 3x1 * 1x2 = 3x2 matrix
My questions are:
1) How do I add this del Sij to triangle vertex position? IIRC dp was a 3x1 column vector so it was directly added to positions of each vertex?
1) Why are there two formulas for calculating Strain? What will be the difference in output if i take the first or the second one?
Last edited by mobeen on Wed May 13, 2015 1:44 pm, edited 3 times in total.
-
- Posts: 122
- Joined: Thu May 05, 2011 11:47 am
Re: Regarding strain based dynamics
Nevermind, I think I know the answer. The del S that we calculate is a 3x2 matrix because it contains two column vectors of dp, one for point p1 and second for point p2. For p0, we just calculate it using the dp of p1 and p2. So in the end these are the corrections
calculate S using Eq. 28 or Eq. 11. Using either case, S is a 2x2 matrix
calculate del S using Eq. 29. del S is a 3x2 matrix
dp1 = getColumn1(del S);
dp2 = getColumn2(del S);
dp3 = -dp1-dp2 as given in eq. 31
So now I think it should work. Thanks anyway for reading my post
problem solved.
calculate S using Eq. 28 or Eq. 11. Using either case, S is a 2x2 matrix
calculate del S using Eq. 29. del S is a 3x2 matrix
dp1 = getColumn1(del S);
dp2 = getColumn2(del S);
dp3 = -dp1-dp2 as given in eq. 31
So now I think it should work. Thanks anyway for reading my post
![Smile :)](./images/smilies/icon_smile.gif)
-
- Posts: 122
- Joined: Thu May 05, 2011 11:47 am
Re: Regarding strain based dynamics
Dear all,
I am at it again. I am still finding difficulties in implementing strain based dynamics paper. I am using it to simulate a piece of cloth hence I need to worry about the triangle constraint. The bending modes work perfectly fine for my cloth but the triangle constraint fails. At the moment, i modify my existing pbd implementation by replacing distance constraints with triangle constraints. Here is how the StepPhysics function is implemented
The IntegrateExplicitWithDampingSBD function replaces IntegrateExplicitWithDamping function of PBD that uses center of mass based formulation. For SBD, it is implemented using Eq. 23 of paper as follows
The ComputeStrains function is defined as follows:
The matrix Q is precalculated at initialization using the uv coordinates as given in Appendix eq. 37, 38 as follows:
Finally, the UpdateTriangleConstraint is implemented as follows:
I do hope there is someone who can help. I have emailed the author of the paper too but I have not received any response from him.
Thanks,
Mobeen
I am at it again. I am still finding difficulties in implementing strain based dynamics paper. I am using it to simulate a piece of cloth hence I need to worry about the triangle constraint. The bending modes work perfectly fine for my cloth but the triangle constraint fails. At the moment, i modify my existing pbd implementation by replacing distance constraints with triangle constraints. Here is how the StepPhysics function is implemented
Code: Select all
void StepPhysics(float dt ) {
ComputeForces();
IntegrateExplicitWithDampingSBD(dt); //for strain based dynamics
ComputeStrains();
// for collision constraints
UpdateInternalConstraints();
UpdateExternalConstraints();
Integrate(dt);
}
Code: Select all
void IntegrateExplicitWithDampingSBD(float dt) {
size_t i=0;
float sumVN=0;
for(i=0;i<total_points;i++) {
V[i] *= global_dampening; //global velocity dampening !!!
V[i] = V[i] + (F[i]*dt)*W[i];
float lenDP = glm::length(dp[i]);
if(lenDP > EPSILON)
sumVN = sumVN + glm::dot(V[i], dp[i]/lenDP);
}
for(i=0;i<total_points;i++) {
float lenDP = glm::length(dp[i]);
if(lenDP > EPSILON)
V[i] = V[i] - kDamp*sumVN* dp[i]/lenDP;
}
//calculate predicted position
for(i=0;i<total_points;i++) {
if(W[i] <= 0.0) {
tmp_X[i] = X[i]; //fixed points
} else {
tmp_X[i] = X[i] + (V[i]*dt);
}
}
}
Code: Select all
void ComputeStrains() {
memset(&(S[0][0]),0,S.size()*sizeof(glm::mat2));
int count = 0;
//calculate material coordinates using Eq. 37 page 9
for(size_t i=0;i<indices.size();i+=3) {
int i0 = indices[i];
int i1 = indices[i+1];
int i2 = indices[i+2];
glm::vec3 x0 = X[i0];
glm::vec3 x1 = X[i1];
glm::vec3 x2 = X[i2];
glm::mat2x3 P = glm::mat2x3((x1-x0),(x2-x0));
//calculate strains using Eq. 11 page 4
glm::mat2x2 Qinv = Q[count];
glm::mat2x2 QinvT = glm::transpose(Qinv);
glm::mat3x2 P_T = glm::transpose(P);
glm::mat2x2 St = QinvT*P_T*P*Qinv;
S[count++] = St;
}
}
Code: Select all
int count = 0;
//calculate material coordinates using Eq. 37 page 9
for(i=0;i<indices.size();i+=3) {
int i0 = indices[i];
int i1 = indices[i+1];
int i2 = indices[i+2];
glm::vec3 x0 = X[i0];
glm::vec3 x1 = X[i1];
glm::vec3 x2 = X[i2];
glm::vec2 u0 = UV[i0];
glm::vec2 u1 = UV[i1];
glm::vec2 u2 = UV[i2];
glm::mat2x3 P = glm::mat2x3((x1-x0),(x2-x0));
glm::mat2x2 U = glm::mat2x2((u1-u0),(u2-u0));
glm::mat2x2 Uinv = glm::inverse(U);
glm::mat2x3 T = P*Uinv;
glm::vec3 n1(T[0]);//get first column of T
glm::vec3 n2(T[1]);//get second column of T
n1 = glm::normalize(n1);
n2 = glm::normalize(n2);
glm::mat2x2 C = glm::transpose(glm::mat2x3(n1,n2))*P;
Q[count++] = C;
}
Code: Select all
void UpdateTriangleConstraint(int index) {
TriangleConstraint c = t_constraints[index];
float w1 = W[c.p1];
float w2 = W[c.p2];
float w3 = W[c.p3];
float invMass = c.w;
if(invMass <= EPSILON)
return;
glm::mat2x2 St = S[index]; //contains sqrt(Strain)
glm::mat2x2 C = Q[index]; //contains Qinv
glm::vec3 x0 = tmp_X[c.p1];
glm::vec3 x1 = tmp_X[c.p2];
glm::vec3 x2 = tmp_X[c.p3];
glm::vec2 Ci = C[0];//get first column of C
glm::vec2 Cj = C[1];//get first column of C
glm::mat2x3 P = glm::mat2x3((x1-x0),(x2-x0));
glm::mat2x3 F = P*C;
// glm::vec3 PCi = P*Ci; //fi
// glm::vec3 PCj = P*Cj; //fj
glm::vec3 fi = F[0]; //get first column of F
glm::vec3 fj = F[1]; //get second column of F
// std::cout << "PCi.PCj [" << glm::dot(PCi, PCj) << "], fi.fj [" << glm::dot(fi, fj) << std::endl;
glm::mat2x3 fjci = outerProduct(fj,Ci);
glm::mat2x3 ficj = outerProduct(fi,Cj);
glm::mat2x3 dS = fjci + ficj;
glm::vec3 dp2 = dS[0];//(dS[0][0],dS[0][1],dS[0][2]);
glm::vec3 dp3 = dS[1];//(dS[1][0],dS[1][1],dS[1][2]);
glm::vec3 dp1 = -dp2-dp3;
//using eq. 35 and 36
/*
float lamda1 = 2.0f * ((St[0][0]-1.0f)/(invMass*glm::dot(dp1, dp1)))*(St[0][0]);
float lamda2 = (St[0][1] )/(invMass*glm::dot(dp2, dp2));
float lamda3 = 2.0f * ((St[1][1]-1.0f)/(invMass*glm::dot(dp3, dp3)))*(St[1][1]);
*/
//using eq. 33 and 34
float lamda1 = (St[0][0]-1.0f) / ( invMass*glm::dot(dp1, dp1) );
float lamda3 = (St[1][1]-1.0f) / ( invMass*glm::dot(dp3, dp3) );
float lamda2 = (St[0][1] ) / ( invMass*glm::dot(dp2, dp2) );
//std::cout << "l1 " << lamda1 << " l2: " << lamda2 << " l3: " << lamda3 << std::endl;
//std::cout << "Strain: [0][0] -> " << St[0][0] << " [1][1] -> " << St[1][1] << std::endl;
if(w1 > EPSILON)
dp[c.p1] += -lamda1*w1*dp1*c.k_prime;
if(w2 > EPSILON)
dp[c.p2] += -lamda2*w2*dp2*c.k_prime;
if(w3 > EPSILON)
dp[c.p3] += -lamda3*w3*dp3*c.k_prime;
}
Thanks,
Mobeen
-
- Posts: 4
- Joined: Fri Dec 30, 2011 1:51 pm
Re: Regarding strain based dynamics
Hi Mobeen,
Is this part correct?
"
//using eq. 33 and 34
float lamda1 = (St[0][0]-1.0f) / ( invMass*glm::dot(dp1, dp1) );
float lamda3 = (St[1][1]-1.0f) / ( invMass*glm::dot(dp3, dp3) );
float lamda2 = (St[0][1] ) / ( invMass*glm::dot(dp2, dp2) );
"
The denominators should be weighted sum of all the related terms, instead of a dot product, aren't they?
All the best,
Jianping
Is this part correct?
"
//using eq. 33 and 34
float lamda1 = (St[0][0]-1.0f) / ( invMass*glm::dot(dp1, dp1) );
float lamda3 = (St[1][1]-1.0f) / ( invMass*glm::dot(dp3, dp3) );
float lamda2 = (St[0][1] ) / ( invMass*glm::dot(dp2, dp2) );
"
The denominators should be weighted sum of all the related terms, instead of a dot product, aren't they?
All the best,
Jianping
-
- Posts: 122
- Joined: Thu May 05, 2011 11:47 am
Re: Regarding strain based dynamics
Hi,
The denominator is the weighted sum of squared length of delpkS_ii. I calculated the sum of weights (invSum) then did a dot product of delpkS_ii with itself since that should give me the squared length of delpkS_ii.
The denominator is the weighted sum of squared length of delpkS_ii. I calculated the sum of weights (invSum) then did a dot product of delpkS_ii with itself since that should give me the squared length of delpkS_ii.
-
- Posts: 4
- Joined: Fri Dec 30, 2011 1:51 pm
Re: Regarding strain based dynamics
Hi,mobeen wrote:Hi,
The denominator is the weighted sum of squared length of delpkS_ii. I calculated the sum of weights (invSum) then did a dot product of delpkS_ii with itself since that should give me the squared length of delpkS_ii.
I mean that, e.g., should lamda1 be as
float lamda1 = (St[0][0]-1.0f) / [ invMass1*(glm::dot(dp1, dp1)+invMass2*glm::dot(dp2, dp2)+invMass3*glm::dot(dp3, dp3) ) ;
?
-
- Posts: 30
- Joined: Thu Dec 19, 2013 12:13 pm
Re: Regarding strain based dynamics
Hi Mobeen,
are you planning to add your Strain Based Dynamics implementation to OpenCloth? That would be great.
I am looking at implementing the tetrahedral constraints presented in this paper for my medical simulator. Maybe you have got ideas for this as well?
Cheers
are you planning to add your Strain Based Dynamics implementation to OpenCloth? That would be great.
I am looking at implementing the tetrahedral constraints presented in this paper for my medical simulator. Maybe you have got ideas for this as well?
Cheers
-
- Posts: 122
- Joined: Thu May 05, 2011 11:47 am
Re: Regarding strain based dynamics
Hi korzen,Hi Mobeen,
are you planning to add your Strain Based Dynamics implementation to OpenCloth? That would be great.
I am looking at implementing the tetrahedral constraints presented in this paper for my medical simulator. Maybe you have got ideas for this as well?
Cheers
Yeah this will be integrated in the opencloth project. Moreover there have been some really good progress in the other integrators for instance, the IMEX integrator has been completely rewritten and it is performing really well. I will update the repository soon but before that i need to finish the strain based dynamics approach.
I have already emailed the SBD paper author regarding this error but he has not replied yet. Currently I have managed to finish the bending constraint but the triangle constraint is not done yet. Tetrahedral will be the next logical step after this is finished.
-
- Posts: 4
- Joined: Fri Dec 30, 2011 1:51 pm
Re: Regarding strain based dynamics
Hi Mobeen,mobeen wrote:Hi,
The denominator is the weighted sum of squared length of delpkS_ii. I calculated the sum of weights (invSum) then did a dot product of delpkS_ii with itself since that should give me the squared length of delpkS_ii.
Have you solved your problem?
Recently I implemented a simple 2D version (deformation in 2D plane only, so there is no 'material coordinates' part and bending constraints) and a 3D version as well.
And both are working well using the formulas in the paper.
Sincerely,
Jianping
-
- Posts: 122
- Joined: Thu May 05, 2011 11:47 am
Re: Regarding strain based dynamics
Hi JP,
Nope I have not been able to finish it yet. If possible could you please share your sources by emailing them on my NTU email address.
Thanks,
Mobeen
Nope I have not been able to finish it yet. If possible could you please share your sources by emailing them on my NTU email address.
Thanks,
Mobeen
-
- Posts: 122
- Joined: Thu May 05, 2011 11:47 am
Re: Regarding strain based dynamics [DONE]
OK thanks to help by JP, I was able to implement strain based dynamics both for triangular and tetrahedral objects.
Here is the snapshot of using strain based dynamics with tetrahedral mesh And here is the snapshot of using strain based dynamics with triangular mesh I will release the sources with the next OpenCloth update.
Here is the snapshot of using strain based dynamics with tetrahedral mesh And here is the snapshot of using strain based dynamics with triangular mesh I will release the sources with the next OpenCloth update.
You do not have the required permissions to view the files attached to this post.
-
- Posts: 30
- Joined: Thu Dec 19, 2013 12:13 pm
Re: Regarding strain based dynamics [DONE]
Well done Mobeen, looks really good!
Looking forward to grab the sources. I have partially ported the OpenCloth to Unity3D.
I will release the sources (+ Position Based Fluids on the GPU (DX11) when I am done, probably in the summer.
Bw
Korzen
Looking forward to grab the sources. I have partially ported the OpenCloth to Unity3D.
I will release the sources (+ Position Based Fluids on the GPU (DX11) when I am done, probably in the summer.
Bw
Korzen
-
- Posts: 122
- Joined: Thu May 05, 2011 11:47 am
Re: Regarding strain based dynamics [DONE]
That's nice korzen303. Keep me informed with ur developments.
-
- Posts: 122
- Joined: Thu May 05, 2011 11:47 am
Re: Regarding strain based dynamics [DONE]
Some more update. I have also managed to implement Barycentric mapping on top of strain based dynamics to apply the tetrahedral deformation to the surface mesh.
You do not have the required permissions to view the files attached to this post.
-
- Posts: 73
- Joined: Fri May 01, 2015 8:23 pm
Re: Regarding strain based dynamics [DONE]
Can you post a update or the files you used to strain based dynamics many of us would love to use it ![Smile :)](./images/smilies/icon_smile.gif)
![Smile :)](./images/smilies/icon_smile.gif)