how to Implement linear tetrahedral elements for finite element computations?

Question

I am trying to implement 3D tetrahedral elements in my finite element code (which works fine for linear triangles and quadrangles in 2D). But my simulations are crashing with tetrahedral elements. My implementation details in C++ are as follows:

Shape functions and their derivatives w.r.t. $\xi$, $\eta$ and $\zeta$

$N_0 = \xi$
$N_1 = \eta$
$N_2 = \zeta$
$N_3 = 1-\xi-\eta-\zeta$

$N_{0,\xi}$ = 1.0
$N_{1,\xi}$ = 0.0
$N_{2,\xi}$ = 0.0
$N_{3,\xi}$ = -1.0

$N_{0,\eta}$ = 0.0
$N_{1,\eta}$ = 1.0
$N_{2,\eta}$ = 0.0
$N_{3,\eta}$ = -1.0

$N_{0,\zeta}$ = 0.0
$N_{1,\zeta}$ = 0.0
$N_{2,\zeta}$ = 1.0
$N_{3,\zeta}$ = -1.0

For volume integrals

Gauss quadrature points for numerical integration:

(0.5854101966249685, 0.1381966011250105, 0.1381966011250105)
(0.1381966011250105, 0.5854101966249685, 0.1381966011250105)
(0.1381966011250105, 0.1381966011250105, 0.5854101966249685)
(0.1381966011250105, 0.1381966011250105, 0.1381966011250105)

Weight of integration points

$$W = \frac{1}{24}$$

For surface integrals

face_[0] = {1,3,2} // $\xi=0$
face_[1] = {0,2,3} // $\eta=0$
face_[2] = {0,3,1} // $\zeta=0$
face_[3] = {0,1,2} // $\xi + \eta + \zeta = 1$

Gauss quadrature points

for face 0 $(0, \frac{2}{3}, \frac{1}{6}), (0, \frac{1}{6}, \frac{1}{6}), (0, \frac{1}{6}, \frac{2}{3})$

for face 1 $(\frac{2}{3}, 0, \frac{1}{6}), (\frac{1}{6}, 0, \frac{1}{6}), (\frac{1}{6}, 0, \frac{2}{3})$

for face 2 $(\frac{2}{3}, \frac{1}{6}, 0), (\frac{1}{6}, \frac{1}{6}, 0), (\frac{1}{6}, \frac{2}{3}, 0)$

for face 3 $(\frac{2}{3}, \frac{1}{6}, \frac{1}{6}), (\frac{1}{6}, \frac{2}{3}, \frac{1}{6}), (\frac{1}{6}, \frac{1}{6}, \frac{2}{3})$

Weight of integration points

$$W = \frac{1}{6}$$

Area vectors for each face

if($\xi$ == 0.0)
{
det_normal_[0] = (y3-y4)(z2-z4) - (y2-y4)(z3-z4)
det_normal_[1] = (z3-z4)(x2-x4) - (z2-z4)(x3-x4)
det_normal_[2] = (x3-x4)(y2-y4) - (x2-x4)(y3-y4)
}
else if($\eta$ == 0.0)
{
det_normal_[0] = (y1-y4)(z3-z4) - (y3-y4)(z1-z4)
det_normal_[1] = (z1-z4)(x3-x4) - (z3-z4)(x1-x4)
det_normal_[2] = (x1-x4)(y3-y4) - (x3-x4)(y1-y4)
}
else if($\zeta$ == 0.0)
{
det_normal_[0] = (z1-z4)(y2-y4) - (y1-y4)(z2-z4)
det_normal_[1] = (x1-x4)(z2-z4) - (z1-z4)(x2-x4)
det_normal_[2] = (y1-y4)(x2-x4) - (x1-x4)(y2-y4)
}
else if($\xi + \eta + \zeta$ == 1.0)
{
det_normal_[0] = (y2-y1)(z3-z1) - (y3-y1)(z2-z1)
det_normal_[1] = (z2-z1)(x3-x1) - (z3-z1)(x2-x1)
det_normal_[2] = (x2-x1)(y3-y1) - (x3-x1)(y2-y1)
}

My questions are:

1. Are Gauss integration points for face 3 correct?
2. Are my computations for Area vector are correct?

Any suggestion or comment will be a great help.

Answer 1

I figured it out. The above formulation is correct. I had a typo in my implementation. In area vector computations the last else if condition must be replaced with else condition. Mathematically $\xi + \eta + \zeta = 2/3 + 1/6 + 1/6 = 1$, but in computer decimal representation this may not be equal to one.