Simple Finite Element Analysis with Scikit-FEM

About

There is a very useful library that is able to solve Finite Element Analysis in Python. The name is “Scikit-FEM”. This is much more beautiful than I expected. I leave a simple test code here.

Structural Analysis with FEA

Theory – Static Condensation

The one of the most stable solution is “Static Condensation” to find force equilibrium on the equation. Static Condensation is a technique in the FEA that eliminates internal degrees of freedom (DOFs) to reduce the size of the system. Let’s say we have equation of force equilibrium, and we can decompose it as follows:

\begin{bmatrix}
   K_{bb} & K_{ic} \\
   K_{ib} & K_{ii}
   \end{bmatrix}
   \begin{bmatrix}
   u_{b} \\
   u_{i}
   \end{bmatrix}
   =
   \begin{bmatrix}
   f_{b} \\
   f_{i}
   \end{bmatrix}

It partitions the stiffness matrix K into boundary K_{bb} and internal K_{ii} components. And then, solve for the internal DOFs

u_i = K_{ii}^{-1} (f_{ii} - K_{ib} u_b)

And we can construct the reduced system

K_{\text{up}} = K_{bb} - K_{bi} K_{ii}^{-1}K_{ib}\\
f_{\text{up}} = f_b - K_{bi} K_{ii}^{-1} f_i)

Code

import numpy as np
import scipy.sparse.linalg
import scipy
import skfem
from skfem import *
from skfem.helpers import ddot, sym_grad, eye, trace
from skfem.models.elasticity import lame_parameters
from skfem.models.elasticity import linear_elasticity
from skfem.helpers import dot, sym_grad


def get_dofs_in_range(basis, x_range, y_range, z_range):
    return basis.get_dofs(
        lambda x: (x_range[0] <= x[0]) & (x[0] <= x_range[1]) &
                  (y_range[0] <= x[1]) & (x[1] <= y_range[1]) &
                  (z_range[0] <= x[2]) & (x[2] <= z_range[1])
    )

# Parameters
E = 210e9       # Young [Pa]
nu = 0.3        # Poisson
L, W, H = 1.0, 1.0, 1.0  # The size of Box
force_value = 100  # 100 Pa 
displacement_scale = 1  # Scale to Enhance the Displacement

mesh = skfem.MeshHex().refined(3).with_defaults()
e = skfem.ElementVector(skfem.ElementHex1())
basis = skfem.Basis(mesh, e, intorder=3)
lam, mu = lame_parameters(1e3, 0.3)

def C(T):
    return 2. * mu * T + lam * eye(trace(T), T.shape[0])

@skfem.BilinearForm
def stiffness(u, v, w):
    return ddot(C(sym_grad(u)), sym_grad(v))

K = stiffness.assemble(basis)


# Fixes dofs on plane (z=0) 
fixed_dofs = get_dofs_in_range(basis, [0, L], [0, W], [0, 0]).all()

# Get the dofs(z only) on the surface (z=H) 
top_dofs = get_dofs_in_range(basis, [0, L], [0, W], [H, H]).nodal['u^3']

# Distribute force
F = np.zeros(K.shape[0])
F[top_dofs] = force_value / len(top_dofs)

# Boundary Condition
K, F = skfem.enforce(K, F, D=fixed_dofs)

# Solve linear equation
# U = solve(K, F)
U = scipy.sparse.linalg.spsolve(K, F)
# ---------------------------
sf = 1.0
# print(U, np.max(U), np.min(U))
m = mesh.translated(sf * U[basis.nodal_dofs])
m.save('box_deformed.vtk')
mesh.save('box_org.vtk')

Result