Page 335 - The Combined Finite-Discrete Element Method

318 ALGORITHM IMPLEMENTATION

.............continued from previous listing..............
for(i=0;i<4;i++)/* Nodal Forces */
{ j=i+1; if(j>3)j=0;
k=j+1; if(k>3)k=0;
l=k+1; if(l>3)l=0;
nx=((d1nccy[(i2elto[k][ielem])]-d1nccy[(i2elto[j][ielem])])*
(d1nccz[(i2elto[l][ielem])]-d1nccz[(i2elto[j][ielem])])-
(d1nccy[(i2elto[l][ielem])]-d1nccy[(i2elto[j][ielem])])*
(d1nccz[(i2elto[k][ielem])]-d1nccz[(i2elto[j][ielem])]))/R6;
ny=((d1nccz[(i2elto[k][ielem])]-d1nccz[(i2elto[j][ielem])])*
(d1nccx[(i2elto[l][ielem])]-d1nccx[(i2elto[j][ielem])])-
(d1nccx[(i2elto[k][ielem])]-d1nccx[(i2elto[j][ielem])])*
(d1nccz[(i2elto[l][ielem])]-d1nccz[(i2elto[j][ielem])]))/R6;
nz=((d1nccx[(i2elto[k][ielem])]-d1nccx[(i2elto[j][ielem])])*
(d1nccy[(i2elto[l][ielem])]-d1nccy[(i2elto[j][ielem])])-
(d1nccy[(i2elto[k][ielem])]-d1nccy[(i2elto[j][ielem])])*
(d1nccx[(i2elto[l][ielem])]-d1nccx[(i2elto[j][ielem])]))/R6;
d1nmct[(i2elto[i][ielem])]=d1nmct[(i2elto[i][ielem])]+
dpero*voli/R6;
if((i==0)||(i==2))
{ d1nfcx[(i2elto[i][ielem])]=d1nfcx[(i2elto[i][ielem])]+
(T[0][0]*nx+T[0][1]*ny+T[0][2]*nz);
d1nfcy[(i2elto[i][ielem])]=d1nfcy[(i2elto[i][ielem])]+
(T[1][0]*nx+T[1][1]*ny+T[1][2]*nz);
d1nfcz[(i2elto[i][ielem])]=d1nfcz[(i2elto[i][ielem])]+
(T[2][0]*nx+T[2][1]*ny+T[2][2]*nz);
}
else
{ d1nfcx[(i2elto[i][ielem])]=d1nfcx[(i2elto[i][ielem])]-
(T[0][0]*nx+T[0][1]*ny+T[0][2]*nz);
d1nfcy[(i2elto[i][ielem])]=d1nfcy[(i2elto[i][ielem])]-
(T[1][0]*nx+T[1][1]*ny+T[1][2]*nz);
d1nfcz[(i2elto[i][ielem])]=d1nfcz[(i2elto[i][ielem])]-
(T[2][0]*nx+T[2][1]*ny+T[2][2]*nz);
}} }}}
Listing 10.42 Calculation of nodal forces.

In Listing 10.42, the nodal forces are calculated using surface tractions, as explained in
Chapter 4. The surface traction force on each of the surfaces of the tetrahedron is replaced
by equivalent forces on three nodes belonging to the particular surface. One-third of the
traction force is applied to each node. Surface normals are obtained as a cross product
of the corresponding edges (i.e. base vectors) in the deformed configuration. This cross
product is divided by 2, because the surfaces of the tetrahedron are triangles (half of the
area of the corresponding rectangle). Thus, division by 2 is followed by division by 3
(one-third of the traction force is assigned to each node), which is equivalent to division
by 6. That is why normals are initially divided by R6, which is a constant defined in
frame.h andisequal to6.