0
$\begingroup$
## 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{equation} \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} \end{equation}

\begin{equation} \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} \end{equation}

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{equation} \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} \end{equation}

\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}
$\endgroup$
3
  • 1
    $\begingroup$ 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. $\endgroup$
    – whpowell96
    Commented May 17 at 17:32
  • 1
    $\begingroup$ Adding an introduction and a comment to your code would make the question more understandable $\endgroup$
    – Rigel
    Commented May 17 at 17:37
  • $\begingroup$ I added the mathematical formulas and titles that they were taken from along with an introductory comment of the code. $\endgroup$
    – prusso
    Commented May 17 at 17:55

1 Answer 1

1
$\begingroup$

You have applied constraints at only a single node so your global stiffness matrix is singular. (Your linear equation solver should have reported this.) To prevent rigid body motion of your cube of material, you have to constrain at least six degrees of freedom total at three non-co-linear nodes. One way to achieve this is to constrain the three dofs at node 0 as you are now doing, constrain dofs 2 and 3 at node 1, and dof 3 at node 2.

If, in addition, you constrain dof 3 at node 3 and then apply a z-direction force at the last 4 nodes in the cube (4-7 with your zero-based numbering), you should get a uniform displacement state in the cube, You can easily compute this displacement state from your material properties, cube dimension, and applied force magnitude to compare.

$\endgroup$
5
  • $\begingroup$ Ke = compute_ke(pa,C), Ke[0, 0] += 1e-3, Ke[1, 1] += 1e-3, Ke[2, 2] += 1e-3, Ke += 1e-5 * np.eye(Ke.shape[0]), then if one applies dirichlet of all three degrees of freedom to nodes 0,1 and 2 then solve by u = Ke @ f the rsult is: u: [ 0.e+00 0.e+00 0.e+00 0.e+00 0.e+00 0.e+00 0.e+00 0.e+00 0.e+00 0.e+00 -5.e-05 0.e+00]. OK... seems fairly reasonable, still a bit difficult to resolve the singular matrix which is why the use of @. I would imagine from there a global K can be assembled. It's a type of system I still don't know much about how it works. Reasonable... maybe disaccurate. $\endgroup$
    – prusso
    Commented May 18 at 2:30
  • $\begingroup$ I could follow only a little of your last comment. I added another suggestion to my answer that might help you. $\endgroup$ Commented May 18 at 10:11
  • $\begingroup$ So a Z force of f_global index 23,20,17,14 to cover 4 nodes Z with scalar -5.0, The I used the existing apply_dirichlet_bc to fix node 0, dof 0-2, node 1, dof 1,2, node 2, dof 2 and node 3 dof 2 (I believe dofs are 0 indexed). As far as a vanilla solve the matrix is still singular. The @ method is producing erroneous results still for those dirichlets. There must be something wrong still. $\endgroup$
    – prusso
    Commented May 18 at 22:15
  • $\begingroup$ Calculate the eigenvalues of an element matrix. There should be six and only six essentially zero. $\endgroup$ Commented May 19 at 9:51
  • $\begingroup$ eigens Ke of i 0,[3.99623730e+13+0.j 2.21215286e+12+0.j 1.22876136e+13+0.j 5.49036383e-05+0.00118848j 5.49036383e-05-0.00118848j 1.12680288e+13+0.j 1.34224255e-03+0.j 0.00000000e+00+0.j 0.00000000e+00+0.j 0.00000000e+00+0.j 0.00000000e+00+0.j 0.00000000e+00+0.j ] $\endgroup$
    – prusso
    Commented May 19 at 17:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.