Finite Element Analysis using Condensation

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{align} \begin{bmatrix} K_{ff} & K_{fb} \\ K_{bf} & K_{bb} \end{bmatrix} \begin{bmatrix} U_{f} \\ U_{b} \end{bmatrix} = \begin{bmatrix} f_{f} \\ f_{b} \end{bmatrix} \end{align}

It partitions the stiffness matrix K into $K_{ff}$ and internal $K_{bb}$ components. There are 2 ways to reduce degree of freedom, directory substitutes Boundary Condition, and Schur Compliment.

1: Substitutes Boundary Condition (Enforce)

Since we know the boundary conditions, we know what value that of displacement take $U_b = G_b$ .

In the most of the cases, $K_{bb}$ is set to unit matrix or quite big value and then, $f_b$ is set to $G_b$ so that $u_b = G_b$ $u_b = 0$ after the equation is solved. This means that it directly inputs BCs while condensing the dimension. It is computationally efficient, but it is not able to compute react force.

2: Schur Compliment (Condense)

Another way is to utilize Schur Compliment. As we learned in the Algebra class, it can be decomposed to inverse matrix with block matrices. So, if we know a part of matrix, we can use efficient computation. Let’s think of general case of Schur Compliment (theorem),

\mathbf{M} = \begin{bmatrix} \mathbf{A} & \mathbf{B} \\ \mathbf{C} & \mathbf{D} \end{bmatrix}

If M is not singular matrix, there is overall inverse matrix M is gotten by

\begin{align} \mathbf{M}^{-1} = \begin{bmatrix} \mathbf{A} - \mathbf{B} \mathbf{D}^{-1} \mathbf{C} & \mathbf{B} \mathbf{D}^{-1} \\ \mathbf{D}^{-1} \mathbf{C} & \mathbf{D}^{-1} \end{bmatrix}^{-1}\\ \mathbf{S} = \mathbf{A} - \mathbf{B} \mathbf{D}^{-1} \mathbf{C} \end{align}

this $\mathbf{S}$ is so called Schur Compliment. By condensing block matrix D, we get an new system (2). Let’s go back to Force Balance Equilibrium.

If $K_{dd}$ is invertible, we will obtaion condensed system by eliminating $u_b$ with Schur Compliment.

U_b = K_{bb}^{-1} (f_b - K_{bf} U_f)

(This is just brought from equation (1)). Substituting this into the degree of freedom equation in the (1) on the first row,

(K_{ff} - K_{fb} K_{bb}^{-1} K_{bf}) U_f = f_f - K_{fb} K_{bb}^{-1} f_b

This $(K_{ff} - K_{fb} K_{bb}^{-1} K_{bf})$ corresponds to Schur Compliment $S$ .　The computation cost is a bit higher than 1., but we can compute react force.

r_d = K_{bf} u_f - f_b

https://scikit-fem.readthedocs.io/en/latest/api.html#function-condense

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 = 10000  # 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(E, nu)

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]
    basis, [L/3, 2*L/3], [L/3, 2*W/3], [H, H]
).nodal['u^3']

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

K_e, F_e = skfem.enforce(K, F, D=fixed_dofs)
K_c, F_c, U_c, I = skfem.condense(K, F, D=fixed_dofs)
# K_e, F_e (2187, 2187) (2187,)
# K_c, F_c (1944, 1944) (1944,)


# Solve linear equation
# U_e = solve(K_e, F_e)
U_e = scipy.sparse.linalg.spsolve(K_e, F_e)
# U_c = solve(A_c, f_c)
U_c = scipy.sparse.linalg.spsolve(K_c, F_c)

U_c_full = np.zeros(K.shape[0])
free_dofs = np.setdiff1d(np.arange(K.shape[0]), fixed_dofs)
U_c_full[free_dofs] = U_c
reaction_force = K[fixed_dofs, :] @ U_c_full- F[fixed_dofs]
print("reaction_force:", np.sum(reaction_force))
# print(np.sum(U_e - U_c))

# ---------------------------
sf = 1.0
# print(U, np.max(U), np.min(U))
m = mesh.translated(sf * U_e[basis.nodal_dofs])
m.save('box_deformed.vtk')
mesh.save('box_org.vtk')