# Formulation of K for tetrahedral elements for the Finite Element Method

## A code made from the formulations of Jesep Birnic's Engineering ##Structures and Material Behavior and Engineering Computation of ##Structures: The Finite Element Method (Springer).

import meshio
import numpy as np
import symfem

# nominal_d in inches
nom_dim = 5
# E in PSI
E = 2000000
nu = 0.3

C = (E / ((1 + nu) * (1 - 2 * nu))) * np.array([
[1 - nu, nu, nu, 0, 0, 0],
[nu, 1 - nu, nu, 0, 0, 0],
[nu, nu, 1 - nu, 0, 0, 0],
[0, 0, 0, (1 - 2 * nu) / 2, 0, 0],
[0, 0, 0, 0, (1 - 2 * nu) / 2, 0],
[0, 0, 0, 0, 0, (1 - 2 * nu) / 2]
])

global points
points = [
[0.0, 0.0, 0.0],
[nom_dim, 0.0, 0.0],
[nom_dim, nom_dim, 0.0],
[0.0, nom_dim, 0.0],
[0.0, 0.0, nom_dim],
[nom_dim, 0.0, nom_dim],
[nom_dim, nom_dim, nom_dim],
[0.0, nom_dim, nom_dim],
]

global cells
cells = [
[0, 1, 3, 6],
[1, 3, 2, 6],
[7, 6, 3, 0],
[7, 6, 5, 1],
[7, 5, 4, 0]
]

def parray(i, cells, points):
pi1 = cells[i][0]
pi2 = cells[i][1]
pi3 = cells[i][2]
pi4 = cells[i][3]

return [
points[pi1], points[pi2], points[pi3], points[pi4]
]

lagrange = []
for i in range(len(cells)):
e = symfem.create_element("tetrahedron", "Lagrange", 1)
N = e.tabulate_basis(parray(0, cells, points))
lagrange.append([e,N])

def volume(pa):
x1, y1, z1 = pa[0]
x2, y2, z2 = pa[1]
x3, y3, z3 = pa[2]
x4, y4, z4 = pa[3]

V = np.array([
[1, x1, y1, z1],
[1, x2, y2, z2],
[1, x3, y3, z3],
[1, x4, y4, z4]
])
return abs(np.linalg.det(V))

def apply_dirichlet_bc(K, f, node, dof):
index = node * 3 + dof
K[index, :] = 0
K[:, index] = 0
K[index, index] = 1.0
f[index] = 0
return K, f

def get_cmats(j,k,l):
a_i = np.array([
[j[0], j[1], j[2]],
[k[0], k[1], k[2]],
[l[0], l[1], l[2]]
])

b_i = np.array([
[1, j[1], j[2]],
[1, k[1], k[2]],
[1, l[1], l[2]]
])

c_i = np.array([
[j[1], 1, j[2]],
[k[1], 1, k[2]],
[l[1], 1, l[2]],
])

d_i = np.array([
[j[1], j[2], 1],
[k[1], k[2], 1],
[l[1], l[2], 1]
])

return np.linalg.det(a_i), np.linalg.det(b_i), np.linalg.det(c_i), np.linalg.det(d_i)

def coefficients(pa):
cos = []
for i in range(0,4,1):
i = pa[0]
j = pa[1]
k = pa[2]
l = pa[3]

a_i,b_i,c_i,d_i = get_cmats(j,k,l)

pa[3] = i
pa[0] = j
pa[1] = k
pa[2] = l

cos.append([a_i,b_i,c_i,d_i])
return cos

def barycentric_coords(cos, pa):
CO = np.array([
[cos[0][0], cos[1][0], cos[2][0], cos[3][0]],
[cos[1][0], cos[1][1], cos[1][2], cos[1][3]],
[cos[2][0], cos[2][1], cos[2][2], cos[2][3]],
[cos[3][0], cos[3][1], cos[3][2], cos[3][3]]
])
V = volume(pa)
return ((1/6)*V)*CO

def co_r(i, cos):
return cos[i][0], cos[i][1], cos[i][2], cos[i][3]

def get_B(cos, pa):
a1,b1,c1,d1 = co_r(0,cos)
a2,b2,c2,d2 = co_r(1,cos)
a3,b3,c3,d3 = co_r(2,cos)
a4,b4,c4,d4 = co_r(3,cos)

V = volume(pa)

B = ((1/2)*V)*np.array([
[b1,0,0,b2,0,0,b3,0,0,b4,0,0],
[0,c1,0,0,c2,0,0,c3,0,0,c4,0],
[0,0,d1,0,0,d2,0,0,d3,0,0,d4],
[c1,b1,0,c2,b2,0,c3,b3,0,c4,b4,0],
[0,d1,c1,0,d2,c2,0,d3,c2,0,d4,c2],
[d1,0,b1,d2,0,b2,d3,0,b3,d4,0,b4]
])

return B

def compute_ke(pa, C):
cos = coefficients(pa)
#cos = barycentric_coords(cos, pa)
B = get_B(cos, pa)
V = volume(pa)
Ke = V * np.dot(np.dot(B.T, C), B)
return Ke

# Number of nodes
n_nodes = len(points)
# Number of degrees of freedom per node (3 for 3D problems)
dof_per_node = 3
# Total degrees of freedom
total_dof = n_nodes * dof_per_node

# Initialize global stiffness matrix
K_global = np.zeros((total_dof, total_dof))

# Assemble global stiffness matrix
for i, cell in enumerate(cells):
pa = parray(i, cells, points)
Ke = compute_ke(pa, C)

