.. 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 example 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 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 a force-displacement curve .. GENERATED FROM PYTHON SOURCE LINES 20-29 .. code-block:: Python import contique import matplotlib.pyplot as plt import meshio import numpy as np import felupe as fem import felupe.constitution.tensortrax as mat .. GENERATED FROM PYTHON SOURCE LINES 31-33 First, setup a problem as usual (mesh, region, field, boundaries and umat). The unstable material behavior is plotted for uniaxial incompressible tension. .. GENERATED FROM PYTHON SOURCE LINES 33-50 .. 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 boundaries = fem.dof.symmetry(field[0], axes=(True, True, True)) # partition degrees of freedom dof0, dof1 = fem.dof.partition(field, boundaries) # constitutive isotropic hyperelastic material formulation yeoh = mat.Hyperelastic(mat.models.hyperelastic.yeoh, C10=0.5, C20=-0.25, C30=0.025) ax = yeoh.plot(incompressible=True, ux=np.linspace(1, 2.76), bx=None, ps=None) solid = fem.SolidBodyNearlyIncompressible(yeoh, field, bulk=5000) .. image-sg:: /examples/images/sphx_glr_ex09_numeric-continuation_001.png :alt: Hyperelastic (Yeoh), C10=0.5, C20=-0.25, C30=0.025 :srcset: /examples/images/sphx_glr_ex09_numeric-continuation_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 51-55 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 55-63 .. code-block:: Python # external force vector at x=1 right = region.mesh.points[:, 0] == 1 v = 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 64-66 The next step involves the problem definition for contique. For details have a look at `contique's README `_. .. GENERATED FROM PYTHON SOURCE LINES 66-100 .. code-block:: Python def fun(x, lpf, *args): "The system vector of equilibrium equations." # re-create field-values from active degrees of freedom field[0].values.ravel()[dof1] = x load.update(values_load * lpf) return fem.tools.fun(items=[solid, load], x=field)[dof1] def dfundx(x, lpf, *args): """The jacobian of the system vector of equilibrium equations w.r.t. the primary unknowns.""" field[0].values.ravel()[dof1] = x load.update(values_load * lpf) r = fem.tools.fun(items=[solid, load], x=field) K = fem.tools.jac(items=[solid, load], x=field) return fem.solve.partition(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.""" load.update(values_load) return fem.tools.fun(items=[load], x=field)[dof1] .. GENERATED FROM PYTHON SOURCE LINES 101-104 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 104-132 .. 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.06, dlpfmax=0.02, maxsteps=35, rebalance=True, overshoot=1.05, callback=step_to_xdmf, tol=1e-2, ) X = np.array([res.x for res in Res]) # check the final lpf value assert np.isclose(X[-1, 1], -0.3982995) .. rst-class:: sphx-glr-script-out .. code-block:: none |Step,C.| Control Component | Norm (Iter.#) | Message | |-------|-------------------|---------------|-------------| | 1,1 | 12+ => 12+ | 1.1e-03 ( 2#) | | | 2,1 | 12+ => 12+ | 1.3e-03 ( 2#) | | | 3,1 | 12+ => 12+ | 1.5e-03 ( 2#) | | | 4,1 | 12+ => 12+ | 7.8e-03 ( 2#) | | | 5,1 | 12+ => 12+ | 1.5e-04 ( 4#) | | | 6,1 | 12+ => 5- | 5.4e+02 ( 8#) |Failed | | 7,1 | 12+ => 12+ | 1.1e-01 ( 8#) |Failed | | 8,1 | 12+ => 2+ | 3.2e-03 ( 2#) | => re-Cycle | | 2 | 2+ => 2+ | 1.4e-03 ( 2#) | | | 9,1 | 2+ => 9+ | 1.3e-03 ( 2#) |tol.Overshoot| | 10,1 | 9+ => 2+ | 1.1e-03 ( 2#) |tol.Overshoot| | 11,1 | 2+ => 5+ | 3.6e-03 ( 2#) |tol.Overshoot| | 12,1 | 5+ => 9+ | 1.8e-08 ( 3#) |tol.Overshoot| | 13,1 | 9+ => 2+ | 1.1e-07 ( 3#) |tol.Overshoot| | 14,1 | 2+ => 12- | 6.1e-07 ( 3#) | => re-Cycle | | 2 | 12- => 12- | 3.4e-03 ( 4#) | | | 15,1 | 12- => 12- | 1.5e-04 ( 4#) | | | 16,1 | 12- => 12- | 8.9e-06 ( 4#) | | | 17,1 | 12- => 12- | 4.0e-06 ( 4#) | | | 18,1 | 12- => 12- | 5.4e-03 ( 3#) | | | 19,1 | 12- => 12- | 8.5e-03 ( 4#) | | | 20,1 | 12- => 12- | 1.0e+02 ( 8#) |Failed | | 21,1 | 12- => 5- | 8.5e+07 ( 8#) |Failed | | 22,1 | 12- => 2+ | 6.3e-01 ( 8#) |Failed | | 23,1 | 12- => 12- | 8.3e-03 ( 3#) | | | 24,1 | 12- => 12- | 1.4e+01 ( 8#) |Failed | | 25,1 | 12- => 0+ | 2.6e-03 ( 3#) | => re-Cycle | | 2 | 0+ => 0+ | 6.0e-05 ( 2#) | | | 26,1 | 0+ => 5+ | 5.7e-05 ( 2#) |tol.Overshoot| | 27,1 | 5+ => 2+ | 5.5e-05 ( 2#) |tol.Overshoot| | 28,1 | 2+ => 2+ | 1.9e-04 ( 2#) | | | 29,1 | 2+ => 9+ | 6.3e-04 ( 2#) |tol.Overshoot| | 30,1 | 9+ => 12+ | 2.1e-03 ( 2#) | => re-Cycle | | 2 | 12+ => 12+ | 1.9e-05 ( 4#) | | | 31,1 | 12+ => 12+ | 6.6e-04 ( 3#) | | | 32,1 | 12+ => 12+ | 1.8e-04 ( 3#) | | | 33,1 | 12+ => 12+ | 8.3e-05 ( 3#) | | | 34,1 | 12+ => 12+ | 5.4e-05 ( 3#) | | | 35,1 | 12+ => 12+ | 4.2e-05 ( 3#) | | .. GENERATED FROM PYTHON SOURCE LINES 133-136 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. .. GENERATED FROM PYTHON SOURCE LINES 136-143 .. 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/stable/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 2.416 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 ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: ex09_numeric-continuation.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_