.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/ex09_numeric-continuation.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_examples_ex09_numeric-continuation.py: 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. .. topic:: Numeric continuation of a hyperelastic cube. * use FElupe with contique * on-the-fly XDMF-file export * plot force-displacement curve .. admonition:: This example requires external packages. :class: hint .. code-block:: pip install contique matadi .. GENERATED FROM PYTHON SOURCE LINES 22-30 .. code-block:: Python import contique import matadi as mat import matplotlib.pyplot as plt import meshio import numpy as np import felupe as fem .. GENERATED FROM PYTHON SOURCE LINES 32-35 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. .. GENERATED FROM PYTHON SOURCE LINES 35-63 .. code-block:: Python # 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) .. image-sg:: /examples/images/sphx_glr_ex09_numeric-continuation_001.png :alt: Material Formulation: Yeoh :srcset: /examples/images/sphx_glr_ex09_numeric-continuation_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 64-68 An external normal force is applied at :math:`x=1` on a quarter model of a cube with symmetry planes at :math:`x=y=z=0`. Therefore, we have to define an external load vector which will be scaled with the load-proportionality factor :math:`\lambda` during numeric continuation. .. GENERATED FROM PYTHON SOURCE LINES 68-76 .. code-block:: Python # 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) .. GENERATED FROM PYTHON SOURCE LINES 77-79 The next step involves the problem definition for contique. For details have a look at `contique's README `_. .. GENERATED FROM PYTHON SOURCE LINES 79-117 .. code-block:: Python 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] .. GENERATED FROM PYTHON SOURCE LINES 118-121 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. .. GENERATED FROM PYTHON SOURCE LINES 121-145 .. code-block:: Python # 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]) .. rst-class:: sphx-glr-script-out .. code-block:: none |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#) | | .. GENERATED FROM PYTHON SOURCE LINES 146-149 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. .. GENERATED FROM PYTHON SOURCE LINES 149-156 .. code-block:: Python 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() .. tab-set:: .. tab-item:: Static Scene .. image-sg:: /examples/images/sphx_glr_ex09_numeric-continuation_002.png :alt: ex09 numeric continuation :srcset: /examples/images/sphx_glr_ex09_numeric-continuation_002.png :class: sphx-glr-single-img .. tab-item:: Interactive Scene .. offlineviewer:: /home/docs/checkouts/readthedocs.org/user_builds/felupe/checkouts/latest/docs/examples/images/sphx_glr_ex09_numeric-continuation_002.vtksz .. image-sg:: /examples/images/sphx_glr_ex09_numeric-continuation_003.png :alt: ex09 numeric continuation :srcset: /examples/images/sphx_glr_ex09_numeric-continuation_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 4.956 seconds) .. _sphx_glr_download_examples_ex09_numeric-continuation.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: ex09_numeric-continuation.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: ex09_numeric-continuation.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_