10

The Quadratic
Tetrahedron

101

Chapter 10: THE QUADRATIC TETRAHEDRON

TABLE OF CONTENTS
Page

10.1
10.2

10.3

10.4

10.5
10.6
10.7
10.
10.

Introduction
. . . . . . . . . . . . . . . . . . . . .
The Ten-Node Tetrahedron
. . . . . . . . . . . . . . .
10.2.1
Element Definition
. . . . . . . . . . . . . . .
10.2.2
Constant Versus Variable Metric
. . . . . . . . . .
10.2.3
Hierarchical Geometric Formulation . . . . . . . . . .
10.2.4
*Geometric Visualization . . . . . . . . . . . . .
Partial Derivative Calculations
. . . . . . . . . . . . . .
10.3.1
Arbitrary Iso-P Tetrahedron . . . . . . . . . . . .
10.3.2
Specialization To Quadratic Tetrahedron
. . . . . . . .
10.3.3
Shape Function Module
. . . . . . . . . . . . .
The Quadratic Tetrahedron As Iso-P Elasticity Element
. . . . .
10.4.1
Displacements, Strains, Stresses . . . . . . . . . . .
10.4.2
Numerically Integrated Stiffness Matrix
. . . . . . . .
Element Stiness Matrix
. . . . . . . . . . . . . . . .
*The Jacobian Matrix Paradox
. . . . . . . . . . . . . .
*Gauss Integration Rules For Tetrahedra
. . . . . . . . . .
Notes and Bibliography
. . . . . . . . . . . . . . . . .
Exercises . . . . . . . . . . . . . . . . . . . . . .

102

103
103
103
104
105
106
109
109
1012
1013
1013
1013
1015
1016
1017
1020
1024
1025

10.2

THE TEN-NODE TETRAHEDRON

10.1. Introduction
The four-node, linear tetrahedron introduced in the previous Chapter is easy to formulate and
implement. It is the workhorse of most 3D mesh generators. It performs poorly, however, for stress
analysis in structures and solid mechanics. The quadratic tetrahedron covered in this Chapter is
significantly better in that regard. In addition, it retains the geometry favored in 3D mesh generation,
and allows for curved faces and sides. The element has several drawbacks: more complicated
formulation, increased formation time (about 15 to 30 times that of the linear tetrahedron) and
increasing connectivity in the master stiffness matrix.
10.2. The Ten-Node Tetrahedron
The ten-node tetrahedron, also called quadratic tetrahedron, is shown in Figure 10.1. In programming contexts it will be called Tet10. This is the second complete-polynomial member of the
isoparametric tetrahedron family. As noted above, for stress calculations it behaves significantly
better than the four-node (linear) tetrahedron.1 It takes advantage of the existence of fully automatic
tetrahedral meshers, something still missing for bricks.
The element has four corner nodes with local numbers 1 through 4, which must be traversed
following the same convention explained for the four node tetrahedron in the previous Chapter. It
has six side nodes, with local numbers 5 through 10. Nodes 5,6,7 are located on sides 12, 23 and
31, respectively, while nodes 8,9,10 are located on sides 14, 24, and 34, respectively. The side
nodes are not necessarily placed at the midpoints of the sides but may deviate from those locations,
subjected to positive-Jacobian-determinant constraints. Each element face is defined by six nodes.
These do not necessarily lie on a plane, but they should not deviate too much from it. This freedom
allows the element to have curved sides and faces; see Figure 10.1(b).
The tetrahedral natural coordinates 1 through 4 are introduced in a fashion similar to that described
in the previous Chapter for the linear tetrahedron. If the element has variable metric (VM), as
discussed below, there are some differences:
1.

i = constant is not necessarily the equation of a plane;

2.

Geometric definitions in terms of either distance or volume ratios are no longer valid.

10.2.1. Element Definition

The definition of the quadratic tetrahedron as an isoparametric element is

1
1
x x1

y y1

z = z1

u x u x1

uy
u y1
uz
u z1
1

1
x2
y2
z2
u x2
u y2
u z2

1
x3
y3
z3
u x3
u y3
u z3

1
x4
y4
z4
u x4
u y4
u z4

1
x5
y5
z5
u x5
u y5
u z5

N1e
...
1
N2e
. . . x10

N3e
. . . y10
N4e
. . . z 10
N5e
. . . u x10

. . . u y10 ...
. . . u z10
Ne

(10.1)

10

A mesh of quadratic tetrahedra exhibits wider connectivity than a corresponding mesh of linear tetrahedra with the same
number of nodes. Consequently the solution of the master stiffness equations will be costlier for the former.

103

Chapter 10: THE QUADRATIC TETRAHEDRON

(a)

(b)

10
5

10
5

1
4

1
4
3

Figure 10.1. The quadratic (ten-node) tetrahedron. (a): element with planar faces and side
nodes located at side midpoints; (b): element with curved faces and sides.

The conventional (non-hierarchical) shape functions are given by

N1e = 1 (21 1), N2e = 2 (22 1), N3e = 3 (23 1), N4e = 4 (24 1),
e
= 43 4 .
N5e = 41 2 , N6e = 42 3 , N7e = 43 1 , N8e = 41 4 , N9e = 42 4 , N10
(10.2)
These shape functions can be readily built by inspection, using the product guessing technique
explained in Chapter 18 of IFEM. They closely resemble those of the six node quadratic triangle
discussed in that Chapter. If the element is curved (that is, the six nodes that define each face are
not on a plane and/or the side nodes are not located at the side midpoints), constant tetrahedron
coordinates no longer define planes, but form a curvilinear system. More on that topic later.
10.2.2. Constant Versus Variable Metric
A constant metric (abbreviation: CM) tetrahedron is mathematically defined as one in which the
Jacobian determinant J over the element domain is constant.2 Otherwise the element is said to
be of variable metric (abbreviation: VM). The linear tetrahedron is always a CM element because
J = 6V is constant over it. The quadratic tetrahedron is CM if and only if the six midside modes
are collocated at the midpoints between adjacent corners. Mathematically:
x5 = 12 (x1 + x2 ),

x6 = 12 (x2 + x3 ),

...

x10 = 12 (x3 + x4 ),

y5 = 12 (y1 + y2 ),

y6 = 12 (y2 + y3 ),

...

y10 = 12 (y3 + y4 ),

z 5 = 12 (z 1 + z 2 ),

z 6 = 12 (z 2 + z 3 ),

...

z 10 = 12 (z 3 + z 4 ).

(10.3)

If the conditions (10.3) are imposed upon the first four rows of the iso-P definition (10.1), that
portion reduces to


1
1
1 1 1 1
x
x
x

x
x
1
2
3
4 2
(10.4)
.
=
y1 y2 y3 y4
3
y
z1 z2 z3 z4
4
z
2

The expression of J for an arbitrary iso-P tetrahedron is worked out in 10.3.1.

104

10.2
4

THE TEN-NODE TETRAHEDRON

4
8

9 10

8
10

7
1

Figure 10.2. Illustrating differences between constant metric (CM) and variable metric (VM) tetrahedra.
(a) CM quadratic tetrahedron; (b) VM quadratic tetrahedron with the same corner locations as the CM
one; (c) the isocoordinate surfaces 4 = 0, 4 = 14 , 4 = 12 , 4 = 34 and 4 = 49
50 for the VM tetrahedron.
Plots produced by module PlotTet10Element listed in Figure 10.4.

This is exactly the geometric definition (9.10) of the linear tetrahedron. Therefore J is the same at
all points, and expressed in terms of the corner coordinates by (9.3) or (9.4). Consequences: the
four faces are planar, and the six sides straight.3 If at least one midside node deviates from (10.3),
J will be a function of position and the element is VM.
Figure 10.2(a,b) visually contrasts two quadratic tetrahedra of constant and variable metric, respectively. Both tetrahedra have the same corner locations. For the VM one, Figure 10.2(c) pictures
the five isocoordinate surfaces 4 = 0, 4 = 14 , 4 = 12 , 4 = 34 , and 4 = 49
(a plot of 4 = 1 would
50
be invisible, as it collapses to a point). This plot illustrates the fact that the tetrahedral coordinates
i form a system of curvilinear coordinates if the element has VM geometry. These plots were
produced by the module PlotTet10Element described in the next subsection.
Why is CM versus VM important? Implementation side effects. If the tetrahedron is CM,
1.

Computation of partial derivatives such as i / x is easy because (9.18) can be reused.

2.

Analytical integration (over volume, faces or sides) is straightforward, as long as integrands

are polynomials in tetrahedral coordinates. See 9.1.10 about how to do it.

For VM elements the computation of partial derivatives is considerably more laborious.4 In addition,
numerical integration is necessary. So a key implementation question is: are VM tetrahedra
allowed? In the sequel we shall assume that to be the case. Consequently the more general
formulation is presented.
Remark 10.1. Elements whose geometric definition is of lower order (in the sense of polynomial variation

in natural coordinates) than the displacement interpolation are called subparametric in the FEM literature.
Pre-1967 triangular and tetrahedral elements were subparametric until isoparametric models came along. The
opposite case: elements whose geometry is of higher order (again, polynomial-wise) than the displacement
field interpolation are called superparametric. These are comparatively rare.
3

The converse is not true: even if a quadratic tetrahedron has planar faces and straight sides, it is not necessarily CM, since
deviations from (10.3) could be such that midside nodes remain on the straight line passing through adjacent corners.

Reason: if one tries to directly express the i in terms of {x, y, z} from the first four rows of the Iso-P definition (10.1),
it would entail solving 4 quadratic polynomial equations for 4 coupled unknowns. A closed form solution comparable
to (9.11) does not generally exist. But such a solution is possible for the differentials, as shown in 10.3.1.

105

Chapter 10: THE QUADRATIC TETRAHEDRON

10.2.3. Hierarchical Geometric Formulation

An alternative geometric descrition of the tetrahedral element uses deviations from midpoint locations to establish the position of the midside nodes:
x5 = 12 (x1 + x2 ) + x5 ,

