Regarding strain based dynamics [DONE]

Please don't post Bullet support questions here, use the above forums instead.
mobeen
Posts: 122
Joined: Thu May 05, 2011 11:47 am

Regarding strain based dynamics [DONE]

Post by mobeen »

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?
Last edited by mobeen on Wed May 13, 2015 1:44 pm, edited 3 times in total.
mobeen
Posts: 122
Joined: Thu May 05, 2011 11:47 am

Re: Regarding strain based dynamics

Post by mobeen »

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.
mobeen
Posts: 122
Joined: Thu May 05, 2011 11:47 am

Re: Regarding strain based dynamics

Post by mobeen »

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

Code: Select all

void StepPhysics(float dt ) {
	ComputeForces();	 
	IntegrateExplicitWithDampingSBD(dt); //for strain based dynamics
	
	ComputeStrains();

	// for collision constraints	
	UpdateInternalConstraints();	
	UpdateExternalConstraints();
	
	Integrate(dt);
}
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

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);				 
		}
	} 
}
The ComputeStrains function is defined as follows:

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;
	}
}
The matrix Q is precalculated at initialization using the uv coordinates as given in Appendix eq. 37, 38 as follows:

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;  
}
Finally, the UpdateTriangleConstraint is implemented as follows:

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;
}
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
kxiaocai
Posts: 4
Joined: Fri Dec 30, 2011 1:51 pm

Re: Regarding strain based dynamics

Post by kxiaocai »

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
mobeen
Posts: 122
Joined: Thu May 05, 2011 11:47 am

Re: Regarding strain based dynamics

Post by mobeen »

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.
kxiaocai
Posts: 4
Joined: Fri Dec 30, 2011 1:51 pm

Re: Regarding strain based dynamics

Post by kxiaocai »

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.
Hi,

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) ) ;
?
korzen303
Posts: 30
Joined: Thu Dec 19, 2013 12:13 pm

Re: Regarding strain based dynamics

Post by korzen303 »

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
mobeen
Posts: 122
Joined: Thu May 05, 2011 11:47 am

Re: Regarding strain based dynamics

Post by mobeen »

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
Hi korzen,
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.
kxiaocai
Posts: 4
Joined: Fri Dec 30, 2011 1:51 pm

Re: Regarding strain based dynamics

Post by kxiaocai »

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.
Hi Mobeen,

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
mobeen
Posts: 122
Joined: Thu May 05, 2011 11:47 am

Re: Regarding strain based dynamics

Post by mobeen »

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
mobeen
Posts: 122
Joined: Thu May 05, 2011 11:47 am

Re: Regarding strain based dynamics [DONE]

Post by mobeen »

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
sbd_tetra.png
sbd_tetra.png (22.91 KiB) Viewed 24351 times
And here is the snapshot of using strain based dynamics with triangular mesh
sbd_cloth.png
sbd_cloth.png (29.44 KiB) Viewed 24351 times
I will release the sources with the next OpenCloth update.
korzen303
Posts: 30
Joined: Thu Dec 19, 2013 12:13 pm

Re: Regarding strain based dynamics [DONE]

Post by korzen303 »

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
mobeen
Posts: 122
Joined: Thu May 05, 2011 11:47 am

Re: Regarding strain based dynamics [DONE]

Post by mobeen »

That's nice korzen303. Keep me informed with ur developments.
mobeen
Posts: 122
Joined: Thu May 05, 2011 11:47 am

Re: Regarding strain based dynamics [DONE]

Post by mobeen »

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.
sbd_mapped.png
sbd_mapped.png (24.47 KiB) Viewed 24308 times
sbd_mapped_2.png
sbd_mapped_2.png (49.35 KiB) Viewed 24308 times
johnsonalpha
Posts: 73
Joined: Fri May 01, 2015 8:23 pm

Re: Regarding strain based dynamics [DONE]

Post by johnsonalpha »

Can you post a update or the files you used to strain based dynamics many of us would love to use it :)
Post Reply