Implementing Ex.36 from scikit-fem in FElupe

[adtzlr]

Here are some possible implementations of Ex.36, taken from the docs of scikit-fem. These examples require pip install felupe[all] matadi to be installed.

import felupe as fem

mesh = fem.Cube(n=5).triangulate().add_midpoints_edges()
region = fem.RegionQuadraticTetra(mesh)
field = fem.FieldsMixed(region, n=2)

The mixed-field container automatically creates the Taylor-Hood element. The region of the first field is the one we created (RegionQuadraticTetra)

field[0].region

<felupe Region object>
  Element formulation: QuadraticTetra
  Quadrature rule: Tetrahedron
  Gradient evaluated: True

and the second field is based on a region one order lower, RegionTetra - but with the quadrature scheme from RegionQuadraticTetra.

field[1].region

<felupe Region object>
  Element formulation: Tetra
  Quadrature rule: Tetrahedron
  Gradient evaluated: False

Boundary conditions are created with a template for uniaxial tension. For more details on boundary conditions please see Boundary.

boundaries, loadcase = fem.dof.uniaxial(field, clamped=False)

In this first implementation, we'll use the gradients of the strain energy function, carried out by automatic differentiation using matadi.

import matadi as mat
from matadi.math import det, gradient, log, sqrt, trace

F = mat.Variable("F", 3, 3)
p = mat.Variable("p", 1, 1)
ζ = mat.Variable("ζ", 1, 1)


def ψ(x, μ, λ):
    F, p, ζ = x[0], x[1], x[2]

    J = det(F)
    Js = (λ + p + sqrt((λ + p) ** 2 + 4 * λ * μ)) / (2 * λ)

    C = F.T @ F
    I1 = trace(C)

    ψs = p * (J - Js) + μ / 2 * (I1 - 3) - μ * log(Js) + λ / 2 * (Js - 1) ** 2
    ζ_new = ζ

    return gradient(ψs, F), gradient(ψs, p), ζ_new


umat = mat.MaterialTensor(x=[F, p, ζ], fun=ψ, kwargs=dict(μ=1, λ=1e4), statevars=1)
solid = fem.SolidBody(umat, field)

move = fem.math.linsteps([0, 1], num=5)
step = fem.Step(items=[solid], ramp={boundaries["move"]: move}, boundaries=boundaries)

job = fem.CharacteristicCurve(steps=[step], boundary=boundaries["move"])
job.evaluate()
fig, ax = job.plot(
    xlabel="Displacement $u$ in mm $\longrightarrow$",
    ylabel="Normal Force $F$ in N $\longrightarrow$",
)

ax2 = field.imshow("Principal Values of Logarithmic Strain")

[adtzlr]

FElupe also has a very similar built-in ThreeFieldVariation (u,p,J). It basically slices out the volumetric part of your strain energy function, adds the constraints between (u,p,J) and eliminates volumetric locking. This time we'll use the Hyperelastic material definition from FElupe. Once again with automatic differentiation but with tensortrax as backend.

Here's the code:

import felupe as fem

mesh = fem.Cube(n=5).triangulate().add_midpoints_edges()
region = fem.RegionQuadraticTetra(mesh)
field = fem.FieldsMixed(region, n=3)

boundaries, loadcase = fem.dof.uniaxial(field, clamped=False)

import tensortrax.math as tm

def ψ(C, μ, λ):
    
    I1 = tm.trace(C)
    J = tm.sqrt(tm.linalg.det(C))
    lnJ = tm.log(J)
    
    return μ / 2 * (I1 - 3) - μ * lnJ + λ / 2 * (J - 1) ** 2

umat = fem.ThreeFieldVariation(fem.Hyperelastic(ψ, μ=1, λ=1e4, parallel=True))
solid = fem.SolidBody(umat, field)

move = fem.math.linsteps([0, 1], num=5)
step = fem.Step(items=[solid], ramp={boundaries["move"]: move}, boundaries=boundaries)