y5 = 12 (y1 + y2 ) + y5 ,

z 5 = 12 (z 1 + z 2 ) + z 5 ,

x6 = 12 (x2 + x3 ) + x6 ,

y6 = 12 (y2 + y3 ) + y6 ,

z 6 = 12 (z 2 + z 3 ) + z 6 ,

...
x10 = 12 (x3 + x4 ) + x10 ,

y10 = 12 (y3 + y4 ) + y10 ,

(10.5)

z 10 = 12 (z 3 + z 4 ) + z 10 .

in which xi , yi , and z i are the midpoint coordinate deviations. Rewriting the geometry
definition from (10.1) in terms of these deviations yields

1
1
x x1
=
y
y1
z1
z

1
x2
y2
z2

1
x3
y3
z3

1
x4
y4
z4

0
x5
y5
z 5

2
...
0
3
. . . x10
4e
. . . y10
N5
. . . z 10
..
.

(10.6)

e
N10

The four corner shape functions change from i (2i 1) to just i , which are the shape functions
of the linear tetrahedron. The shape functions associated with the coordinate deviations remain
unchanged. This reformulation is called a geometry-hierarchical element. It has certain implementation merits. The chief advantage: for a CM element all midpoint deviations vanish, allowing
immediate regression to linear tetrahedron formulas. In addition the expression of the Jacobian
derived in 10.3 is simpler.
If midside nodal displacements are similarly adjusted, the model is called a iso-P hierarchical
element. It offers certain advantages (as well as shortcomings) that are not elaborated upon here.
10.2.4. *Geometric Visualization
Picturing VM tetrahedra is useful to visually understand concepts such as curved sides and faces, as well as
tetrahedral coordinate isosurfaces. The plots shown in Figure 10.3 were produced by the Mathematica module
PlotTet10Element listed in Figure 10.4. That module is invoked as
PlotTet10Element[xyztet,plotwhat,Nsub,aspect,view,box]

(10.7)

The arguments are:

xyztet

Tetrahedon {x, y, z} nodal coordinates stored node-by-node as a two-dimensional list:

{ { x1,y1,z1 },{ x2,y2,z2 },{ x3,y3,z3 }, ... { x10,y10,z10 } }.

plotwhat

A list specifying what is to be plotted. List items must be of one of three forms:
{ i,c }, where i = 1, 2, 3, 4, and c a real number c [0, 1]: plot surface i = c.
10*i+j, where i, j = 1, 2, 3, 4: plot side joining corners i and j,
n, where n = 1, 2, . . . 10: plot element node n as dark dot and write its number nearby.
The plotwhat configurations used to generate the plots of Figure 10.4 are given below.

106

10.2 THE TEN-NODE TETRAHEDRON

(a)

(b)

(c)
4

8
10
5

1
7

6
3

(d)

(e)

(f)

8
10
5

7
6

Figure 10.3. Sample plots for a variable metric (VM) quadratic tetrahedron, produced by Mathematica module
PlotTet10Element. (a): sides and faces; (b) sides only; (c) sides and nodes; (d) sides, faces and nodes; (e) the
five surfaces 2 = 0 (face 2), 2 = 14 , 2 = 12 , 2 = 34 , and 2 = .98 (2 = 1 would produce only a point as it
reduces to node 2.); (e) the four surfaces 1 = 1/4 for i = 1, 2, 3, 4 showing intersection at ( 14 , 14 , 14 , 14 ) (which is
not necessarily the centroid if tetrahedron is VM). For plot input data see text.

Nsub

Number of subdivisions along sides used to triangulate faces. If tetrahedron is constant

metric (CM), Nsub=1 is enough, but if it has curved faces a pleasing value of Nsub should be
determined by trial and error. The plots of Figure 10.4 use Nsub=16.

aspect

Specifies the z aspect ratio of the plot box frame. Normally 1.

view

The coordinates of a point in space from which the objects plotted are to be viewed. Used to
set the option ViewPoint->view in Mathematica. Typically determined by trial and error.
The plots of Figure 10.4 were done with view = { -2,-2,1 }.

box

A logical flag. Set box to True to get a coordinates aligned box-frame drawn around the
tetrahedron, as illustrated by the plots of Figure 10.4. If set to False, no box is drawn.

The function PlotTet10Element does not return a value.

The tetrahedron geometry plotted in Figure 10.3 is defined by the node coordinates
xyztet = { { 1.,0.,0. },{ 0.,1.,0. }, { 0.,0.,0. },{ 0.5,0.5,1. },{ 0.5,0.65,0. },
{ -0.1,0.4,0. },{ 0.5,0.1,0. }, { 0.85,0.25,0.6 },{ 0.35,0.85,0.5 },{ 0.15,0.25,0.6 } };
Further Nsub=16, view={ -2,-2,1 }, aspect=1, and box=True for all six plots, which are produced by
(a) PlotTet10Element[xyztet,{{1,0},{2,0},{3,0},{4,0},12,23,31,14,24,34},
Nsub,aspect,view,box];
(b) PlotTet10Element[xyztet,{12,23,31,14,24,34},Nsub,aspect,view,box];

107

Chapter 10: THE QUADRATIC TETRAHEDRON

PlotTet10Element[xyztet_,plotwhat_,Nsub_,aspect_,view_,box_]:=
Module[{Ne,m,i,k,x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,y1,y2,y3,y4,y5,
y6,y7,y8,y9,y10,z1,z2,z3,z4,z5,z6,z7,z8,z9,z10, 1, 2, 3, 4,c,d,
c1,c2,c3,Ne1,Ne2,Ne3,xyz1,xyz2,xyz3,i1,i2,i3,nplot,iplot,what,
i,Ni=3*Nsub,xyzn,exyz,nlab,style, lines={},polys={},sides={},
nodlab={},nodcir={},DisplayChannel,modname="PlotTet10Element: "},
Ne[{ 1_, 2_, 3_, 4_}]:={ 1*(2* 1-1), 2*(2* 2-1), 3*(2* 3-1),
4*(2* 4-1),4* 1* 2,4* 2* 3,4* 3* 1,4* 1* 4,4* 2* 4,4* 3* 4};
{{x1,y1,z1},{x2,y2,z2},{x3,y3,z3},{x4,y4,z4},{x5,y5,z5},{x6,y6,z6},
{x7,y7,z7},{x8,y8,z8},{x9,y9,z9},{x10,y10,z10}}=Take[xyztet,10];
xn={x1,x2,x3,x4,x5,x6,x7,x8,x9,x10};
yn={y1,y2,y3,y4,y5,y6,y7,y8,y9,y10};
zn={z1,z2,z3,z4,z5,z6,z7,z8,z9,z10};
style=TextStyle->{FontFamily->"Times",FontSize->15,
FontWeight->"Bold"};
nplot=Length[plotwhat]; If [nplot<=0,Return[]];
For [iplot=1,iplot<=nplot,iplot++, what=plotwhat[[iplot]];
If [Length[what]==2, {i, i}=what;
If [Count[{1,2,3,4},i]<=0, Print[modname,
"Bad face index in plotwhat"]; Return[]]; c=(1- i)/Ni;
For [i1=1,i1<=Ni,i1++, For [i2=1,i2<=Ni-i1,i2++,
i3=Ni-i1-i2; If [i3<=0, Continue[]]; d=0;
If [Mod[i1+2,3]==0&&Mod[i2-1,3]==0, d= 1];
If [Mod[i1-2,3]==0&&Mod[i2+1,3]==0, d=-1];
If [d==0, Continue[]]; k=i-1;
c1=RotateRight[N[{ i,(i1+d+d)*c,(i2-d)*c,(i3-d)*c}],k];
c2=RotateRight[N[{ i,(i1-d)*c,(i2+d+d)*c,(i3-d)*c}],k];
c3=RotateRight[N[{ i,(i1-d)*c,(i2-d)*c,(i3+d+d)*c}],k];
Ne1=Ne[c1]; xyz1={xn.Ne1,yn.Ne1,zn.Ne1};
Ne2=Ne[c2]; xyz2={xn.Ne2,yn.Ne2,zn.Ne2};
Ne3=Ne[c3]; xyz3={xn.Ne3,yn.Ne3,zn.Ne3};
AppendTo[polys,Polygon[{xyz1,xyz2,xyz3}]];
AppendTo[lines,Line[{xyz1,xyz2,xyz3,xyz1}]];
]]; Continue[];
]; m=what;
If [m>0&&m<=10, nlab=ToString[m];
xyzn={xn[[m]],yn[[m]],zn[[m]]}; exyz={0.02,0,.06};
AppendTo[nodlab,Graphics3D[Text[nlab,xyzn+exyz,style]]];
AppendTo[nodcir,Graphics3D[Point[xyzn]]]; Continue[]];
If [Count[{12,21,13,31,14,41,23,32,24,42,34,43},m]<=0,
Print[modname," Bad side id in plotwhat"]; Return[]];
For [i1=1,i1<=Nsub,i1++, i2=i1+1;
1=N[(i1-1)/Nsub]; 2=N[(i2-1)/Nsub];
If [m==12||m==21,c1={ 1,1- 1,0,0}; c2={ 2,1- 2,0,0}];
If [m==13||m==31,c1={ 1,0,1- 1,0}; c2={ 2,0,1- 2,0}];
If [m==14||m==41,c1={ 1,0,0,1- 1}; c2={ 2,0,0,1- 2}];
If [m==23||m==32,c1={0, 1,1- 1,0}; c2={0, 2,1- 2,0}];
If [m==24||m==42,c1={0, 1,0,1- 1}; c2={0, 2,0,1- 2}];
If [m==34||m==43,c1={0,0, 1,1- 1}; c2={0,0, 2,1- 2}];
Ne1=Ne[c1]; xyz1={xn.Ne1,yn.Ne1,zn.Ne1};
Ne2=Ne[c2]; xyz2={xn.Ne2,yn.Ne2,zn.Ne2};
AppendTo[sides,Line[{xyz1,xyz2}]]];
];
DisplayChannel=$DisplayFunction;
If [$VersionNumber>=6.0, DisplayChannel=Print];
Show[Graphics3D[RGBColor[1,0,0]],Graphics3D[polys],
Graphics3D[AbsoluteThickness[1]],Graphics3D[lines],
Graphics3D[RGBColor[0,0,0]],Graphics3D[AbsoluteThickness[2]],
Graphics3D[sides],
Graphics3D[AbsolutePointSize[6]],nodcir,nodlab,
PlotRange->All,ViewPoint->view,
AmbientLight->GrayLevel[0],BoxRatios->{1,1,aspect},
Boxed->box,DisplayFunction->DisplayChannel];
ClearAll[polys,lines,sides,nodlab,nodcir]];

