Corotated SPH for deformable solids
Posted: Fri Nov 18, 2011 3:11 am
Hi all,
I am currently trying to implement this paper titled Corotated SPH for deformable solids. Currently I am having problem understanding section 5.1. The authors extract the rotation matrix from the point neighborhood. Now in this calculation first the polar decomposition is discussed to extract R matrix similar to how meshless shapematching does. This requires calculating the COM(center of mass) of initial and deformed config. Pseudocode follows,
p_i = x_i - x_com
q_i =x0_i - x0_com
A = sum_i(m_i p_i q_i^T) sum_i(m_i q_i q_i^T) = A_(pq) A_(qq)
Since A_(qq) is symmetric it is ignored and A_(pq) is polar decomposed to get the symmetric part S=sqrt(A_(pq)^T A_(pq)) and rotational part R = A_(pq)S^(-1). All of this makes sense infact this is plane old meshless shape matching. This translates to this if i understood correctly,
Now the interesting bit comes in the second paragraph section 5.1 where they compute individual orientation Ri for each particle using the SPH formulation as follows,
A_(pq)_i = sum_j {m_j W(x0_(ij),h) ((x_j-x_i)(x0_j-x0_i)^T)}
This translates to the following code
Now the values of Apq are completely different by both methods. Would they be so?
Anyone who has any idea what may be the reason for this. Perhaps I am doing something wrong.
I am free for a discussion on this. Any ideas?
Regards,
Mobeen
I am currently trying to implement this paper titled Corotated SPH for deformable solids. Currently I am having problem understanding section 5.1. The authors extract the rotation matrix from the point neighborhood. Now in this calculation first the polar decomposition is discussed to extract R matrix similar to how meshless shapematching does. This requires calculating the COM(center of mass) of initial and deformed config. Pseudocode follows,
p_i = x_i - x_com
q_i =x0_i - x0_com
A = sum_i(m_i p_i q_i^T) sum_i(m_i q_i q_i^T) = A_(pq) A_(qq)
Since A_(qq) is symmetric it is ignored and A_(pq) is polar decomposed to get the symmetric part S=sqrt(A_(pq)^T A_(pq)) and rotational part R = A_(pq)S^(-1). All of this makes sense infact this is plane old meshless shape matching. This translates to this if i understood correctly,
Code: Select all
//Calculate center of mass (COM)
glm::vec3 Xcom, X0com;
for(int i=0;i<total_points;i++) {
X0com += Xi[i];
Xcom += X[i];
}
Xcom = Xcom/float(total_points);
X0com = X0com/float(total_points);
A = glm::mat3(0);
for(int i=0;i<total_points;i++) {
glm::vec3 pi = X[i]-Xcom;
glm::vec3 qi = Xi[i]-X0com;
A += M[i]* glm::outerProduct(pi, qi);
}
Apq[i] =A;
A_(pq)_i = sum_j {m_j W(x0_(ij),h) ((x_j-x_i)(x0_j-x0_i)^T)}
This translates to the following code
Code: Select all
glm::mat3 A = glm::mat3(0);
vector<neighbor>& pNeighbor = neighbors[index];
for(size_t i=0;i<pNeighbor.size();i++)
{
A += M[i]*pNeighbor[i].w * glm::outerProduct(pNeighbor[i].rdist, pNeighbor[i].rdist);
}
Apq[index] = A;
Anyone who has any idea what may be the reason for this. Perhaps I am doing something wrong.
I am free for a discussion on this. Any ideas?
Regards,
Mobeen