Usage in scikit-fem

[adtzlr]

This is a little tutorial how to use matADi in sckit-fem. Let's set up a simple example - stretch of a clamped hyperelastic cube. We'll start with some necessary imports.

import numpy as np

import skfem as fem
from skfem.assembly import LinearForm, BilinearForm
from skfem.helpers import grad, identity

from matadi import MaterialHyperelastic
from matadi.models import neo_hooke

In the next step we define some variables, our material and create the linear and bilinear forms.

right = 0.5
tol = 1e-10

mat = MaterialHyperelastic(neo_hooke, C10=0.5, bulk=2.0)

m  = fem.MeshHex().refined(2)
e1 = fem.ElementHex1()
e  = fem.ElementVectorH1(e1)
ib = fem.Basis(m, e, intorder=1)

def P(w):
    H = grad(w["displacement"])
    F = H + identity(H)
    return mat.gradient([F])[0]

def A(w):
    H = grad(w["displacement"])
    F = H + identity(H)
    return mat.hessian([F])[0]

@LinearForm(nthreads=4)
def L(v, w):
    return np.einsum("ij...,ij...", P(w), grad(v))

@BilinearForm(nthreads=4)
def a(u, v, w):
    return np.einsum("ijkl...,ij...,kl...", A(w), grad(u), grad(v))

This will be followed by the definition of the displacement boundaries.

u  = ib.zeros()
du = ib.zeros()
uv = ib.interpolate(u)

dofs = {
    "left": ib.get_dofs(lambda x: x[0] == 0.0),
    "right": ib.get_dofs(lambda x: x[0] == 1.0),
}

I = ib.complement_dofs(dofs)

A newton-rhapson procedure is used to obtain equilibrium.

f = fem.asm(L, ib, displacement=uv)
K = fem.asm(a, ib, ib, displacement=uv)

for iteration in range(8):
    
    du_D = u.copy()
    du_D[dofs["right"].nodal['u^1']] = right - du_D[dofs["right"].nodal['u^1']]
    
    system = fem.condense(K, -f, x=du_D, I=I)
    du = fem.solve(*system)
    norm_du = np.linalg.norm(du)
    
    u += du
    
    uv = ib.interpolate(u)
    f = fem.asm(L, ib, displacement=uv)
    K = fem.asm(a, ib, ib, displacement=uv)
    
    print(1 + iteration, norm_du)
    
    if norm_du < tol:
        break

Finally the model is saved as VTK file.

savedata = {"u": u[ib.nodal_dofs].T}
m.save("result.vtk", savedata)

The ouput of this script shows the residuals of iterative displacements - convergence is fine.

0 3.4810410954530506
1 0.13063064781292275
2 0.009208944826176955
3 4.989804896048235e-05
4 2.3408069100383187e-09
5 3.2247002135208875e-16

[adtzlr]

Warning: The code could be outdated because this example was not tested for more than a year.