Figure 10.4. Quadratic tetrahedron plotting module. See text for usage.

108

10.3

PARTIAL DERIVATIVE CALCULATIONS

(c) PlotTet10Element[xyztet,{12,23,31,14,24,34,1,2,3,4,5,6,7,8,9,10},
Nsub,aspect,view,box];
(d) PlotTet10Element[xyztet,{{1,0},{2,0},{3,0},{4,0},12,23,31,14,24,34,
1,2,3,4,5,6,7,8,9,10},Nsub,aspect,view,box];
(e) PlotTet10Element[xyztet,{{1,0},{1,1/4},{1,1/2},{1,3/4},{1,.98},
1,12,23,31,14,24,34},Nsub,aspect,view,box];
(f) PlotTet10Element[xyztet,{{1,1/4},{2,1/4},{3,1/4},{4,1/4},
12,23,31,14,24,34},Nsub,aspect,view,box];
See legend of Figure 10.3 for description of what the plots show.
Implementation Considerations. Several plot options and details are built in the module; for example lighting,
face colors, node number fonts, line thicknesses, etc. Those could be externally controlled with additional
arguments, but it is just as easy to go in and modify the module. The Show statement contains a specification,
namely, DisplayFunction->DisplayChannel, that is Mathematica version dependent.

10.3. Partial Derivative Calculations

The problem considered in this section is: given a smooth scalar function F(1 , 2 , 3 , 4 ) expressed
in tetrahedral coordinates and interpolated by the shape functions in terms of nodal values, find
F/ x, F/ y, and F/z at any point in the element. This requires finding the partial derivatives
i / x, i / y, and i /z in order to apply the chain rule. The following derivation assumes that
the tetrahedron has variable metric (VM). Simplifications for the constant metric (CM) case are
noted as appropriate. All stated derivatives are assumed to exist.
The derivation given in 10.3.1 is general in the sense that it applies to any VM, isoparametric
tetrahedral element with n nodes. The results are specialized to the quadratic tetrahedron in 10.3.2.
10.3.1. Arbitrary Iso-P Tetrahedron
Consider a generic scalar function F(1 , 2 , 3 , 4 ) that is interpolated as
Ne
1

F = [ F1

F2

F3

F4

...

N2e
e

Fn ]
.. = Fk Nk .
.
Nne

(10.8)

Here Fk denotes the value of F(1 , 2 , 3 , 4 ) at the k th node, while Nke = Nke (i ) are the appropriate
shape functions. The summation convention applies to the last expression over k = 1, 2, , . . . n.
Symbol F may stand for 1, x, y, z, u x , u y or u z in the isoparametric representation, or other
element-varying quantities such as temperature or body force components. On taking partials with
respect to x, y and z, and applying the chain rule twice we get


Nk
Nk i
Nk 2
Nk 3
Nk 4
Nk 1
F
= Fk
= Fk
+
+
+
= Fk
,
x
x
1 x
2 x
3 x
4 x
i x


Nk
Nk i
Nk 1
Nk 2
Nk 3
Nk 4
F
= Fk
= Fk
+
+
+
= Fk
, (10.9)
y
y
1 y
2 y
3 y
4 x
i y


Nk
Nk i
Nk 2
Nk 3
Nk 4
F
Nk 1
= Fk
= Fk
+
+
+
= Fk
,
z
y
1 z
2 z
3 z
4 x
i z
109

Chapter 10: THE QUADRATIC TETRAHEDRON

in which sums run over k = 1, 2, . . . n and i = 1, 2, 3, 4, and element superscripts on the shape
functions have been suppressed for clarity. In matrix form:

F
x
F

y
F
z

1
x

= 1
y
1
z

2
x
2
y
2
z

3
x
3
y
3
z

4
x
4
y
4
z

Nk
Fk

1
Nk
Fk

F Nk
k 3

Nk
Fk
4

(10.10)

Transpose both sides of (10.10), and switch the left and right hand sides to get

Nk
Fk

Nk
Fk

1
x

 2
Nk x
Fk

4
3
x

4
x

Nk
Fk

1
y
2
y
3
y
4
y

1
z
2
z
3
z
4
z

F
= x

F
y

F . (10.11)
z

Next, set F to x, y and z, stacking results row-wise. To make the coefficient matrix square,
differentiate both sides of 1 +2 +3 +4 = 1 with respect to x, y and z, and insert as first row:

x Nk
k 1

N
k
yk
1

Nk
z k

Nk
xk
2

N
yk k
2

N
z k k
2

Nk
xk
3

N
yk k
3

N
z k k
3

Nk
xk
4

N
yk k
4

N
z k k
4

1
x
2
x

x
4
x

1
y
2
y
3
y
4
y

1
z
2
z
3
z
4
z

1
x
x

x
=
y

x
z
x

1
y
x
y
y
y
z
y

1
z
x
z
y
z
z
z

. (10.12)

But x/ x = y/ y = z/z = 1, and 1/ x = 1/ y = 1z = x/ y = x/z = y/ x =

y/z = z/ x = z/ y = 0, because x, y and z are independent coordinates. Thus

1
x Nk
k 1

N
k
yk
1

Nk
z k
1

1
Nk
xk
2
Nk
yk
2
Nk
z k
2

1
Nk
xk
3
Nk
yk
3
Nk
z k
3

1
Nk
xk
4
Nk
yk
4
Nk
z k
4

1
x
2
x

x
4
x

1
y
2
y
3
y
4
y

1
z
2
z
3
z
4
z

1
=

0
0
1
0

0
0
.
0
1

(10.13)

This can be expressed in compact matrix notation as

J P = Iaug .
1010

(10.14)

10.3

PARTIAL DERIVATIVE CALCULATIONS

Here P is the 4 3 matrix of i Cartesian partial derivatives, Iaug is the order 3 identity matrix
augmented with a zero first row, and
1

J = Jx1
Jy1
Jz1

1
Jx2
Jy2
Jz2

1
Jx3
Jy3
Jz3

1
x k Nk

Jx4 = 1
Nk
yk
Jy4
1
Jz4
Nk
z k
1

1
Nk
xk
2

N
yk k
2

N
z k k
2

1
Nk
xk
3

N
yk k
3

N
z k k
3

1
Nk
xk
4

N
yk k
4

N
z k k
4

(10.15)

is called the Jacobian matrix. To express its determinant and inverse it is convenient to introduce
the abbreviations
Jxi j = Jxi Jx j ,

Jyi j = Jyi Jy j ,

Jzi j = Jzi Jz j ,

i, j = 1, . . . 4.

(10.16)

With the abbreviations (10.16), the Jacobian determinant J = det(J) can be compactly stated as
J = Jx21 (Jy23 Jz34 Jy34 Jz23 )+ Jx32 (Jy34 Jz12 Jy12 Jz34 )+ Jx43 (Jy12 Jz23 Jy23 Jz12 ). (10.17)
If J = 0, explicit inversion gives

J1

6V01
1 6V
= 02
J 6V03
6V04

The solution P = J1 Iaug

1
x
2
x
P=
3

x
4
x

1
y
2
y
3
y
4
y

1
z
2
z
3
z
4
z

Jy42
Jy13
.
Jy24
Jy31
(10.18)
of (10.14) for the i partials is supplied by the last 3 columns of J1 :

Jy42
Jy31
Jy24
Jy13

Jz32
Jz43
Jz14
Jz21

Jy32
Jy34
Jy14
Jy12

Jz42
Jz13
Jz24
Jz31

Jx32
Jx43
Jx14
Jx21

Jz42
Jz31
Jz24
Jz13

Jx42
Jx13
Jx24
Jx31

Jz32
Jz34
Jz14
Jz12

Jx42
Jx31
Jx24
Jx13

Jy32
Jy43
Jy14
Jy21

Jx32
Jx34
Jx14
Jx12

Jy42

Jy31
= 1

J Jy24

Jy13

Jz32 Jy32
Jz43 Jy34
Jz14 Jy14
Jz21 Jy12

Jz42
Jz13
Jz24
Jz31

Jx32
Jx43
Jx14
Jx21

Jz42 Jx42
Jz31 Jx13
Jz24 Jx24
Jz13 Jx31

Jz32
Jz34
Jz14
Jz12

Jx42 Jy32 Jx32 Jy42

Jx31 Jy43 Jx34 Jy13
.
Jx24 Jy14 Jx14 Jy24
Jx13 Jy21 Jx12 Jy31

(10.19)
With (10.19) available, the Cartesian partial derivatives of any function F(1 , 2 , 3 , 4 ) follow from
the chain rule expression (10.9). As in the previous Chapter, to facilitate use of the summation
convention it is convenient to represent entries of P using a more compact notation:
1
x
2
x
P=
3

x
4
x

1
y
2
y
3
y
4
y

1
z
2
z
3
z
4
z

a1

a2
= 1

J a3

a4

1011

b1
b2
b3
b4

c1
c2

c3
c4

(10.20)

Chapter 10: THE QUADRATIC TETRAHEDRON

Thus i / x = ai /J , i / y = bi /J , and i /z = ci /J , and from (10.9) and the abbreviations

(10.20) we get
F
Nk ai
= Fk
,
x
i J

Nk bi
F
= Fk
,
y
i J

Nk ci
F
= Fk
.
z
i J

(10.21)