job = fem.CharacteristicCurve(steps=[step], boundary=boundaries["move"])
job.evaluate(tol=1e-2)
fig, ax = job.plot(
    xlabel="Displacement $u$ in mm $\longrightarrow$",
    ylabel="Normal Force $F$ in N $\longrightarrow$",
)

ax2 = field.imshow("Principal Values of Logarithmic Strain")

[adtzlr]

FYI, there is already a code-block for finite-strain-viscoelasticity, please see the docstring of Hyperelastic.

[adtzlr]

The strain energy function from Ex.36 is implemented in FElupe as NeoHookeCompressible. Instead of defining a custom strain energy function, one could also use this class instead.

import felupe as fem

mesh = fem.Cube(n=5).triangulate().add_midpoints_edges()
region = fem.RegionQuadraticTetra(mesh)
field = fem.FieldsMixed(region, n=3)

boundaries, loadcase = fem.dof.uniaxial(field, clamped=False)

umat = fem.ThreeFieldVariation(fem.NeoHookeCompressible(mu=1, lmbda=1e4))
solid = fem.SolidBody(umat, field)

move = fem.math.linsteps([0, 1], num=5)
step = fem.Step(items=[solid], ramp={boundaries["move"]: move}, boundaries=boundaries)

job = fem.CharacteristicCurve(steps=[step], boundary=boundaries["move"])
job.evaluate(tol=1e-2)
fig, ax = job.plot(
    xlabel="Displacement $u$ in mm $\longrightarrow$",
    ylabel="Normal Force $F$ in N $\longrightarrow$",
)

ax2 = field.imshow("Principal Values of Logarithmic Strain")

[adtzlr]

a) For the most general approach and without all the fancy automatic differentiation stuff, please have a look at Material. Unfortunately, I only have a minimal documentation docstring ready for mixed-fields yet.

b) For a more scikit-fem approach, there is also Form along with FormItem.

This should be really enough to read for the beginning 😄 .

[bhaveshshrimali]

Great @adtzlr ! This is super helpful!! Many thanks

[adtzlr]

And here you can find more or less the 1:1 translation of the original code. I replaced only the math-helpers. I used pypardiso as a faster alternative to the scipy-sparse solver (pip install pypardiso). Instead of Material, it is also possible to create the class by hand. There is nothing special except for the shape of the state variables, see Constitution.

import numpy as np
import felupe as fem
from felupe.math import det, inv, transpose, identity
import pypardiso


def F1(F, p, mu, lmbda):
    J = det(F)
    Finv = inv(F)
    return p * J * transpose(Finv) + mu * F


def F2(F, p, mu, lmbda):
    J = det(F)
    Js = 0.5 * (lmbda + p + 2.0 * np.sqrt(lmbda * mu + 0.25 * (lmbda + p) ** 2)) / lmbda
    dJsdp = (
        0.25 * lmbda + 0.25 * p + 0.5 * np.sqrt(lmbda * mu + 0.25 * (lmbda + p) ** 2)
    ) / (lmbda * np.sqrt(lmbda * mu + 0.25 * (lmbda + p) ** 2))
    return J - (Js + (p + mu / Js - lmbda * (Js - 1)) * dJsdp)


def A11(F, p, mu, lmbda):
    eye = identity(F)
    J = det(F)
    Finv = inv(F)
    L = (
        p * J * np.einsum("lk...,ji...->ijkl...", Finv, Finv)
        - p * J * np.einsum("jk...,li...->ijkl...", Finv, Finv)
        + mu * np.einsum("ik...,jl...->ijkl...", eye, eye)
    )
    return L


def A12(F, p, mu, lmbda):
    J = det(F)
    Finv = inv(F)
    return J * transpose(Finv)


