# how to Implement linear tetrahedral elements for finite element computations?

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$$

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.

• Please provide more details about where your simulation crashes. Is it due to negative Jacobian or segmentation fault or something else? If it is seg fault, then it is easy to fix. If it is not seg fault, then compute the Jacobian for one element. It should match the volume of the element. By the way, you don't need four quadrature points for computing the Jacobian (or stiffness matrix) for 4-noded Tet element. One is enough. – Chenna K Oct 22 '20 at 10:59
• Some of those details can be in the form of mathematics. This site uses MathJax to render equations. – nicoguaro Oct 22 '20 at 17:22
• @ChennaK thank you for replying. I do not get segmentation fault. It is simply the simulation diverges as soon as I launch it. The same simulation works just fine with hexahedral elements. So surely the mistake is in the implementation of tetrahedral element. I have checked the Jacobean matrix for an element. It is positive and equal to the 6*volume. My doubt is specifically for computation of surface integral on face made by nodes 012; as the triangles 013, 123 and 230 are themselves master triangles while 012 is not. What will be the integration points and weights for the face 012? – Deepak Garg Oct 24 '20 at 15:27
• @nicoguaro I have put the question with mathematical notations now. – Deepak Garg Oct 24 '20 at 15:28

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.