where sums run over i = 1, 2, 3, 4 and k = 1, 2, . . . n. This shows that these partials are rational
functions in the triangular coordinates.
Remark 10.2. Entries of the first column of J1 are rarely of interest, but if necessary they can be obtained

explicitly from
6V01 = Jx2 (Jy3 Jz4 Jy4 Jz3 ) + Jx3 (Jy4 Jz2 Jy2 Jz4 ) + Jx4 (Jy2 Jz3 Jy3 Jz2 ),
6V02 = Jx1 (Jy4 Jz3 Jy3 Jz4 ) + Jx3 (Jy1 Jz4 Jy4 Jz1 ) + Jx4 (Jy3 Jz1 Jy1 Jz3 ),
6V03 = Jx1 (Jy2 Jz4 Jy4 Jz2 ) + Jx2 (Jy4 Jz1 Jy1 Jz4 ) + Jx4 (Jy1 Jz2 Jy2 Jz1 ),

(10.22)

6V04 = Jx1 (Jy3 Jz2 Jy2 Jz3 ) + Jx2 (Jy1 Jz3 Jy3 Jz1 ) + Jx3 (Jy2 Jz1 Jy1 Jz2 ).
These are the analog of (?) for the linear tetrahedron, and become identical for a CM element.
Remark 10.3. The Jacobian matrix J given in (10.15) is valid for any variable metric (VM) iso-P tetrahedron.
If the element is of constant metric (CM), J may be expected to collapse to that of the linear tetrahedron,
given in (9.3), which depends only on corner coordinates. As discussed in 10.6, those expectations are
not necessarily realized. The outcome depends on geometric definition details. What is true is that the CM
Jacobian determinant, which is explicitly given by (9.4), is always obtained regardless of definition.
Remark 10.4. The Jacobian determinant appears in the transformation of element volume integrals from

Cartesian to tetrahedral coordinates through the replacement d x d y dz = 16 J d1 d2 d3 d4 .. Thus

F(x, y, z) d
=

F(x, y, z) d x d y dz =

F(1 , 2 , 3 , 4 ) 16 J d1 d2 d3 d4 .

(10.23)

10.3.2. Specialization To Quadratic Tetrahedron

We now specialize the foregoing results to the quadratic tetrahedron, for which n = 10 and the
shape functions are given by (10.2). It is convenient to built a table of shape function derivatives

0
0
0
42 0 43 44 0
0
41 1
0
0
41 43 0
0 44 0
42 1
0
N =
(10.24)
0
0
43 1
0
0 42 41 0
0 44
0
0
0
44 1 0
0
0 41 42 43
where the {i, j} entry is Nk /i . Taking dot products with the node coordinates {xk , yk , z k } yields

/
J=4
/
/

x1 1 +x5 2 +x7 3 +x8 4

x5 1 +x2 2 +x6 3 +x9 4
x7 1 +x6 2 +x3 3 +x10 4
x8 1 +x9 2 +x10 3 +x4 4

y1 1 +y5 2 +y7 3 +y8 4

y5 1 +y2 2 +y6 3 +y9 4
y7 1 +y6 2 +y3 3 +y10 4
y8 1 +y9 2 +y10 3 +y4 4
1012

z 1 1 +z 5 2 +z 7 3 +z 8 4 T
z 5 1 +z 2 2 +z 6 3 +z 9 4

z 7 1 +z 6 2 +z 3 3 +z 10 4
z 8 1 +z 9 2 +z 10 3 +z 4 4
(10.25)

10.4

THE QUADRATIC TETRAHEDRON AS ISO-P ELASTICITY ELEMENT

in which 1 = 1 /, 2 = 2 /, 3 = 3 /, and 4 = 4 /. (Matrix J is displayed in

transposed form to fit within page width.) The entries of P follow from (10.19).
The shape function Cartesian derivatives are explicitly given by
Nn
an
bn
cn
Nn
Nn
= (4n 1) ,
= (4n 1) ,
= (4n 1) ,
x
J
y
J
z
J
ai j + a j i
bi j + b j i
ci j + c j i
Nm
Nm
Nm
=4
,
=4
,
=4
,
x
J
y
J
z
J

(10.26)

For the corner node shape functions, n = 1, 2, 3, 4. For the side node shape functions, m =
5, 6, 7, 8, 9, 10; the corresponding RHS indices being i = 1, 2, 3, 1, 2, 3, and j = 2, 3, 1, 4, 4, 4.
10.3.3. Shape Function Module
A shape function module called IsoTet10ShapeFunDer, which returns shape function values and
their derivatives with respect to {x, y, z}, is listed in Figure 10.5, The module is referenced as
{ Nfx,Nfy,Nfz,Jdet }=IsoTet10CarDer[xyztet,tcoor,numer]

(10.27)

The arguments are:

xyztet

Tetrahedral node coordinates, supplied as the two-dimensional list

{ { x1,y1,z1 },{ x2,y2,z2 },{ x3,y3,z3 }, ... { x10,y10,z10 } }.

tcoor

Tetrahedral coordinates of the point at which shape functions and their derivatives will
be evaluated, supplied as a list: {1 , 2 , 3 , 4 }

numer

A logical flag. True to carry out floating-point numeric computations, else False.

The function returns are:

Nfx

List of scaled shape function partial derivatives J Ni / x, i = 1, 2, . . . 10

Nfy

List of scaled shape function partial derivatives J Ni / y, i = 1, 2, . . . 10

Nfz

List of scaled shape function partial derivatives J Ni /z, i = 1, 2, . . . 10

Jdet

Jacobian matrix determinant J at the evaluation point.

Note that it is convenient to return, say, J Ni / x and not Ni / x to avoid zero-J exceptions inside
IsoTet10ShapeFunDer; expressions in (10.26) make this clear. Consequently the module always
returns the above data, even if Jdet is zero or negative. It is up to the calling program to proceed
upon testing whether Jdet 0.
10.4. The Quadratic Tetrahedron As Iso-P Elasticity Element
The derivation of the quadratic tetrahedron as an iso-P displacement model based on the Total
Potential Energy (TPE) principle as source variational form, largely follows the description of 9.2
in the previous Chapter. Only several implementation differences are highlighted here.
1013

Chapter 10: THE QUADRATIC TETRAHEDRON

IsoTet10ShapeFunDer[xyztet_,tcoor_,numer_]:= Module[{
x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,
x5,y5,z5,x6,y6,z6,x7,y7,z7,x8,y8,z8,
x9,y9,z9,x10,y10,z10, 1, 2, 3, 4,ncoor=xyztet,
Jx1,Jx2,Jx3,Jx4,Jy1,Jy2,Jy3,Jy4,Jz1,Jz2,Jz3,Jz4,
Jx12,Jx13,Jx14,Jx23,Jx24,Jx34,Jx21,Jx41,Jx31,Jx32,Jx42,Jx43,
Jy12,Jy13,Jy14,Jy23,Jy24,Jy34,Jy21,Jy41,Jy31,Jy32,Jy42,Jy43,
Jz12,Jz13,Jz14,Jz23,Jz24,Jz34,Jz21,Jz41,Jz31,Jz32,Jz42,Jz43,
a1,a2,a3,a4,b1,b2,b3,b4,c1,c2,c3,c4,Nfx,Nfy,Nfz,Jdet},
If [numer, ncoor=N[xyztet]];
{{x1,y1,z1},{x2,y2,z2},{x3,y3,z3},{x4,y4,z4},
{x5,y5,z5},{x6,y6,z6},{x7,y7,z7},{x8,y8,z8},
{x9,y9,z9},{x10,y10,z10}}=ncoor; { 1, 2, 3, 4}=tcoor;
Jx1=4*(x1*( 1-1/4)+x5* 2+x7* 3+x8* 4);
Jy1=4*(y1*( 1-1/4)+y5* 2+y7* 3+y8* 4);
Jz1=4*(z1*( 1-1/4)+z5* 2+z7* 3+z8* 4);
Jx2=4*(x5* 1+x2*( 2-1/4)+x6* 3+x9* 4);
Jy2=4*(y5* 1+y2*( 2-1/4)+y6* 3+y9* 4);
Jz2=4*(z5* 1+z2*( 2-1/4)+z6* 3+z9* 4);
Jx3=4*(x7* 1+x6* 2+x3*( 3-1/4)+x10* 4);
Jy3=4*(y7* 1+y6* 2+y3*( 3-1/4)+y10* 4);
Jz3=4*(z7* 1+z6* 2+z3*( 3-1/4)+z10* 4);
Jx4=4*(x8* 1+x9* 2+x10* 3+x4*( 4-1/4));
Jy4=4*(y8* 1+y9* 2+y10* 3+y4*( 4-1/4));
Jz4=4*(z8* 1+z9* 2+z10* 3+z4*( 4-1/4));
Jx12=Jx1-Jx2; Jx13=Jx1-Jx3; Jx14=Jx1-Jx4; Jx23=Jx2-Jx3;
Jx24=Jx2-Jx4; Jx34=Jx3-Jx4; Jy12=Jy1-Jy2; Jy13=Jy1-Jy3;
Jy14=Jy1-Jy4; Jy23=Jy2-Jy3; Jy24=Jy2-Jy4; Jy34=Jy3-Jy4;
Jz12=Jz1-Jz2; Jz13=Jz1-Jz3; Jz14=Jz1-Jz4;
Jz23=Jz2-Jz3; Jz24=Jz2-Jz4; Jz34=Jz3-Jz4;
Jx21=-Jx12; Jx31=-Jx13; Jx41=-Jx14; Jx32=-Jx23; Jx42=-Jx24;
Jx43=-Jx34; Jy21=-Jy12; Jy31=-Jy13; Jy41=-Jy14; Jy32=-Jy23;
Jy42=-Jy24; Jy43=-Jy34; Jz21=-Jz12; Jz31=-Jz13; Jz41=-Jz14;
Jz32=-Jz23; Jz42=-Jz24; Jz43=-Jz34;
Jdet=Jx21*(Jy23*Jz34-Jy34*Jz23)+Jx32*(Jy34*Jz12Jy12*Jz34)+Jx43*(Jy12*Jz23-Jy23*Jz12);
a1=Jy42*Jz32-Jy32*Jz42; a2=Jy31*Jz43-Jy34*Jz13;
a3=Jy24*Jz14-Jy14*Jz24; a4=Jy13*Jz21-Jy12*Jz31;
b1=Jx32*Jz42-Jx42*Jz32; b2=Jx43*Jz31-Jx13*Jz34;
b3=Jx14*Jz24-Jx24*Jz14; b4=Jx21*Jz13-Jx31*Jz12;
c1=Jx42*Jy32-Jx32*Jy42; c2=Jx31*Jy43-Jx34*Jy13;
c3=Jx24*Jy14-Jx14*Jy24; c4=Jx13*Jy21-Jx12*Jy31;
Nfx={(4* 1-1)*a1, (4* 2-1)*a2, (4* 3-1)*a3, (4* 4-1)*a4,
4*( 1*a2+ 2*a1), 4*( 2*a3+ 3*a2), 4*( 3*a1+ 1*a3),
4*( 1*a4+ 4*a1), 4*( 2*a4+ 4*a2), 4*( 3*a4+ 4*a3)};
Nfy={(4* 1-1)*b1, (4* 2-1)*b2, (4* 3-1)*b3, (4* 4-1)*b4,
4*( 1*b2+ 2*b1), 4*( 2*b3+ 3*b2), 4*( 3*b1+ 1*b3),
4*( 1*b4+ 4*b1), 4*( 2*b4+ 4*b2), 4*( 3*b4+ 4*b3)};
Nfz={(4* 1-1)*c1, (4* 2-1)*c2, (4* 3-1)*c3, (4* 4-1)*c4,
4*( 1*c2+ 2*c1), 4*( 2*c3+ 3*c2), 4*( 3*c1+ 1*c3),
4*( 1*c4+ 4*c1), 4*( 2*c4+ 4*c2), 4*( 3*c4+ 4*c3)};
Return[Simplify[{Nfx,Nfy,Nfz,Jdet}]]];