def A22(F, p, mu, lmbda):
    Js = 0.5 * (lmbda + p + 2.0 * np.sqrt(lmbda * mu + 0.25 * (lmbda + p) ** 2)) / lmbda
    dJsdp = (
        0.25 * lmbda + 0.25 * p + 0.5 * np.sqrt(lmbda * mu + 0.25 * (lmbda + p) ** 2)
    ) / (lmbda * np.sqrt(lmbda * mu + 0.25 * (lmbda + p) ** 2))
    d2Jdp2 = 0.25 * mu / (lmbda * mu + 0.25 * (lmbda + p) ** 2) ** (3 / 2)
    L = (
        -2.0 * dJsdp
        - p * d2Jdp2
        + mu / Js**2 * dJsdp**2
        - mu / Js * d2Jdp2
        + lmbda * (Js - 1.0) * d2Jdp2
        + lmbda * dJsdp**2
    )
    return L


def stress(x, mu, lmbda):
    F, p, statevars = x[0], x[1], x[-1]
    return [F1(F, p, mu, lmbda), F2(F, p, mu, lmbda), statevars]


def elasticity(x, mu, lmbda):
    F, p = x[0], x[1]
    return [A11(F, p, mu, lmbda), A12(F, p, mu, lmbda), A22(F, p, mu, lmbda)]


mesh = fem.Cube(n=5).triangulate().add_midpoints_edges()

region_u = fem.RegionQuadraticTetra(mesh=mesh)
region_p = fem.RegionTetra(mesh=fem.mesh.dual(mesh, disconnect=False), quadrature=region_u.quadrature, grad=False)

displacement = fem.Field(region_u, dim=3)
pressure = fem.Field(region_p)
field = fem.FieldContainer([displacement, pressure])

boundaries, loadcase = fem.dof.uniaxial(field, clamped=False)
umat = fem.Material(stress, elasticity, mu=1, lmbda=1e4)
solid = fem.SolidBody(umat, field)

move = fem.math.linsteps([0, 1], num=5)
step = fem.Step(items=[solid], ramp={boundaries["move"]: move}, boundaries=boundaries)

job = fem.CharacteristicCurve(steps=[step], boundary=boundaries["move"])
job.evaluate(tol=1e-2, solver=pypardiso.spsolve)
fig, ax = job.plot(
    xlabel="Displacement $u$ in mm $\longrightarrow$",
    ylabel="Normal Force $F$ in N $\longrightarrow$",
)

ax2 = field.imshow("Principal Values of Logarithmic Strain")

[adtzlr]

FieldsMixed(region, n=2) for RegionQuadraticTetra is just a template for the following lines of code:

mesh = fem.Cube(n=5).triangulate().add_midpoints_edges()

region_u = fem.RegionQuadraticTetra(mesh=mesh)
region_p = fem.RegionTetra(mesh=fem.mesh.dual(mesh, disconnect=False), quadrature=region_u.quadrature, grad=False)

displacement = fem.Field(region_u, dim=3)
pressure = fem.Field(region_p)
field = fem.FieldContainer([displacement, pressure])

Note: The dual mesh for the secondary field $p$ may be disconnect-ed or not, this is your choice. It seems to be more stable if it is not disconnected here.

A german Wikipedia article contains a note on the instability when dealing with a disconnected mesh/field for the pressure. @bhaveshshrimali Do you know more on that topic?

[bhaveshshrimali]

As a side note, the assembly is remarkably fast in felupe. I will open separate discussions as I grok through the code, but running example 36 for a similar mesh was super-fast (even with no numba etc.). Great work on this. Thanks!

[adtzlr]

As a side note, the assembly is remarkably fast in felupe. I will open separate discussions as I grok through the code, but running example 36 for a similar mesh was super-fast (even with no numba etc.). Great work on this. Thanks!

Sounds great. But please compare mesh points /dof, PyVista shows all internal edges too and the mesh could look finer than it is!

See pyvista/pyvista#867 (comment)

[bhaveshshrimali]

I used mesh.ncells as a metric to compare number of cells, and visualized in Paraview.

[adtzlr]

That are indeed good news 👍🏻. I'm quite happy with assembly runtimes now. However, the overall memory footprint is still quite massive and lots of substeps take their time to complete because of array (re-)creations / allocations. I try to use the out-argument of numpy-functions when possible but there is still a lot to improve.

[adtzlr]

Q1-Q0 does not satisfy inf-sup, right ?

Yes, at least from what I read and understand 😄 on that topic, you're right.