Ten node tetrahedron shape funtion

by Yi Zhang

A piece of Mathematica script. The Cartesian coordinates representation is transformed into natural coordinates of tetrahedron, in which a standard one is used for isoparametric form of FE. The Lagrangian quadratic shape function corresponds to the ten node elements generated in TetGen using -o2 switch.

tet = {{0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
v[{p1_, p2_, p3_, p4_}] :=
 Abs[(p2 -
     p1).((p3 - p1)\[Cross](p4 - p1))]/6;(*volumn of a tetrahedron*)
nt[{x_, y_, z_}, tet_] :=
 1/v[tet] {v[{{x, y, z}, tet[[2]], tet[[3]], tet[[4]]}],
   v[{{x, y, z}, tet[[1]], tet[[3]], tet[[4]]}],
   v[{{x, y, z}, tet[[1]], tet[[2]], tet[[4]]}],
   v[{{x, y, z}, tet[[1]], tet[[2]], tet[[3]]}]}
shape10[{x_, y_, z_}, tet_] :=
 Module[{L = nt[{x, y, z}, tet], n1, n2, n3, n4, n5, n6, n7, n8, n9,
   n10},
  n1 = (2 L[[1]] - 1) L[[1]]; n2 = (2 L[[2]] - 1) L[[2]];
  n3 = (2 L[[3]] - 1) L[[3]]; n4 = (2 L[[4]] - 1) L[[4]];
  n5 = 4 L[[1]] L[[2]]; n6 = 4 L[[3]] L[[1]]; n7 = 4 L[[4]] L[[1]];
  n8 = 4 L[[3]] L[[2]]; n9 = 4 L[[3]] L[[4]]; n10 = 4 L[[4]] L[[2]];
  {n1, n2, n3, n4, n5, n6, n7, n8, n9, n10}]

Update:another version, involving no ABS function, better suited for symbolic differentiation, results from different representation of natural coordinates.

tet0 = {{0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
nt[{x_, y_, z_}, tet] := {
  (tet[[2]] - tet[[4]])\[Cross](tet[[2]] - tet[[3]]).(tet[[2]] - {x,
      y, z})/(tet[[2]] - tet[[4]])\[Cross](tet[[2]] -
      tet[[3]]).(tet[[2]] -
     tet[[1]]), (tet[[3]] - tet[[1]])\[Cross](tet[[3]] -
      tet[[4]]).(tet[[3]] - {x, y, z})/(tet[[3]] -
      tet[[1]])\[Cross](tet[[3]] - tet[[4]]).(tet[[3]] - tet[[2]]),
  (tet[[4]] - tet[[2]])\[Cross](tet[[4]] - tet[[1]]).(tet[[4]] - {x,
      y, z})/(tet[[4]] - tet[[2]])\[Cross](tet[[4]] -
      tet[[1]]).(tet[[4]] - tet[[3]]),
  (tet[[1]] - tet[[3]])\[Cross](tet[[1]] - tet[[2]]).(tet[[1]] - {x,
      y, z})/(tet[[1]] - tet[[3]])\[Cross](tet[[1]] -
      tet[[2]]).(tet[[1]] - tet[[4]])
  }
nt[{x_, y_, z_}] := nt[{x, y, z}, tet0];
shape4[{x_, y_, z_}, tet_] := nt[{x, y, z}, tet];
shape4[{x_, y_, z_}] := shape4[{x, y, z}, tet0];
shape10[{x_, y_, z_}, tet_] := Module[
  {L = nt[{x, y, z}, tet], n1, n2, n3, n4, n5, n6, n7, n8, n9, n10},
  n1 = (2 L[[1]] - 1) L[[1]];
  n2 = (2 L[[2]] - 1) L[[2]];
  n3 = (2 L[[3]] - 1) L[[3]];
  n4 = (2 L[[4]] - 1) L[[4]];
  n5 = 4 L[[1]] L[[2]];
  n6 = 4 L[[3]] L[[1]];
  n7 = 4 L[[4]] L[[1]];
  n8 = 4 L[[3]] L[[2]];
  n9 = 4 L[[3]] L[[4]];
  n10 = 4 L[[4]] L[[2]];
  {n1, n2, n3, n4, n5, n6, n7, n8, n9, n10}
  ]
shape10[{x_, y_, z_}] := shape10[{x, y, z}, tet0]
Advertisements