Code:

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);        }    }}

Code:

//Compute internal forcesint 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 displacementfor(i=0;i<total_points;i++) {   U[i] = X[i]-Xi[i]; }glm::mat3 delU = glm::mat3(1);  glm::mat3 e = glm::mat3(1);             //strainglm::mat3 s = glm::mat3(1);//stressglm::mat3 B = glm::mat3(1);//Calculate derivativesfor( 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);            }        }    }}