Figure 10.5. Quadratic tetrahedron shape function module IsoTet10ShapeFunDer.

10.4.1. Displacements, Strains, Stresses

The element displacement field is defined by the three components u x , u y and u z . These are
quadratically interpolated from their node values as shown in the last three rows of the element
definition (10.1).
1014

10.4

THE QUADRATIC TETRAHEDRON AS ISO-P ELASTICITY ELEMENT

The 30 1 node displacement vector is configured node-wise as

ue = [ u x1

u y1

u z1

u x2

u y2

u z2

. . . u z10 ]T .

(10.28)

It is easily shown that the resulting element is C 0 inter-element conforming. It is also complete for
the TPE functional, which has a variational index of one. Consequently, the element is consistent
in the sense discussed in Chapter 19 of IFEM.
The element strain field is strongly connected to the displacements by the strain-displacement
equations, which are listed as (9.30), (9.31), and (9.32) in the previous Chapter. The B matrix is
now 6 30: and takes the following configuration:
N1
x
0

1
0
B = N
1
J
y

N1
z

0
N1
y
0
N1
x
N1
z
0

N2
x
0

0
0
N1
z
0
N1
y
N1
x

0
N2
y
0
N2
x
N2
z
0

0
N2
y
0
N2
z

...

0
N2
z
0
N2
y
N2
x

...
...
...
...
...

N10
x
0
0
N10
y
0
N10
z

0
N10
y
0
N10
x
N10
z
0

0
0
N10
z
0
N10
y
N10
x

(10.29)

Note that this matrix varies over the element.

The stress field is defined by the stress vector (9.35). This is linked to the strains by the matrix
elasticity constitutive equations (9.34). For a general anisotropic material the detailed stress-strain
equations are (9.37), which reduce to (9.38) for the isotropic case.
10.4.2. Numerically Integrated Stiffness Matrix
Introducing e = Bu and = Ee into the TPE functional restricted to the element volume and
rendering the resulting algebraic form stationary with respect to the node displacements ue we get
the usual expression for the element stiffness matrix

K =
e

BT E B d
e .

(10.30)

Unlike the linear tetrahedron, the integrand varies over the element, and numerical integration by
Gauss rules is recommended. To apply these rules it is necessary to express (10.30) in terms of
tetrahedral coordinates. Since B is already a function of the i while E is taken to be constant, the
only adjustment needed is onthe volume differential d
e and the integration limits. Using (10.23)
we get

K =

BT E B
0

1015

1
6

J d1 d2 d3 d4 .

(10.31)

Chapter 10: THE QUADRATIC TETRAHEDRON

Denote the 30 30 matrix integrand by F(1 , 2 , 3 , 4 ) = BT EB 16 J . Application of a p-point

Gauss integration rule with abcissas {1k 2k 3k 4k } and weights wk with k = 1, 2, . . . p, gives
Ke =

p


wk F(1k 2k 3k 4k ).

(10.32)

k=1

The target rank of Ke is 30 6 = 24. Since each Gauss point adds 6 to the rank up to a maximum
of 24, the number of Gauss points should be 4 or higher. There is actually a tetrahedron Gauss rule
with 4 points and degree 2,5 which is the default one. For additional details on Gauss integration
rules over tetrahedra, see 10.7.
10.5. Element Stiffness Matrix
An implementation of the quadratic tetrahedron stiffness matrix computations as a Mathematica
module is listed in Figure 10.6. The module is invoked as
Ke=IsoTet10Stiffness[xyztet,Emat,{ },options];

(10.33)

The arguments are

xyztet

Element node coordinates, arranged as a two-dimensional list:

{ { x1,y1,z1 },{ x2,y2,z2 },{ x3,y3,z3 },...{ x10,y10,z10 } }.

Emat

The 6 6 matrix of elastic moduli (?) provided as a two-dimensional list:

{ { E11,E12,E13,E14,E15,E16 }, ... { E16,E26,E36,E46,E56,E66 } }.

options

A list containing optional additional information. For this element:

options={ numer,p,e }, in which
numer

A logical flag: True to specify floating-point numeric work, False to

request exact calculations. If omitted, False is assumed.

Gauss integration rule identifier. See 10.7 for instructions on how to

designate rules. If omitted, p=4 is assumed.

Element number. Only used in error messages. If omitted, 0 is assumed.

The third argument is a placeholder and should be set to the empty list { }.
The stiffness module calls IsoTet10ShapeFunDer, which was described in 10.3.3, and listed in
Figure 10.5, to get the Cartesian partial derivatives of the shape functions as well as the Jacobian
determinant. Note that in the implementation of Figure 10.6, Be is J B and not B. The correct
scaling is restored in the Ke+=(wk/(6*Jdet))*Transpose[Be].(Emat.Be) statement. (For the
provenance of the / factor, see (10.31).)
Module IsoTet10Stiffness is exercised by the Mathematica script listed in Figure 10.7, which
forms the stiffness matrix of a linear tetrahedron with corner coordinates
xyztet={ { 2,3,4 },{ 6,3,2 },{ 2,5,1 },{ 4,3,6 } }.
5

(10.34)

Meaning that it integrates exactly polynomials of degree 2 or higher in tetrahedral coordinates if the tetrahedron is CM.

1016

10.6 *THE JACOBIAN MATRIX PARADOX

IsoTet10Stiffness[xyztet_,Emat_,{},options_]:= Module[{e=0,
k,ip,p=4,tcoork,wk,numer=False,Nfx,Nfy,Nfz,
Jdet,Be,Ke}, Ke=Table[0,{30},{30}];
If [Length[options]>=1, numer=options[[1]]];
If [Length[options]>=2, p=options[[2]]];
If [Length[options]>=3, e=options[[3]]];
For [k=1,k<=Abs[p],k++,
{tcoork,wk}=TetrGaussRuleInfo[{p,numer},k];
{Nfx,Nfy,Nfz,Jdet}=IsoTet10ShapeFunDer[xyztet,tcoork,numer];
If [numer&&(Jdet<=0), Print["IsoTet10Stiffness: Neg "
"or zero Jacobian, element," e]; Return[Null]];
Be= {Flatten[Table[{Nfx[[i]],0,0},{i,1,10}]],
Flatten[Table[{0,Nfy[[i]],0},{i,1,10}]],
Flatten[Table[{0,0,Nfz[[i]]},{i,1,10}]],
Flatten[Table[{Nfy[[i]],Nfx[[i]],0},{i,1,10}]],
Flatten[Table[{0,Nfz[[i]],Nfy[[i]]},{i,1,10}]],
Flatten[Table[{Nfz[[i]],0,Nfx[[i]]},{i,1,10}]]};
Ke+=(wk/(6*Jdet))*Transpose[Be].(Emat.Be)];
If [!numer,Ke=Simplify[Ke]]; Return[Ke] ];

Figure 10.6. Quadratic tetrahedron stiffness matrix module.

ClearAll[Em, ,p,numer]; Em=480; =1/3; p=4; numer=False;

