### Ten node tetrahedron shape funtion

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]
```