Meshless FEM [SOLVED]
Posted: Sun Sep 25, 2011 6:02 am
Hi all,
I am trying to implement the following paper "Point Based Animation of Elastic, Plastic and Melting Objects" available:http://www.matthiasmueller.info/publications/sca04.pdf however, currently the output I am getting is wrong my points just simply fall down due to gravity and then explode. This is how I have done it so far if someone here has any idea of where I might be wrong please help.
From what I have understood, basically, there are four steps in the algorithm:
1) Estimation of the moment matrix for each point (sum (Xj-Xi)(Xj-Xi)^T Wij) and then inverting it.
2) Calculate external forces (gravity/collision/wind etc.)
3) Calculate the internal forces due to stresses between the current point and neighbour.
4) Integration of the velocities and then obtaining new positions from them.
For 1) this is how I am calculating it.
2) is obvious
3) is more involved and this is how i calculate the internal forces. D is a vec3 containing the three coefficients of the elasiticity matrix.
4) Currently I am using explicit Euler integration from an existing code so this is fine.
Thanks,
Mobeen
I am trying to implement the following paper "Point Based Animation of Elastic, Plastic and Melting Objects" available:http://www.matthiasmueller.info/publications/sca04.pdf however, currently the output I am getting is wrong my points just simply fall down due to gravity and then explode. This is how I have done it so far if someone here has any idea of where I might be wrong please help.
From what I have understood, basically, there are four steps in the algorithm:
1) Estimation of the moment matrix for each point (sum (Xj-Xi)(Xj-Xi)^T Wij) and then inverting it.
2) Calculate external forces (gravity/collision/wind etc.)
3) Calculate the internal forces due to stresses between the current point and neighbour.
4) Integration of the velocities and then obtaining new positions from them.
For 1) this is how I am calculating it.
Code: Select all
for( k=0;k<numZ;k++) {
for( j=0;j<numY;j++) {
for( i=0;i<numX;i++) {
int current_index =(((k*numY)+j)*numX)+i;
glm::mat3 sumM=glm::mat3(1);
for(int m=0;m<18;m++) {
glm::vec3 n = getNextNeighbour(m);
int ix = int(n.x+i);
int iy = int(n.y+j);
int iz = int(n.z+k);
int index = (((iz*numY)+iy)*numX)+ix;
if (index < 0)
index =0;
if (index >=total_points)
index = total_points-1;
glm::vec3 Xij = X[index]-X[current_index];
float r = glm::length(Xij);
float h = 1;
float k = kernel(r,h);
rho_i[current_index] += k*mass;
glm::mat3 tmp = glm::mat3(1);
tmp[0][0] = Xij.x*Xij.x; tmp[0][1] = Xij.x*Xij.y; tmp[0][2] = Xij.x*Xij.z;
tmp[1][0] = Xij.y*Xij.x; tmp[1][1] = Xij.y*Xij.y; tmp[1][2] = Xij.y*Xij.z;
tmp[2][0] = Xij.z*Xij.x; tmp[2][1] = Xij.z*Xij.y; tmp[2][2] = Xij.z*Xij.z;
sumM += k*tmp;
}
Vol[current_index] = mass/rho_i[current_index];
Minv[current_index] = glm::inverse(sumM);
}
}
}
3) is more involved and this is how i calculate the internal forces. D is a vec3 containing the three coefficients of the elasiticity matrix.
Code: Select all
//Compute internal forces
int j,k;
float d15 = Y / (1.0f + nu) / (1.0f - 2 * nu);
float d16 = (1.0f - nu) * d15;
float d17 = nu * d15;
float d18 = Y / 2 / (1.0f + nu);
glm::vec3 D(d16, d17, d18); //Isotropic elasticity matrix D
//Calculate displacement
for(i=0;i<total_points;i++) {
U[i] = X[i]-Xi[i];
}
glm::mat3 delU = glm::mat3(1);
glm::mat3 e = glm::mat3(1); //strain
glm::mat3 s = glm::mat3(1); //stress
glm::mat3 B = glm::mat3(1);
//Calculate derivatives
for( k=0;k<numZ;k++) {
for( j=0;j<numY;j++) {
for( i=0;i<numX;i++) {
int current_index =(((k*numY)+j)*numX)+i;
glm::vec3 sumMx = glm::vec3(0);
glm::vec3 sumMy = glm::vec3(0);
glm::vec3 sumMz = glm::vec3(0);
for(int m=0;m<18;m++) {
glm::vec3 n = getNextNeighbour(m);
int ix = int(n.x+i);
int iy = int(n.y+j);
int iz = int(n.z+k);
int index = (((iz*numY)+iy)*numX)+ix;
if (index < 0)
index =0;
if (index >=total_points)
index = total_points-1;
glm::vec3 Xij = X[index]-X[current_index];
float r = glm::length(Xij);
float h = 3*r;
glm::vec3 Uij = U[index]-U[current_index];
sumMx += Uij.x*Xij*kernel(r,h);
sumMy += Uij.y*Xij*kernel(r,h);
sumMz += Uij.z*Xij*kernel(r,h);
}
glm::vec3 delUx = Minv[current_index]*sumMx;
glm::vec3 delUy = Minv[current_index]*sumMy;
glm::vec3 delUz = Minv[current_index]*sumMz;
delU = glm::mat3(delUx, delUy, delUz);
glm::mat3 delUT = glm::transpose(delU);
e = 0.5*( delU + delUT + delUT*delU);
s[0][0] = D.x*e[0][0]+D.y*e[1][1]+D.y*e[2][2];
s[1][1] = D.y*e[0][0]+D.x*e[1][1]+D.y*e[2][2];
s[2][2] = D.y*e[0][0]+D.y*e[1][1]+D.x*e[2][2];
s[0][1] = D.z*e[0][1];
s[1][2] = D.z*e[1][2];
s[2][0] = D.z*e[2][0];
s[0][2] = s[2][0];
s[1][0] = s[0][1];
s[2][1] = s[1][2];
B = - 2*Vol[current_index]*(delU + I)*s*Minv[current_index];
for(int m=0;m<18;m++) {
gm::vec3 n = getNextNeighbour(m);
int ix = int(n.x+i);
int iy = int(n.y+j);
int iz = int(n.z+k);
int index = (((iz*numY)+iy)*numX)+ix;
if (index < 0)
index =0;
if (index >=total_points)
index = total_points-1;
glm::vec3 Xij = X[index]-X[current_index];
float r = glm::length(Xij);
float h = 3;
F[index] += B*Xij*kernel(r,h);
F[current_index] -= B*Xij*kernel(r,h);
}
}
}
}Thanks,
Mobeen