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,
Code:
//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;
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
Code:
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;
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
Login
Register
FAQ
Search
