How to build K for a linear elasticity example of a tetrahedral mesh?

[prusso630]

import numpy as np
import meshio
import skfem

# Load points from CSV file
points = np.loadtxt('points.csv', delimiter=',')

# Load cells from CSV file
cells = np.loadtxt('cells.csv', delimiter=',', dtype=int)

# Create a mesh object
skfem_mesh = skfem.MeshTet(points.T, cells.T)

# Define the element type
element_type = skfem.ElementTetP1()

# Create basis
basis = skfem.InteriorBasis(skfem_mesh, element_type)

# 2. Define material properties
E = 2000000
nu = 0.3
mu = E / (2 * (1 + nu))
lmbda = E * nu / ((1 + nu) * (1 - 2 * nu))
C = np.array([[lmbda + 2 * mu, lmbda, lmbda, 0, 0, 0],
              [lmbda, lmbda + 2 * mu, lmbda, 0, 0, 0],
              [lmbda, lmbda, lmbda + 2 * mu, 0, 0, 0],
              [0, 0, 0, mu, 0, 0],
              [0, 0, 0, 0, mu, 0],
              [0, 0, 0, 0, 0, mu]])

# 3. Define boundary conditions
dofs = np.arange(points.shape[0] * 3)  # All degrees of freedom
dirichlet_dofs = np.array([
    0, 1, 2,  # Node 0: [ux, uy, uz]
    3, 4, 5,  # Node 1: [ux, uy, uz]
    6, 7, 8,  # Node 2: [ux, uy, uz]
    9, 10, 11  # Node 3: [ux, uy, uz]
])


f = np.zeros(basis.N)
# Force in pounds...
#f_global[23] = -5.0
#f_global[20] = -5.0
#f_global[17] = -5.0
#f_global[14] = -5.0

# 4. Assemble stiffness matrix and load vector
K = skfem.asm(form=basis.dim * [[None]], fun=skfem.asm.elasticity_stiffness, args=(C,), basis=basis)

# 5. Apply boundary conditions
for dof in dirichlet_dofs:
    K[dof, :] = 0
    K[:, dof] = 0
    K[dof, dof] = 1.0
    f[dof] = 0

# 6. Solve the system
u = np.linalg.solve(K, f)

# Compare results
print('Displacements:', u)

I am trying to get this example working to verify the results of an example I put together in the manual sense. I would like to see if my results are similar to what scikit-fem says. In this case points and cells are from a PyVista sphere that was tetrahedralized. So far my manual putting together says: U: [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 ... -2.67398613e-07
-5.13549340e-07 7.47276076e-08].

How can I get scikit-fem to start to reproduce these results. So far I was not yet able to build K from the provided information. How to get this sample attached above printing displacements U of the global displacements vector U?

[kinnala]

It seems that the code does not define any governing equation. In short, the idea in FEM is to turn a governing equation into a matrix equation. I assume you are interested in linear elasticity, then the governing equation consists of

$-\mathrm{div}\, \sigma(u)=f$

and

$\sigma(u) = 2 \mu \epsilon(u) + \lambda \mathrm{tr}\, \epsilon(u) I$

and

$\epsilon(u) = \frac{1}{2}(\nabla u + \nabla u^T).$

Your first step is to write these equations as code.

[kinnala]

In scikit-fem, you need to write the left hand side and the right hand side of the governing equation separately because the left hand side becomes a stiffness matrix and the right hand side becomes a load vector.

The left hand side becomes, after multiplying by the test function/variation and integration by parts, the following bilinear form:

$\int_\Omega \sigma(u) : \epsilon(v) \,dx$

This must be implemented using @BilinearForm decorator.

[kinnala]

This is something like:

from skfem import *
from skfem.helpers import *

@BilinearForm
def bilinf(u, v, _):
    epsu = 0.5 * (grad(u) + transpose(grad(u)))
    epsv = 0.5 * (grad(v) + transpose(grad(v)))
    sigu = 2 * mu * epsu + lmbda * eye(trace(epsu), u.shape[0])
    return ddot(sigu, epsv)