# Map local stiffness matrix to global matrix
for a in range(4):
for b in range(4):
for i in range(dof_per_node):
for j in range(dof_per_node):
K_global[cell[a] * dof_per_node + i, cell[b] * dof_per_node + j] += Ke[a * dof_per_node + i, b * dof_per_node + j]

# Force vector
f_global = np.zeros(total_dof)

# Force in pounds...
f_global[4] = -5.0

# Apply Dirichlet BC to node 0 (fix all DOF)
for dof in range(dof_per_node):
K_global, f_global = apply_dirichlet_bc(K_global, f_global, 0, dof)

# Solve for displacements
displacements = K_global @ f_global
print('displacements:',displacements)


The displacements look far too large. It seems something is wrong with this formulation so far. Why are the displacements so large for only a few pounds of force being input here?

displacements: [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  2.93438251e+15
-6.45564153e+15  2.93438251e+15  0.00000000e+00  0.00000000e+00
0.00000000e+00  0.00000000e+00  0.00000000e+00  1.17375300e+15
0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
1.17375300e+15  0.00000000e+00  0.00000000e+00  0.00000000e+00]


These are the mathematical formulas I used: LaTex's (mathematical formulation): \input $$$$\begin{bmatrix} 1 \\ x_1 \\ x_2 \\ x_3 \end{bmatrix} = \begin{bmatrix} 1 & 1 & 1 & 1 \\ x_1 & x_2^2 & x_3^3 & x_4^4 \\ x_2 & x_2^2 & x_3^3 & x_4^4 \\ x_3 & x_3^2 & x_3^3 & x_4^4 \end{bmatrix} \begin{bmatrix} \xi_{1} \\ \xi_{2} \\ \xi_{3} \\ \xi_{4} \end{bmatrix}$$$$

$$$$\begin{bmatrix} \xi_1 \\ \xi_2 \\ \xi_3 \\ \xi_4 \end{bmatrix} = \frac{1}{6V} \begin{bmatrix} a_1 & b_1 & c_1 & d_1 \\ a_2 & b_2 & c_2 & d_2 \\ a_3 & b_3 & c_3 & d_3 \\ a_4 & b_4 & c_4 & d_4 \end{bmatrix}$$$$

B = frac{1}{2V}$$\begin{vmatrix} b_{1} & 0 & 0 & b_{2} & 0 & 0 & b_{3} & 0 & 0 & b_{4} & 0 & 0 \\ 0 & c_{1} & 0 & 0 & c_{2} & 0 & 0 & c_{3} & 0 & 0 & c_{4} & 0 \\ 0 & 0 & d_{1} & 0 & 0 & d_{2} & 0 & 0 & d_{3} & 0 & 0 & d_{4} \\ c_{1} & b_{1} & 0 & c_{2} & b_{2} & 0 & c_{3} & b_{3} & 0 & c_{4} & b_{4} & 0 \\ 0 & d_{1} & c_{1} & 0 & d+{2} & c_{2} & 0 & d_{3} & c_{2} & 0 & d_{4} & c_{2} \\ d_{1} & 0 & b_{1} & d_{2} & 0 & b_{2} & d_{3} & 0 & b_{3} & d_{4} & 0 & b_{4} \end{vmatrix}$$

$$$$\begin{bmatrix} \xi_1 \\ \xi_2 \\ \xi_3 \\ \xi_4 \end{bmatrix} = \frac{1}{6V} \begin{bmatrix} a_1 & b_1 & c_1 & d_1 \\ a_2 & b_2 & c_2 & d_2 \\ a_3 & b_3 & c_3 & d_3 \\ a_4 & b_4 & c_4 & d_4 \end{bmatrix}$$$$

\begin{align*} x_1 &= \delta_{1}x_1^1 + \delta_{2}x_1^2 + \delta_{3}x_1^3 + \delta_{4}x_1^4 \\ x_2 &= \delta_{1}x_2^1 + \delta_{2}x_2^2 + \delta_{3}x_2^3 + \delta_{4}5x_2^4 \\ x_3 &= \delta_{1}x_3^1 + \delta+{2}x_3^2 + \delta_{3}x_3^3 + \delta_{4}x_3^4 \end{align*}


a_i = $$\begin{vmatrix} x_{1}^{j} & x_{2}^{j} & x_{3}^{j} \\ x_{1}^{k} & x_{2}^{k} & x_{3}^{k} \\ x_{1}^{l} & x_{2}^{l} & x_{3}^{l} \end{vmatrix}$$,
b_i = $$\begin{vmatrix} 1 & x_{2}^{j} & x_{3}^{j} \\ 1 & x_{2}^{k} & x_{3}^{k} \\ 1 & x_{2}^{l} & x_{3}^{l} \end{vmatrix}$$,
c_i = -$$\begin{vmatrix} x_{2}^{j} & 1 & x_{3}^{j} \\ x_{2}^{k} & 1 & x_{3}^{k} \\ x_{2}^{l} & 1 & x_{3}^{l} \end{vmatrix}$$,
d_i = -\begin{vmatrix}
x_{2}^{j} & x_{3}^{j}    & 1\\
x_{2}^{k} & x_{3}^{k}    & 1\\
x_{2}^{l} & x_{3}^{l}    & 1\\

\end{vmatrix}

• You will receive much better feedback if you give a more mathematical description of your problem, what you have implemented, and how it differs from the expected results. Commented May 17 at 17:32
• Adding an introduction and a comment to your code would make the question more understandable Commented May 17 at 17:37
• I added the mathematical formulas and titles that they were taken from along with an introductory comment of the code. Commented May 17 at 17:55