Emat=Em/((1+)*(1-2*))*{{1-,,,0,0,0},
{,1-,,0,0,0},{,,1-,0,0,0},{0,0,0,1/2-,0,0},
{0,0,0,0,1/2-,0},{0,0,0,0,0, 1/2-}};
{{x1,y1,z1},{x2,y2,z2},{x3,y3,z3},{x4,y4,z4}}=
{{2,3,4},{6,3,2},{2,5,1},{4,3,6}};
{{x5,y5,z5},{x6,y6,z6},{x7,y7,z7},{x8,y8,z8},{x9,y9,z9},{x10,y10,z10}}=
{{x1+x2,y1+y2,z1+z2},{x2+x3,y2+y3,z2+z3},{x3+x1,y3+y1,z3+z1},
{x1+x4,y1+y4,z1+z4},{x2+x4,y2+y4,z2+z4},{x3+x4,y3+y4,z3+z4}}/2;
xyztet={{x1,y1,z1},{x2,y2,z2},{x3,y3,z3},{x4,y4,z4},{x5,y5,z5},
{x6,y6,z6},{x7,y7,z7},{x8,y8,z8},{x9,y9,z9},{x10,y10,z10}};
Ke=IsoTet10Stiffness[xyztet,Emat,{},{numer,p}];
KeL=Table[Table[Ke[[i,j]],{j, 1,15}],{i,1,30}];
KeR=Table[Table[Ke[[i,j]],{j,16,30}],{i,1,30}];
Print["Ke=",KeL//MatrixForm]; Print[KeR//MatrixForm];
Print["eigs of Ke=",Chop[Eigenvalues[N[Ke]]]];

Figure 10.7. Test script for stiffness matrix module.

Midpoint node coordinates are set by averaging coordinates of adjacent corner nodes; consequently
this test element has constant metric. Its volume is +24. The material is isotropic with E = 480 and
= 1/3. The results are shown in Figure 10.8. The computation of stiffness matrix eigenvalues
is always a good programming test, since six eigenvalues (associated with rigid body modes) must
be exactly zero while the other six must be real and positive. This is verified by the results shown
at the bottom of Figure 10.8.
10.6.

*The Jacobian Matrix Paradox

There are some dark alleys in the Jacobian matrix computations covered in (10.3). Those are elucidated
here in some detail, since they are not discussed in the FEM literature. They are manifested when the general
tetrahedral geometry is specialized to the constant metric (CM) case. If so the Jacobian matrix J should be
constant and reduce to that of the linear tetrahedron. Which will not necessarily happen. This mismatch may
cause grief to element implementors, and prompt a fruitless search for bugs. But it is not a bug: only a side

1017

Chapter 10: THE QUADRATIC TETRAHEDRON

Ke=

447
324
72
1
6
12
54
48
0
94
66
36
152
90
12
55
42
12
311
252
24
431
306
132
95
60
24
148
114
36
55
48
0
83
54
12
90
72
0
26
30
12
392
312
48
376
0
96
152
192
0
116
72
48
292
144
0
136
288
96

324
1032
162
24
104
42
24
216
12
60
232
84
180
32
72
48
112
30
180
992
90
288
1040
306
84
128
42
84
448
96
42
112
6
90
260
90
36
360
36
24
68
12
336
1424
96
0
928
0
96
256
96
72
176
72
216
736
216
144
256
144

72
162
339
0
30
35
0
24
54
24
60
94
24
36
8
0
6
19
24
126
275
96
234
395
24
30
59
24
84
148
12
30
19
12
54
83
0
72
90
0
12
10
0
96
248
96
0
376
96
192
136
48
72
116
48
72
148
0
288
152

1
24
0
87
54
36
18
24
0
10
18
12
32
54
12
83
90
12
19
0
0
11
6
12
59
72
48
28
42
12
311
180
24
19
18
12
198
144
0
58
54
36
232
456
144
152
96
96
1048
576
192
292
168
48
308
144
96
680
672
288

6
104
30
54
132
54
12
72
12
0
76
36
36
268
72
54
260
54
18
32
18
6
28
6
18
272
126
12
148
48
252
992
126
0
32
18
72
792
36
36
88
36
336
352
96
192
256
192
576
2176
288
192
736
168
144
224
72
480
1568
528

12
42
35
36
54
87
0
24
18
0
36
46
48
108
76
12
90
83
12
18
17
12
6
11
12
126
167
0
60
64
24
90
275
0
18
17
0
72
198
24
36
58
96
192
232
0
96
136
192
288
760
0
120
148
96
72
164
192
480
680

54
24
0
18
12
0
108
0
0
36
12
0
72
12
0
90
36
0
198
72
0
18
12
0
18
24
0
72
36
0
431
288
96
11
6
12
18
24
0
350
234
132
136
216
48
116
72
48
292
192
0
984
648
144
152
72
48
392
408
48

48
216
24
24
72
24
0
432
0
24
144
48
24
288
48
72
360
72
144
792
72
24
72
24
48
72
24
72
288
144

0
12
54
0
12
18
0
0
108
0
24
36
0
24
72
0
36
90
0
36
198
0
12
18
0
12
18
0
72
72

306
1040
234
6
28
6
12
72
12
216
860
324
216
256
144
72
176
72
168
736
120
648
2208
432
216
256
144
240
1424
48

132
306
395
12
6
11
0
24
18
96
252
386
48
288
152
48
72
116
48
168
148
144
432
984
48
144
136
0
48
248

94
60
24
10
0
0
36
24
0
204
108
72
104
60
24
26
24
0
58
36
24
350
216
96
98
36
24
40
36
24
95
84
24
59
18
12
18
48
0
98
18
12
680
648
240
292
216
48
308
144
96
152
216
48
696
216
144
232
504
144

66
232
60
18
76
36
12
144
24
108
492
216
48
308
96
30
68
12
54
88
36
234
860
252
18
392
180
0
268
0

36
84
94
12
36
46
0
48
36
72
216
312
24
120
140
12
12
10
36
36
58
132
324
386
12
180
242
24
72
4

60
128
30
72
272
126
24
72
12
36
392
180
504
1568
432
144
736
72
144
224
72
72
256
144
216
1056
432
288
352
144

152
180
24
32
36
48
72
24
0
104
48
24
1416
648
144
392
336
0
232
336
96
136
216
48
680
504
240
704
288
96

24
42
59
48
126
167
0
24
18
24
180
242
240
576
680
0
216
148
96
72
164
48
144
136
144
432
696
96
144
232

90
32
36
54
268
108
12
288
24
60
308
120
648
3936
864
312
1424
96
456
352
192
216
256
288
648
1568
576
288
2384
576

148
84
24
28
12
0
72
72
0
40
0
24
704
288
96
136
144
0
680
480
192
392
240
0
232
288
96
1120
432
192

12
72
8
12
72
76
0
48
72
24
96
140
144
864
1416
48
96
248
144
96
232
48
144
152
240
432
680
96
576
848

114
448
84
42
148
60
36
288
72
36
268
72
288
2384
576
288
256
288
672
1568
480
408
1424
48
504
352
144
432
3616
864

36
96
148
12
48
64
0
144
72
24
0
4
96
576
848
96
144
152
288
528
680
48
48
248
144
144
232
192
864
1408

eigs of Ke = {88809.45, 49365.01, 2880.56, 2491.66, 2004.85, 1632.49, 1264.32, 1212.42, 817.907, 745.755, 651.034, 517.441,
255.100, 210.955, 195.832, 104.008, 72.7562, 64.4376, 53.8515, 23.8417, 16.6354, 9.54682, 6.93361, 2.22099, 0, 0, 0, 0, 0, 0}

Figure 10.8. Result from stiffness matrix test.

1018

10.6 *THE JACOBIAN MATRIX PARADOX

(a)

1 (1=1,2 =0) 3 (1=2 =1/2) 2 (1=0,2 =1)

x1

x3

x2 = x1+L

(b)
N1 = (21 1) 1
1

N3 = 4 1 2

N2 = (22 1) 2

2 1

2 1

Figure 10.9. Three-node bar element for case study in ?.

effect of the use of non-independent natural coordinates. In fact, the phenomenon may occur in all higher
order iso-P elements developed with that kind of coordinates. For instance, the six-node plane stress triangle.
Covering this topic for the quadratic tetrahedron leads to an algebraic maze. It is more readily illustrated with
the simplest quadratic iso-P element: the three-node bar introduced in Chapter 16 of IFEM. See Figure 10.9.
Instead of the single iso-P natural coordinate employed there, we will use two natural coordinates: 1 and 2 ,
defined as shown in Figure 10.9(a). Both coordinates vary from 0 to 1, and are linked by 1 + 2 = 1. Quick
preview of what may happen: F = 1 and G = 1 (1 + 2 ), say, are obviously the same function. But the
formal partial derivatives F/1 = 1 and G/1 = 21 + 2 = 1 + 1 are plainly different. The discrepancy
comes from the constraint 1 + 2 = 1 not being properly accounted for in the differentiation.
For the element of Figure 10.9(a), define for convenience
x0 =

x1 + x2
,
2

x3 =

x1 + x2
+ (x2 x1 ) = x0 + L .
2

(10.35)

The last relation defines the position of node 3 hierarchically, in terms of the midpoint deviation x3 = L.
Here / < < /. The element is of constant metric (CM) if = 0. The shape functions, depicted in
Figure 10.9(b), are N1 = (21 1)1 , N2 = (22 1)2 , and N3 = 41 2 . Two iso-P definitions for the axial
coordinate x are examined. The conventional (C) definition is
xC = x1 N1 + x2 N2 + x3 N3 = x1 (21 1)1 + x2 (22 1)2 + 4x3 1 2 ,

(10.36)

The hierarchical (H) definition is

x H = x1 (N1 + 12 N3 ) + x2 (N2 + 12 N3 ) + (x2 x1 ) N3 = x1 1 + x2 2 + 4 L 1 2 .

(10.37)

The difference is
xC x H = C H (1 + 2 1),

C H = 2(x1 1 + x2 2 ).

(10.38)

in which C H may be viewed as a Lagrange multiplier operating on the constraint 1 + 2 = 1.

The Jacobian entries associated with (10.36) are
JC x1 =

xC
= x1 (41 1) + 4x3 2 ,
1

JC x2 =

xC
= x2 (42 1) + 4x3 1 .
2

(10.39)

On introducing x3 = (x1 + x2 )/2 + L, and replacing 1 1 and 1 2 as appropriate, the Jacobian matrix
emerges as
JC =

JC x1

JC x2

1
x1 (3 22 ) + 2(x2 + L)2

1019

1
2(x1 + L)1 + x2 (3 21 )

(10.40)

Chapter 10: THE QUADRATIC TETRAHEDRON

Setting = 0, the resulting JC differs substantially from the expected constant-metric result:
JC |0 =

1
x1 (3 + 22 ) + 2x2 2

1
1
= JC M =
2x1 1 + x2 (3 + 21 )
x1

1
.
x2

(10.41)

On the other hand, going through the same process with x H as defined in (10.37) yields
JH =

J H x1

J H x2

1
x1 + 4 L 2

1
.
x2 + 4 L 1

(10.42)

On setting = 0, J H reduces to JC M , as expected. Lesson: application of metric constraints on shape functions

and differentiation do not commute when natural coordinates are not independent. However, the Jacobian
determinants do agree:

JC = det(JC ) = L 1 + (1 2 ) ,

J H = det(J H ) = JC ,

(10.43)

so we may simply denote J = JC = J H . Taking inverses gives

JC1 =

1
2(x1 + 2L)1 + x2 (3 21 )
J x1 (3 22 ) + 2(x2 + 2L)2

1
,
1

J1
H =

1
x2 + 4 L 1
J x1 4 L2

1
1

(10.44)

The crucial portion of J1 is the last column, which stores 1 / x = 1/J and 2 / x = 1/J . As can
be seen, this portion is the same for both definitions. This coincidence is not accidental. Jacobian matrices
produced by different iso-P geometric definitions in any number of dimensions obey a common configuration:
1D: J =

1
1
,
Jx1 Jx Jx2 Jx


2D: J =

1
1
1
Jx1 Jx Jx2 Jx Jx3 Jx
Jy1 Jy Jy2 Jy Jy3 Jy


,

. . . etc.

(10.45)

where the J s come from differentiating Lagrange multiplier terms such as that shown in (10.38). These
terms do not change the shape functions; only their derivatives. It is readily verified through matrix algebra
that in k space dimensions, the Jacobian determinant J and the last k columns of the inverse Jacobian matrix
J1 do not depend on the J s. Hence all definitions coalesce. Conclusions for element work:
1.

If only J and/or the i -partials portion of J1 are needed, the iso-P geometric definition: conventional
or hierarchical (or any combination thereof) makes no difference.

2.

If J is needed, and getting the expected JC M for a CM element is important, use of the hierarchical
definition is recommended. However, the direct use of J (as opposed to its determinant) in element
development is rare.

10.7.

*Gauss Integration Rules For Tetrahedra

For convenience we restate here two key properties that a Gauss integration rule must have to be useful in
finite element work:
1.

Full Symmetry. The same result must be obtained if the local element numbers are cyclically renumbered,
which modifies the natural coordinates.6

2.

Positivity. All integration weights must be positive.

As in the case of the triangle discussed in IFEM, fully symmetric integration rules over tetrahedra must be of
non-product type.
6

Stated mathematically: the numerically evaluated integral must remain invariant under all affine transformations of the
element domain onto itself.

1020

10.7 *GAUSS INTEGRATION RULES FOR TETRAHEDRA

Table 10.1. Tetrahedron Integration Formulas
Ident

Stars

Points Degree Comments

S4

Centroid formula, useful for Tet4 stiffness

S31

Useful for Tet4 mass and Tet10 stiffness

2S31

2S31

Has corners and face centers as sample points

14

2S31 + S22

14

Useful for Tet10 mass; exact form unavailable

14

2S31 + S22

14

Has edge midpoints as sample points

15

S4 + 2S31 + S22

15

Useful for Tet21 stiffness

15

S4 + 2S31 + S22

15

Less accurate than above one

3S21 + S211

24

Useful for Tet21 mass; exact form unavailable

24

The 5-point, degree-3, symmetric rule listed in some FEM textbooks has one negative weight.

To illustrate the first requirement, suppose that the i th sample point has tetrahedral coordinates {1i , 2i , 3i 4i }
linked by 1i + 2i + 3i + 4i = 1. Then all points obtained by permuting the 4 indices must be also sample
points and have the same weight wi . If the four values are different, permutations provide 24 sample points:
{1i , 2i , 3i , 4i }, {1i , 2i , 4i , 3i }, {1i , 3i , 2i , 4i }, {1i , 3i , 4i , 2i }, . . . {4i , 3i , 2i , 1i }.

(10.46)

This set of points is said to form a sample point star or simple star, which is denoted by S1111 . If two values
are equal, the set (10.46) coalesces to 12 different points, and the star is denoted by S211 . If two point pairs
coalesce, the set (10.46) reduces to 6 points, and the star is denoted by S22 . If three values coalesce, the set
(10.46) reduces to 4 points and the star is denoted by S31 . Finally if the four values coalesce, which can only
happen for 1i = 2i = 3i = 4i = 14 , the set (10.46) coalesces to one point and the star is denoted by S4
Sample point stars S4 , S31 , S22 , S211 and S1111 have 1, 4, 6, 12 or 24 points, respectively. Thus possible rules
can have i + 4 j + 6k + 12l + 24m points, where i, j, k, l, m are nonegative integers and i is 0 or 1. This
restriction exclude rules with 2, 3, 8 and 11 points. Table 10.1 lists nine FEM-useful rules for the tetrahedral
geometry, all of which comply with the requirements listed above. Several of these are well known while
others were derived by the author.
The nine rules listed in Table 10.1 are implemented in a self-contained Mathematica module called
TetrGaussRuleInfo. Because of its length the module logic is split in two Figures: 10.10 and 10.11.
It is invoked as
{ { zeta1,zeta2,zeta3,zeta4 },w } = TetrGaussRuleInfo[{ rule,numer },i]

(10.47)

The module has three arguments: rule, numer and i. The first two are grouped in a two-item list.
Argument rule, which can be 1, 4, 8, 8, 14, 14, 15, 15 or 24, identifies the integration formula as follows.
Abs[rule] is the number of sample points. If there are two useful rules with the same number of points, the
most accurate one is identified with a positive value and the other one with a minus value. For example, there
are two useful 8-point rules. If rule=8 a formula with all interior points is chosen. If rule=-8 a formula
with sample points at the 4 corners and the 4 face centers (less accurate but simpler and easy to remember) is
picked.
Logical flag numer is set to True or False to request floating-point or exact information, respectively, for
rules other than +14 or +24. For the latter see below.

1021

Chapter 10: THE QUADRATIC TETRAHEDRON

TetrGaussRuleInfo[{rule_,numer_},point_]:= Module[{
jk6= {{1,2},{1,3},{1,4},{2,3},{2,4},{3,4}},
jk12={{1,2},{1,3},{1,4},{2,3},{2,4},{3,4},
{2,1},{3,1},{4,1},{3,2},{4,2},{4,3}},
i=point,j,k,g1,g2,g3,g4,h1,w1,w2,w3,eps=10.^(-16),
info={{Null,Null,Null,Null},0} },
If [rule==1, info={{1/4,1/4,1/4,1/4},1}];
If [rule==4, g1=(5-Sqrt[5])/20; h1=(5+3*Sqrt[5])/20;
info={{g1,g1,g1,g1},1/4}; info[[1,i]]=h1];
If [rule==8, j=i-4;
g1=(55-3*Sqrt[17]+Sqrt[1022-134*Sqrt[17]])/196;
g2=(55-3*Sqrt[17]-Sqrt[1022-134*Sqrt[17]])/196;
w1=1/8+Sqrt[(1715161837-406006699*Sqrt[17])/23101]/3120;
w2=1/8-Sqrt[(1715161837-406006699*Sqrt[17])/23101]/3120;
If [j<=0,info={{g1,g1,g1,g1},w1}; info[[1,i]]=1-3*g1];
If [j> 0,info={{g2,g2,g2,g2},w2}; info[[1,j]]=1-3*g2]];
If [rule==-8, j=i-4;
If [j<=0,info={{0,0,0,0}, 1/40}; info[[1,i]]=1];
If [j> 0,info={{1,1,1,1}/3,9/40}; info[[1,j]]=0]];
If [rule==14,
(* g1,g2 +roots of P(g)=0, P=9+96*g1712*g^2-30464*g^3-127232*g^4+86016*g^5+1060864*g^6 *)
g1=0.09273525031089122640232391373703060;
g2=0.31088591926330060979734573376345783;
g3=0.45449629587435035050811947372066056;
If [!numer,{g1,g2,g3}=Rationalize[{g1,g2,g3},eps]];
w1=(-1+6*g2*(2+g2*(-7+8*g2))+14*g3-60*g2*(3+4*g2*
(-3+4*g2))*g3+4*(-7+30*g2*(3+4*g2*(-3+4*g2)))*g3^2)/
(120*(g1-g2)*(g2*(-3+8*g2)+6*g3+8*g2*(-3+4*g2)*g3-4*
(3+4*g2*(-3+4*g2))*g3^2+8*g1^2*(1+12*g2*
(-1+2*g2)+4*g3-8*g3^2)+g1*(-3-96*g2^2+24*g3*(-1+2*g3)+
g2*(44+32*(1-2*g3)*g3))));
w2=(-1-20*(1+12*g1*(2*g1-1))*w1+20*g3*(2*g3-1)*(4*w1-1))/
(20*(1+12*g2*(2*g2-1)+4*g3-8*g3^2));
If [i<5,
info={{g1,g1,g1,g1},w1};info[[1,i]]=1-3*g1];
If [i>4&&i<9, info={{g2,g2,g2,g2},w2};info[[1,i-4]]=1-3*g2];
If [i>8,
info={{g3,g3,g3,g3},1/6-2*(w1+w2)/3};
{j,k}=jk6[[i-8]]; info[[1,j]]=info[[1,k]]=1/2-g3] ];
If [rule==-14,
g1=(243-51*Sqrt[11]+2*Sqrt[16486-9723*Sqrt[11]/2])/356;
g2=(243-51*Sqrt[11]-2*Sqrt[16486-9723*Sqrt[11]/2])/356;
w1=31/280+Sqrt[(13686301-3809646*Sqrt[11])/5965]/600;
w2=31/280-Sqrt[(13686301-3809646*Sqrt[11])/5965]/600;
If [i<5,
info={{g1,g1,g1,g1},w1};info[[1,i]]=1-3*g1];
If [i>4&&i<9, info={{g2,g2,g2,g2},w2};info[[1,i-4]]=1-3*g2];
If [i>8&&i<15,info={{0,0,0,0},2/105};
{j,k}=jk6[[i-8]]; info[[1,j]]=info[[1,k]]=1/2]];

Figure 10.10. Tetrahedral integration rule information module, Part 1 of 2.

Argument i is the index of the sample point, which may range from 1 through Abs[rule].
The module returns the list {{1 , 2 , 3 , 4 }, w}, where 1 , 2 , 3 , 4 are the natural coordinates of the sample point, and w is the integration weight. For example, TetrGaussRuleInfo[{ 4,False },2] returns
{ { (5-Sqrt[5])/20,(5+3*Sqrt[5])/20,(5-Sqrt[5])/20,(5-Sqrt[5])/20 },1/4 }.
If rule is not implemented the module returns { { Null,Null,Null,Null },0 }.
Exact information for rules +14 and +24 is either unknown or only partly known. For these the abscissas are
given in floating-point form with 36 exact digits. If flag numer is False, the abscissas are converted to rational
numbers that represent the data to 16 place accuracy, using the built-in function Rationalize. Weights are
recovered from the abscissas. The conversion precision can be changed by adjusting the value of internal
variable eps, which is set in the module preamble.
We specifically mention the two simplest numerical integration rules, which find applications in the evaluation

1022

10.7 *GAUSS INTEGRATION RULES FOR TETRAHEDRA

If [rule==15,
g1=(7-Sqrt[15])/34; g2=7/17-g1; g3=(10-2*Sqrt[15])/40;
w1=(2665+14*Sqrt[15])/37800; w2=(2665-14*Sqrt[15])/37800;
If [i<5,
info={{g1,g1,g1,g1},w1};info[[1,i]]=1-3*g1];
If [i>4&&i<9, info={{g2,g2,g2,g2},w2};info[[1,i-4]]=1-3*g2];
If [i>8&&i<15,info={{g3,g3,g3,g3},10/189};
{j,k}=jk6[[i-8]]; info[[1,j]]=info[[1,k]]=1/2-g3];
If [i==15,info={{1/4,1/4,1/4,1/4},16/135}] ];
If [rule==-15,
g1=(13-Sqrt[91])/52;
If [i<5,
info={{1,1,1,1}/3,81/2240};info[[1,i]]=0];
If [i>4&&i<9, info={{1,1,1,1}/11,161051/2304960};
info[[1,i-4]]=8/11];
If [i>8&&i<15,info={{g1,g1,g1,g1},338/5145};
{j,k}=jk6[[i-8]]; info[[1,j]]=info[[1,k]]=1/2-g1];
If [i==15,info={{1/4,1/4,1/4,1/4},6544/36015}] ];
If [rule==24,
g1=0.214602871259152029288839219386284991;
g2=0.040673958534611353115579448956410059;
g3=0.322337890142275510343994470762492125;
If [!numer,{g1,g2,g3}=Rationalize[{g1,g2,g3},eps]];
w1= (85+2*g2*(-319+9*Sqrt[5]+624*g2)-638*g324*g2*(-229+472*g2)*g3+96*(13+118*g2*(-1+2*g2))*g3^2+
9*Sqrt[5]*(-1+2*g3))/(13440*(g1-g2)*(g1-g3)*(3-8*g2+
8*g1*(-1+2*g2)-8*g3+16*(g1+g2)*g3));
w2= -(85+2*g1*(-319+9*Sqrt[5]+624*g1)-638*g324*g1*(-229+472*g1)*g3+96*(13+118*g1*(-1+2*g1))*g3^2+
9*Sqrt[5]*(-1+2*g3))/(13440*(g1-g2)*(g2-g3)*(3-8*g2+
8*g1*(-1+2*g2)-8*g3+16*(g1+g2)*g3));
w3= (85+2*g1*(-319+9*Sqrt[5]+624*g1)-638*g224*g1*(-229+472*g1)*g2+96*(13+118*g1*(-1+2*g1))*g2^2+
9*Sqrt[5]*(-1+2*g2))/(13440*(g1-g3)*(g2-g3)*(3-8*g2+
8*g1*(-1+2*g2)-8*g3+16*(g1+g2)*g3));
g4=(3-Sqrt[5])/12; h4=(5+Sqrt[5])/12; p4=(1+Sqrt[5])/12;
If [i<5,
info={{g1,g1,g1,g1},w1};info[[1,i]]= 1-3*g1];
If [i>4&&i<9, info={{g2,g2,g2,g2},w2};info[[1,i-4]]=1-3*g2];
If [i>8&&i<13,info={{g3,g3,g3,g3},w3};info[[1,i-8]]=1-3*g3];
If [i>12,info={{g4,g4,g4,g4},27/560};
{j,k}=jk12[[i-12]];info[[1,j]]=h4;info[[1,k]]=p4] ];
If [numer, Return[N[info]], Return[Simplify[info]]];
];

Figure 10.11. Tetrahedral integration rule information module, Part 2 of 2.

of the element stiffness matrix and consistent force vector for the quadratic tetrahedron.
One point rule (exact for constant and linear polynomials over CM tetrahedra):
1
V

F(1 , 2 , 3 , 4 ) d
e F( 14 , 14 , 14 , 14 ).

(10.48)

Four-point rule (exact for constant through quadratic polynomials over CM tetrahedra):
1
V

F(1 , 2 , 3 , 4 ) d
e 14 F(, , , )+ 14 F(, , , )+ 14 F(, , , )+ 14 F(, , , ), (10.49)

in which = (5 + 3 5)/20 = 0.58541020, = (5 5)/20 = 0.13819660.

(More material to be added, especially timing results)

1023

Chapter 10: THE QUADRATIC TETRAHEDRON

Notes and Bibliography

The 10-node quadratic tetrahedron is the natural 3D generalization of the 6-node quadratic triangle. Once the
latter was formulated by Fraeijs de Veubeke [278] the road to high order elements was open. The constant
metric (CM) version was first published by Argyris in 1965 [27], in which natural coordinates and exact
integration were used. Extensions to variable metric (VM) had to wait until the isoparametric formulation was
developed by Bruce Irons within Zienkiewicz group at Swansea. Once that formulation become generally
known by 1966 [398], VM tetrahedra of arbitrary order could be fitted naturally into one element family. The
first comprehensive description of the 2D and 3D iso-P families may be found in [830].
Examination of the publications of this period brings up some interesting quirks, which occasionally reveal
more about non invented here territorial hubris than actual contributions to FEM.
The calculation of partial derivatives of functions expressed in tetrahedral coordinates, covered in 10.3, is
not available in the literature. For the plane stress 6-node triangle it was formulated by the author during a
1975 software project at Lockheed Palo Alto Research Labs, which was published in [224], as well as in the
IFEM notes [257, Ch 24]. The extension from triangles to tetrahedra is immediate. As discussed in 10.6, the
Jacobian matrix that emerges from such calculations is not unique, but if correctly programmed all variants
lead to the same element equations.
The simplest Gauss integration rules for tetrahedra may be found in the monograph by Stroud and Secrest
[713], along with extensive references to original sources. Rules for 14 or more points cited in Table 10.1 and
implemented in the Mathematica module of Figures 10.10 and 10.11 are taken from the compendium [244].

1024

Exercises

Homework Exercises for Chapter 10

The Quadratic Tetrahedron

EXERCISE 10.1 [A:15] The iso-P 10-node tetrahedron element (Tet10) is converted into an 11-point

tetrahedron (abbreviation: Tet11) by adding an interior node point 11 located at the element center
e
. (You do not need to
1 = 2 = 3 = 4 = 1/4. Construct the shape functions N1e , N5e and and N11
write down the full element definition).
Hint: add a correction bubble by trying N1e = N 1e + c1 1 2 3 4 , in which N 1e is the 10-node tetrahedron
e
, think a while.
shape function given in (10.2); you just need to find c1 . Likewise for N5e . For N11
Interesting tidbit: why does the quartic N 1e + c1 1 2 3 4 satisfies interelement compatibility but the obvious
cubic N1e = c1 N 1e (1 14 ) does not? (Think of the polynomial variation along edges.)
EXERCISE 10.2 [A:20] Obtain the 11 shape functions of the Tet11 tetrahedron of the previous Exercise,

and verify that its sum is identically one (use the constraint 1 + 2 + 3 + 4 = 1 when proving that).
EXERCISE 10.3 [A:15] The next full-polynomial, iso-P member of the tetrahedron family is the cubic
tetrahedron (abbreviation: Tet20), which has 20 node points and 60 degrees of freedom. Where do you think
the nodes are located?

Hint: the last four: numbered 17, 18, 19, and 20, are located at the face centers; those are called face nodes.
EXERCISE 10.4 [A:20] Derive the corner-node shape function N1 and the face-node shape function N17 for
the 20-node cubic tetrahedron. Note: node 17 lies at the centroid of face 1-2-3. Assume the element has
constant metric (CM) to facilitate visualization.

Hint for N1 : cross all 20 nodes except 1 with three planes; those planes have equations 1 = C where the
constant C varies for each plane.
e
: find 3 planes that cross all 20 nodes except 17.
Hint for N17

EXERCISE 10.5 [A/C:25] Derive the 20 shape functions for the 20-node cubic tetrahedron, and verify that

their sum is exactly one. (Only possible in reasonable time with a computer algebra system).
EXERCISE 10.6 [A/C:20] Compute fe for a constant-metric (CM) 10-node tetrahedron if the body forces bx

and b y vanish, while bz = g (g is the acceleration of gravity) is constant over the element. Use the 4-point
Gauss rule (10.49) to evaluate the integral that gives the consistent node force expression


fe =

NT b d
e ,

in which

b=

0
0
g


,

(E10.1)

using the exact expressions for = (5 + 3 5)/20 and = (5 5)/20 in the 4-point rule, and V V e .
You need to list only the z force components. Check that the sum of the z force components is g V e . The
result for the corner nodes is physically unexpected.

1025