Twisting column

[adtzlr]

This is the twisting column example from Bonet & Wood.

import felupe as fem
import numpy as np
import pypardiso

vertical = fem.Cube(b=(1, 1, 6), n=(5, 5, 25))
horizontal = fem.Cube(a=(-1, 0, 6), b=(2, 1, 7), n=(13, 5, 5))

mesh = fem.MeshContainer([vertical, horizontal], merge=True).stack()
region = fem.RegionHexahedron(mesh)
field = fem.FieldContainer([fem.Field(region, dim=3)])

boundaries = fem.dof.BoundaryDict(fixed=fem.Boundary(field[0], fz=0))
umat = fem.NeoHookeCompressible(mu=100, lmbda=100)
solid = fem.SolidBody(umat=umat, field=field)

left = np.logical_and.reduce([mesh.x <= 0.5, mesh.y == 0, mesh.z >= 6])
right = np.logical_and.reduce([mesh.x >= 0.5, mesh.y == 1, mesh.z >= 6])
mask = np.logical_or(left, right)

region_boundary = fem.RegionHexahedronBoundary(mesh, mask=mask)
field_boundary = fem.Field(region_boundary, dim=3).as_container()
pressure = fem.SolidBodyPressure(field_boundary)

p = fem.math.linsteps([0, 3.2], num=20)
step = fem.Step(items=[solid, pressure], ramp={pressure: p}, boundaries=boundaries)

job = fem.Job(steps=[step])
job.evaluate(filename="result.xdmf", solver=pypardiso.spsolve)

ax = field.imshow(color="white", nonlinear_subdivision=2)

[adtzlr]

This one allows for much more rotations, using contique.

import felupe as fem
import numpy as np
import contique

vertical = fem.Cube(b=(1, 1, 6), n=(5, 5, 25))
horizontal = fem.Cube(a=(-1, 0, 6), b=(2, 1, 7), n=(13, 5, 5))

mesh = fem.MeshContainer([vertical, horizontal], merge=True).stack()
region = fem.RegionHexahedron(mesh)
field = fem.FieldContainer([fem.Field(region, dim=3)])

boundaries = fem.dof.BoundaryDict(fixed=fem.Boundary(field[0], fz=0))
umat = fem.NeoHookeCompressible(mu=100, lmbda=100)
solid = fem.SolidBody(umat=umat, field=field)

left = np.logical_and.reduce([mesh.x <= 0.5, mesh.y == 0, mesh.z >= 6])
right = np.logical_and.reduce([mesh.x >= 0.5, mesh.y == 1, mesh.z >= 6])
mask = np.logical_or(left, right)

region_boundary = fem.RegionHexahedronBoundary(mesh, mask=mask)
field_boundary = fem.Field(region_boundary, dim=3).as_container()
pressure = fem.SolidBodyPressure(field_boundary)

dof0, dof1 = fem.dof.partition(field, boundaries)

def fun(x, lpf, *args):
    "The system vector of equilibrium equations."

    # re-create field-values from active degrees of freedom
    solid.field[0].values.fill(0)
    solid.field[0].values.ravel()[dof1] += x
    pressure.update(lpf)

    return fem.tools.fun(items=[solid, pressure], x=solid.field)[dof1]


def dfundx(x, lpf, *args):
    """The jacobian of the system vector of equilibrium equations w.r.t. the
    primary unknowns."""

    solid.field[0].values.fill(0)
    solid.field[0].values.ravel()[dof1] += x
    pressure.update(lpf)

    r = fem.tools.fun(items=[solid, pressure], x=solid.field)
    K = fem.tools.jac(items=[solid, pressure], x=solid.field)

    return fem.solve.partition(solid.field, K, dof1, dof0, r)[2]


def dfundl(x, lpf, *args):
    """The jacobian of the system vector of equilibrium equations w.r.t. the
    load proportionality factor."""

    pressure.update(1.0)

    return fem.tools.fun(items=[pressure], x=solid.field)[dof1]


plotter = field.plot(color="white", off_screen=True)
plotter.open_gif("result.gif", fps=25)
plotter.camera.zoom(0.8)


def write_frame(step, res):
    plotter.mesh.points[:] = region.mesh.points + field[0].values
    plotter.write_frame()

# run contique (with rebalanced steps, 0% overshoot and a callback function)
Res = contique.solve(
    fun=fun,
    jac=[dfundx, dfundl],
    x0=field[0][dof1],
    lpf0=0,
    dxmax=0.5,
    dlpfmax=0.5,
    maxsteps=100,
    rebalance=True,
    # overshoot=1.05,
    callback=write_frame,
    tol=1e-6,
)

X = np.array([res.x for res in Res])

plotter.close()