Note
Go to the end to download the full example code.
Numeric Continuation#
With the help of contique (install with pip install contique) it is possible to apply a numerical parameter continuation algorithm on any system of equilibrium equations. This advanced tutorial demonstrates the usage of FElupe in conjunction with contique. An unstable isotropic hyperelastic material formulation is applied on a single hexahedron. The model will be visualized by the XDMF-output (of meshio) and the resulting force - displacement curve will be plotted.
This example requires external packages.
pip install contique matadi
import contique
import matadi as mat
import matplotlib.pyplot as plt
import meshio
import numpy as np
import felupe as fem
First, setup a problem as usual (mesh, region, field, boundaries and umat). For the
material definition we use matADi (pip install madati). We utilize matADi’s lab-
capability to visualize the unstable material behavior in uniaxial tension.
# setup a numeric region on a cube
mesh = fem.Cube(n=2)
region = fem.RegionHexahedron(mesh)
field = fem.FieldContainer([fem.Field(region, dim=3)])
# introduce symmetry planes at x=y=z=0
bounds = fem.dof.symmetry(field[0], axes=(True, True, True))
# partition degrees of freedom
dof0, dof1 = fem.dof.partition(field, bounds)
# constitutive isotropic hyperelastic material formulation
yeoh = mat.MaterialHyperelastic(mat.models.yeoh, C10=0.5, C20=-0.25, C30=0.025, bulk=5)
lab = mat.Lab(yeoh)
data = lab.run(
ux=True,
bx=False,
ps=False,
shear=False,
stretch_min=1,
stretch_max=2.75,
num=100,
)
lab.plot(data)
body = fem.SolidBody(yeoh, field)
An external normal force is applied at \(x=1\) on a quarter model of a cube with symmetry planes at \(x=y=z=0\). Therefore, we have to define an external load vector which will be scaled with the load-proportionality factor \(\lambda\) during numeric continuation.
# external force vector at x=1
right = region.mesh.points[:, 0] == 1
v = 0.01 * region.mesh.cells_per_point[right]
values_load = np.vstack([v, np.zeros_like(v), np.zeros_like(v)]).T
load = fem.PointLoad(field, right, values_load)
The next step involves the problem definition for contique. For details have a look at contique’s README.
def fun(x, lpf, *args):
"The system vector of equilibrium equations."
# re-create field-values from active degrees of freedom
body.field[0].values.fill(0)
body.field[0].values.ravel()[dof1] += x
load.update(values_load * lpf)
return fem.tools.fun([body, load], body.field)[dof1]
def dfundx(x, lpf, *args):
"""The jacobian of the system vector of equilibrium equations w.r.t. the
primary unknowns."""
body.field[0].values.fill(0)
body.field[0].values.ravel()[dof1] += x
load.update(values_load * lpf)
r = fem.tools.fun([body, load], body.field, True)
K = fem.tools.jac([body, load], body.field, True)
return fem.solve.partition(body.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."""
body.field[0].values.fill(0)
body.field[0].values.ravel()[dof1] += x
load.update(values_load)
return load.assemble.vector()[dof1]
Next we have to init the problem and specify the initial values of unknowns (the undeformed configuration). After each completed step of the numeric continuation the XDMF-file will be updated.
# write xdmf file during numeric continuation
with meshio.xdmf.TimeSeriesWriter("result.xdmf") as writer:
writer.write_points_cells(mesh.points, [(mesh.cell_type, mesh.cells)])
def step_to_xdmf(step, res):
writer.write_data(step, point_data={"u": field[0].values})
# run contique (w/ rebalanced steps, 5% overshoot and a callback function)
Res = contique.solve(
fun=fun,
jac=[dfundx, dfundl],
x0=field[0][dof1],
lpf0=0,
dxmax=0.05,
dlpfmax=0.5,
maxsteps=80,
rebalance=True,
overshoot=1.05,
callback=step_to_xdmf,
)
X = np.array([res.x for res in Res])
|Step,C.| Control Component | Norm (Iter.#) | Message |
|-------|-------------------|---------------|-------------|
| 1,1 | 12+ => 12+ | 6.9e-09 ( 2#) | |
| 2,1 | 12+ => 12+ | 7.7e-09 ( 2#) | |
| 3,1 | 12+ => 12+ | 8.7e-09 ( 2#) | |
| 4,1 | 12+ => 12+ | 3.6e-08 ( 2#) | |
| 5,1 | 12+ => 12+ | 1.6e-07 ( 2#) | |
| 6,1 | 12+ => 12+ | 8.1e-07 ( 2#) | |
| 7,1 | 12+ => 12+ | 5.9e-11 ( 3#) | |
| 8,1 | 12+ => 12+ | 2.9e-09 ( 3#) | |
| 9,1 | 12+ => 12+ | 8.1e-07 ( 3#) | |
| 10,1 | 12+ => 12+ | 3.3e-01 ( 8#) |Failed |
| 11,1 | 12+ => 12+ | 7.8e-08 ( 4#) | |
| 12,1 | 12+ => 12+ | 4.1e-02 ( 8#) |Failed |
| 13,1 | 12+ => 2+ | 3.9e-02 ( 8#) |Failed |
| 14,1 | 12+ => 0+ | 7.5e+01 ( 8#) |Failed |
| 15,1 | 12+ => 5+ | 7.3e-07 ( 3#) | => re-Cycle |
| 2 | 5+ => 5+ | 2.4e-08 ( 2#) | |
| 16,1 | 5+ => 5+ | 1.5e-08 ( 2#) | |
| 17,1 | 5+ => 5+ | 8.1e-09 ( 2#) | |
| 18,1 | 5+ => 9+ | 1.0e-08 ( 2#) |tol.Overshoot|
| 19,1 | 9+ => 12- | 1.6e-11 ( 2#) | => re-Cycle |
| 2 | 12- => 12- | 3.3e-11 ( 4#) | |
| 20,1 | 12- => 12- | 2.7e-08 ( 3#) | |
| 21,1 | 12- => 12- | 5.6e-09 ( 3#) | |
| 22,1 | 12- => 12- | 2.4e-09 ( 3#) | |
| 23,1 | 12- => 12- | 1.5e-09 ( 3#) | |
| 24,1 | 12- => 12- | 1.2e-09 ( 3#) | |
| 25,1 | 12- => 12- | 1.1e-09 ( 3#) | |
| 26,1 | 12- => 12- | 1.1e-09 ( 3#) | |
| 27,1 | 12- => 12- | 1.0e-09 ( 3#) | |
| 28,1 | 12- => 12- | 5.2e-10 ( 3#) | |
| 29,1 | 12- => 12- | 3.7e-11 ( 3#) | |
| 30,1 | 12- => 12- | 1.8e-11 ( 3#) | |
| 31,1 | 12- => 12- | 7.6e-11 ( 3#) | |
| 32,1 | 12- => 12- | 3.3e-10 ( 3#) | |
| 33,1 | 12- => 12- | 1.9e-09 ( 3#) | |
| 34,1 | 12- => 12- | 2.0e-08 ( 3#) | |
| 35,1 | 12- => 12- | 3.0e-10 ( 4#) | |
| 36,1 | 12- => 0+ | 2.7e+03 ( 8#) |Failed |
| 37,1 | 12- => 5- | 2.6e-02 ( 8#) |Failed |
| 38,1 | 12- => 12- | 4.6e-03 ( 8#) |Failed |
| 39,1 | 12- => 12- | 5.0e-08 ( 3#) | |
| 40,1 | 12- => 9- | 2.0e-01 ( 8#) |Failed |
| 41,1 | 12- => 12- | 1.3e-07 ( 3#) | |
| 42,1 | 12- => 9+ | 3.6e-02 ( 8#) |Failed |
| 43,1 | 12- => 12- | 2.4e-09 ( 4#) | |
| 44,1 | 12- => 12- | 3.1e-03 ( 8#) |Failed |
| 45,1 | 12- => 12- | 1.4e-03 ( 8#) |Failed |
| 46,1 | 12- => 0+ | 3.7e-03 ( 8#) |Failed |
| 47,1 | 12- => 12- | 1.7e-08 ( 3#) | |
| 48,1 | 12- => 0+ | 3.1e-04 ( 8#) |Failed |
| 49,1 | 12- => 0+ | 7.4e-07 ( 3#) | => re-Cycle |
| 2 | 0+ => 5+ | 5.8e-09 ( 2#) |tol.Overshoot|
| 50,1 | 5+ => 9+ | 5.1e-09 ( 2#) |tol.Overshoot|
| 51,1 | 9+ => 9+ | 4.5e-09 ( 2#) | |
| 52,1 | 9+ => 0+ | 1.4e-08 ( 2#) |tol.Overshoot|
| 53,1 | 0+ => 0+ | 4.2e-08 ( 2#) | |
| 54,1 | 0+ => 0+ | 1.2e-07 ( 2#) | |
| 55,1 | 0+ => 12+ | 3.1e-07 ( 2#) | => re-Cycle |
| 2 | 12+ => 12+ | 2.0e-08 ( 3#) | |
| 56,1 | 12+ => 12+ | 1.6e-09 ( 3#) | |
| 57,1 | 12+ => 12+ | 4.5e-10 ( 3#) | |
| 58,1 | 12+ => 12+ | 2.3e-10 ( 3#) | |
| 59,1 | 12+ => 12+ | 1.6e-10 ( 3#) | |
| 60,1 | 12+ => 12+ | 1.3e-10 ( 3#) | |
| 61,1 | 12+ => 12+ | 1.3e-10 ( 3#) | |
| 62,1 | 12+ => 12+ | 1.3e-10 ( 3#) | |
| 63,1 | 12+ => 12+ | 1.4e-10 ( 3#) | |
| 64,1 | 12+ => 12+ | 1.6e-10 ( 3#) | |
| 65,1 | 12+ => 12+ | 1.9e-10 ( 3#) | |
| 66,1 | 12+ => 12+ | 2.2e-10 ( 3#) | |
| 67,1 | 12+ => 12+ | 2.5e-10 ( 3#) | |
| 68,1 | 12+ => 12+ | 2.9e-10 ( 3#) | |
| 69,1 | 12+ => 12+ | 3.1e-10 ( 3#) | |
| 70,1 | 12+ => 12+ | 3.3e-10 ( 3#) | |
| 71,1 | 12+ => 12+ | 3.3e-10 ( 3#) | |
| 72,1 | 12+ => 12+ | 2.7e-10 ( 3#) | |
| 73,1 | 12+ => 12+ | 9.9e-11 ( 3#) | |
| 74,1 | 12+ => 12+ | 8.2e-11 ( 3#) | |
| 75,1 | 12+ => 12+ | 4.8e-11 ( 3#) | |
| 76,1 | 12+ => 12+ | 2.4e-11 ( 3#) | |
| 77,1 | 12+ => 12+ | 1.2e-11 ( 3#) | |
| 78,1 | 12+ => 12+ | 6.3e-12 ( 3#) | |
| 79,1 | 12+ => 12+ | 3.5e-12 ( 3#) | |
| 80,1 | 12+ => 12+ | 2.1e-12 ( 3#) | |
Finally, the force-displacement curve is plotted. It can be seen that the resulting (unstable) force-controlled equilibrium path is equal to the displacement-controlled load case of matADi’s lab.
plt.figure()
plt.plot(X[:, 0], X[:, -1], "x-")
plt.xlabel(r"displacement $u(x=1)/L$ $\longrightarrow$")
plt.ylabel(r"load-proportionality-factor $\lambda$ $\longrightarrow$")
field.plot("Displacement", component=0).show()
Total running time of the script: (0 minutes 4.677 seconds)