.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "tutorial/extut04_building_blocks.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_tutorial_extut04_building_blocks.py: Building Blocks --------------- .. topic:: Learn the building blocks of FElupe. * create a meshed cube with hexahedron elements * setup your own numeric region with a mesh, an element and a quadrature * add a displacement field to a field container * define your own Neo-Hookean material formulation * apply your own boundary conditions * solve the problem (create your own Newton-Rhapson iteration loop) * export the displaced mesh along with cauchy stress projected to mesh-points Start setting up a problem in FElupe by the creation of a numeric :class:`~felupe.Region` with a geometry :class:`~felupe.Mesh`, a finite **Element** formulation :class:`~felupe.Hexahedron` and a **Quadrature** rule :class:`~felupe.GaussLegendre`. .. image:: examples/extut03_building_blocks_sketch.svg :width: 600px .. GENERATED FROM PYTHON SOURCE LINES 26-35 .. code-block:: Python import felupe as fem mesh = fem.Cube(n=6) element = fem.Hexahedron() quadrature = fem.GaussLegendre(order=1, dim=3) mesh.plot().show() .. tab-set:: .. tab-item:: Static Scene .. image-sg:: /tutorial/images/sphx_glr_extut04_building_blocks_001.png :alt: extut04 building blocks :srcset: /tutorial/images/sphx_glr_extut04_building_blocks_001.png :class: sphx-glr-single-img .. tab-item:: Interactive Scene .. offlineviewer:: /home/docs/checkouts/readthedocs.org/user_builds/felupe/checkouts/stable/docs/tutorial/images/sphx_glr_extut04_building_blocks_001.vtksz .. GENERATED FROM PYTHON SOURCE LINES 37-38 An overview of the mesh-metadata is available in the console. .. GENERATED FROM PYTHON SOURCE LINES 38-40 .. code-block:: Python print(mesh) .. rst-class:: sphx-glr-script-out .. code-block:: none Number of points: 216 Number of cells: hexahedron: 125 .. GENERATED FROM PYTHON SOURCE LINES 41-51 Region ~~~~~~ A :class:`~felupe.Region` essentially pre-calculates element shape/ansatz/basis functions and derivatives evaluated at every quadrature point of every cell w.r.t. the undeformed coordinates (as attribute :attr:`~felupe.Region.dhdX`). .. note:: By using a template region like :class:`~felupe.RegionHexahedron`, only the mesh has to be created. The element formulation and the appropriate quadrature scheme are chosen automatically. .. GENERATED FROM PYTHON SOURCE LINES 51-55 .. code-block:: Python region = fem.Region(mesh, element, quadrature) # region = fem.RegionHexahedron(mesh) print(region) .. rst-class:: sphx-glr-script-out .. code-block:: none Element formulation: Hexahedron Quadrature rule: GaussLegendre Gradient evaluated: True Hessian evaluated: False .. GENERATED FROM PYTHON SOURCE LINES 56-58 The scheme of the region, i.e. the element with its point labels and the integration points of the quadrature rule may be plotted. .. GENERATED FROM PYTHON SOURCE LINES 58-60 .. code-block:: Python region.plot().show() .. tab-set:: .. tab-item:: Static Scene .. image-sg:: /tutorial/images/sphx_glr_extut04_building_blocks_002.png :alt: extut04 building blocks :srcset: /tutorial/images/sphx_glr_extut04_building_blocks_002.png :class: sphx-glr-single-img .. tab-item:: Interactive Scene .. offlineviewer:: /home/docs/checkouts/readthedocs.org/user_builds/felupe/checkouts/stable/docs/tutorial/images/sphx_glr_extut04_building_blocks_002.vtksz .. rst-class:: sphx-glr-script-out .. code-block:: none /home/docs/checkouts/readthedocs.org/user_builds/felupe/envs/stable/lib/python3.12/site-packages/felupe/element/_base.py:164: PyVistaFutureWarning: The default value of `algorithm` for the filter `UnstructuredGrid.extract_surface` will change in the future. It currently defaults to `'dataset_surface'`, but will change to `None`. Explicitly set the `algorithm` keyword to silence this warning. view.mesh = view.mesh.extract_surface().extract_feature_edges() .. GENERATED FROM PYTHON SOURCE LINES 61-65 An array containing products of quadrature weights multiplied by the determinants of the (geometric) jacobians is stored as the array of (undeformed) differential volumes :attr:`~felupe.Region.dV`. The sum of all differential volumes gives the total (undeformed) volume of the region. .. GENERATED FROM PYTHON SOURCE LINES 65-68 .. code-block:: Python dV = region.dV V = dV.sum() .. GENERATED FROM PYTHON SOURCE LINES 69-84 Field ~~~~~ In a next step, a displacement :class:`~felupe.Field` is added to the :class:`~felupe.Region`. This may be either a scalar- or a vector-valued field. The values at mesh-points are obtained with the attribute :attr:`~felupe.Field.values`. Interpolated field values at quadrature points are calculated with the :meth:`~felupe.Field.interpolate` method. Additionally, the displacement gradient w.r.t. the undeformed coordinates is calculated for every quadrature point of every cell in the region with the field method :meth:`~felupe.Field.grad`. A generalized extraction method :meth:`extract(grad=True, add_identity=True, sym=False) ` allows several arguments to be passed. This involves whether the gradient or the values are extracted. If the gradient is extracted, the identity matrix may be added to the gradient (useful for the calculation of the deformation gradient). Optionally, the symmetric part is returned (small strain tensor). .. GENERATED FROM PYTHON SOURCE LINES 84-91 .. code-block:: Python displacement = fem.Field(region, dim=3) u = displacement.values v = displacement.interpolate() dudX = displacement.grad() .. GENERATED FROM PYTHON SOURCE LINES 92-97 Next, the field is added to a :class:`~felupe.FieldContainer`, which holds one or more fields. Like a :class:`~felupe.Field`, the :class:`~felupe.FieldContainer` also provides the :meth:`extract(grad=True, add_identity=True, sym=False) ` method, returning a list of interpolated field values or gradients. .. GENERATED FROM PYTHON SOURCE LINES 97-100 .. code-block:: Python field = fem.FieldContainer([displacement]) print(field) .. rst-class:: sphx-glr-script-out .. code-block:: none Number of fields: 1 Dimension of fields: Field: 3 .. GENERATED FROM PYTHON SOURCE LINES 101-106 The deformation gradient is obtained by a sum of the identity and the displacement gradient. The :meth:`~felupe.FieldContainer.extract`-method generates a list with extracted values at the integration points for each field in the field container. To evaluate the deformation gradient tensor of the first field of a container, the identity matrix is added to the gradient of the field. .. GENERATED FROM PYTHON SOURCE LINES 106-108 .. code-block:: Python F = field.extract(grad=True, sym=False, add_identity=True) .. GENERATED FROM PYTHON SOURCE LINES 109-111 Methods to evaluate the deformation gradient as well as strain measures are provided in :class:`FieldContainer.evaluate `. .. GENERATED FROM PYTHON SOURCE LINES 111-116 .. code-block:: Python import numpy as np log_strain = field.evaluate.strain(fun=lambda stretch: np.log(stretch), tensor=True) principal_stretches = field.evaluate.strain(fun=lambda stretch: stretch, tensor=False) .. GENERATED FROM PYTHON SOURCE LINES 117-144 Constitution ~~~~~~~~~~~~ The material behavior has to be provided by the first Piola-Kirchhoff stress tensor as a function of the deformation gradient. FElupe provides a comprehensive constitutive library (see :ref:`felupe-api-constitution`) including :class:`~felupe.NeoHooke`, :class:`~felupe.LinearElastic` and a generalized Hu-Washizu (u,p,J) :class:`~felupe.ThreeFieldVariation`. By alternative, an isotropic material formulation :class:`~felupe.constitution.tensortrax.Hyperelastic` is defined by a strain energy density function where both variation (stress) and linearization (elasticity) are carried out by automatic differentiation. The latter one is demonstrated here with a nearly-incompressible version of the Neo-Hookean material model formulation. .. note:: FElupe supports different backends for automatic differentiation and uses `tensortrax `_ by default. It is important to use only differentiable math-functions from the backend, e.g. NumPy-like modules `tensortrax.math `_, `tensortrax.math.linalg `_ or `tensortrax.math.special `_. `tensortrax `_ is a dependency of FElupe and hence, is installed along with `FElupe `_. .. math:: \psi = \frac{\mu}{2} \left( J^{-2/3}\ \text{tr}(\boldsymbol{C}) - 3 \right) + \frac{K}{2} \left( J - 1 \right)^2 .. GENERATED FROM PYTHON SOURCE LINES 144-162 .. code-block:: Python import tensortrax.math as tm import felupe.constitution.tensortrax as mat def W(C, mu, bulk): "Isotropic hyperelastic Neo-Hookean material formulation." J = tm.sqrt(tm.linalg.det(C)) return mu / 2 * (J ** (-2 / 3) * tm.trace(C) - 3) + bulk * (J - 1) ** 2 / 2 umat = mat.Hyperelastic(W, mu=1.0, bulk=2.0) P = umat.gradient A = umat.hessian .. GENERATED FROM PYTHON SOURCE LINES 163-171 Boundary Conditions ~~~~~~~~~~~~~~~~~~~ Next we enforce boundary conditions on the displacement field. Boundary conditions are stored as a dictionary of multiple :class:`~felupe.Boundary` instances. First, the left end of the cube is fixed. Displacements on the right end are fixed in directions y and z whereas displacements in direction x are prescribed with a user-defined ``value=0.5``. A :class:`~felupe.Boundary` holds useful attributes like :attr:`~felupe.Boundary.points` or :attr:`~felupe.Boundary.dof`. .. GENERATED FROM PYTHON SOURCE LINES 171-182 .. code-block:: Python f0 = lambda x: np.isclose(x, 0) f1 = lambda x: np.isclose(x, 1) boundaries = fem.BoundaryDict() boundaries["left"] = fem.Boundary(displacement, fx=f0) boundaries["right"] = fem.Boundary(displacement, fx=f1, skip=(1, 0, 0)) boundaries["move"] = fem.Boundary(displacement, fx=f1, skip=(0, 1, 1), value=0.5) boundaries.plot().show() .. tab-set:: .. tab-item:: Static Scene .. image-sg:: /tutorial/images/sphx_glr_extut04_building_blocks_003.png :alt: extut04 building blocks :srcset: /tutorial/images/sphx_glr_extut04_building_blocks_003.png :class: sphx-glr-single-img .. tab-item:: Interactive Scene .. offlineviewer:: /home/docs/checkouts/readthedocs.org/user_builds/felupe/checkouts/stable/docs/tutorial/images/sphx_glr_extut04_building_blocks_003.vtksz .. GENERATED FROM PYTHON SOURCE LINES 183-189 Partition of deegrees of freedom ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The separation of active and inactive degrees of freedom is performed by a so-called :func:`~felupe.dof.partition`. External values of prescribed displacement degrees of freedom are obtained by the application of the boundary values on the displacement field. .. GENERATED FROM PYTHON SOURCE LINES 189-192 .. code-block:: Python dof0, dof1 = fem.dof.partition(field, boundaries) ext0 = fem.dof.apply(field, boundaries, dof0) .. GENERATED FROM PYTHON SOURCE LINES 193-209 Integral forms of equilibrium equations ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The integral (or weak) forms of equilibrium equations are defined by the :class:`~felupe.IntegralForm` class, see Eq. :eq:`eq_integral_form`. The pre-evaluated function of interest has to be passed as the ``fun`` argument and the test field as the ``v`` argument. By setting ``grad_v=[True]`` (default), FElupe passes the gradient of the first test field to the integral form. FElupe assumes a linear form if ``u=None`` (default) and creates a bilinear form if a field is passed to the trial field argument ``u``. .. math:: :label: eq_integral_form \int_V P_{iJ} : \frac{\partial \delta u_i}{\partial X_J} \ dV \qquad \text{and} \qquad \int_V \frac{\partial \delta u_i} {\partial X_J} : \mathbb{A}_{iJkL} : \frac{\partial u_k}{\partial X_L} \ dV .. GENERATED FROM PYTHON SOURCE LINES 209-212 .. code-block:: Python linearform = fem.IntegralForm(P(F)[:-1], field, dV, grad_v=[True]) bilinearform = fem.IntegralForm(A(F), field, dV, u=field, grad_v=[True], grad_u=[True]) .. GENERATED FROM PYTHON SOURCE LINES 213-215 The assembly of both forms lead to the (point-based) internal force vector and the (sparse) stiffness matrix. .. GENERATED FROM PYTHON SOURCE LINES 215-218 .. code-block:: Python r = linearform.assemble() K = bilinearform.assemble() .. GENERATED FROM PYTHON SOURCE LINES 219-260 Prepare (partition) and solve the linearized equation system ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ In order to solve the linearized equation system a :func:`~felupe.solve.partition` into active and inactive degrees of freedom has to be performed, see Eqs. :eq:`eq_partition`. This system may then be passed to the (sparse direct) solver. Given a set of nonlinear equilibrium equations :math:`\boldsymbol{g}` the unknowns :math:`\boldsymbol{u}` are found by linearization at a valid initial state of equilibrium and an iterative Newton-Rhapson solution prodecure. The incremental values of inactive degrees of freedom are given as the difference of external prescribed and current values of unknowns. The (linear) solution is equal to the first result of a Newton-Rhapson iterative solution procedure. The solution ``du`` is finally added to the displacement field. .. math:: :label: eq_partition \boldsymbol{g}_1(\boldsymbol{u}) &= -\boldsymbol{r}_1(\boldsymbol{u}) + \boldsymbol{f}_1 \boldsymbol{g}_1(\boldsymbol{u} + d\boldsymbol{u}) &\approx -\boldsymbol{r}_1 + \boldsymbol{f}_1 - \frac{\partial \boldsymbol{r}_1}{\partial \boldsymbol{u}_1} \ d\boldsymbol{u}_1 - \frac{ \partial \boldsymbol{r}_1}{\partial \boldsymbol{u}_0} \ d\boldsymbol{u}_0 = \boldsymbol{0} d\boldsymbol{u}_0 &= \boldsymbol{u}_0^{(ext)} - \boldsymbol{u}_0 \text{solve} \qquad \boldsymbol{K}_{11}\ d\boldsymbol{u}_1 &= \boldsymbol{g}_1 - \boldsymbol{K}_{10}\ d\boldsymbol{u}_{0} \boldsymbol{u}_0 &+= d\boldsymbol{u}_0 \boldsymbol{u}_1 &+= d\boldsymbol{u}_1 The default solver of FElupe is `SuperLU `_ provided by the sparse package of `SciPy `_. A significantly faster alternative is `pypardiso `_ which may be installed from PyPI with ``pip install pypardiso`` (not included with FElupe). The optional argument ``solver`` of :func:`~felupe.solve.solve` accepts any user-defined solver. .. GENERATED FROM PYTHON SOURCE LINES 260-268 .. code-block:: Python from scipy.sparse.linalg import spsolve # default # from pypardiso import spsolve system = fem.solve.partition(field, K, dof1, dof0, r) dfield = fem.solve.solve(*system, ext0, solver=spsolve) # .reshape(*u.shape) # field += dfield .. GENERATED FROM PYTHON SOURCE LINES 269-270 A very simple Newton-Rhapson code looks like this: .. GENERATED FROM PYTHON SOURCE LINES 270-289 .. code-block:: Python for iteration in range(8): F = field.extract() linearform = fem.IntegralForm(P(F)[:-1], field, dV) bilinearform = fem.IntegralForm(A(F), field, dV, field) r = linearform.assemble() K = bilinearform.assemble() system = fem.solve.partition(field, K, dof1, dof0, r) dfield = fem.solve.solve(*system, ext0, solver=spsolve) norm = np.linalg.norm(dfield) print(iteration, norm) field += dfield if norm < 1e-12: break .. rst-class:: sphx-glr-script-out .. code-block:: none 0 4.525274987383307 1 0.1686563278412268 2 0.012102266767787915 3 6.444069678832645e-05 4 3.407185906351198e-09 5 3.0674161708530823e-16 .. GENERATED FROM PYTHON SOURCE LINES 290-292 By alternative, one may also use the :func:`Newton-Rhapson ` function of FElupe. .. GENERATED FROM PYTHON SOURCE LINES 292-299 .. code-block:: Python field[0].fill(0) solid = fem.SolidBody(umat, field) loadcase = {"dof1": dof1, "dof0": dof0, "ext0": ext0} res = fem.newtonrhapson(items=[solid], verbose=2, tol=1e-12, **loadcase) field = res.x .. rst-class:: sphx-glr-script-out .. code-block:: none Newton-Rhapson solver ===================== | # | norm(fun) | norm(dx) | |---|-----------|-----------| | 1 | 3.188e-01 | 4.525e+00 | | 2 | 2.719e-02 | 1.687e-01 | | 3 | 1.758e-04 | 1.210e-02 | | 4 | 6.337e-09 | 6.444e-05 | | 5 | 8.504e-16 | 3.407e-09 | Converged in 5 iterations (Assembly: 0.06196 s, Solve: 0.01328 s). .. GENERATED FROM PYTHON SOURCE LINES 300-302 All 3x3 components of the deformation gradient of integration point 2 of cell 1 (Python is 0-indexed) are obtained with .. GENERATED FROM PYTHON SOURCE LINES 302-305 .. code-block:: Python defgrad = field.evaluate.deformation_gradient() print(defgrad[:, :, 1, 0]) .. rst-class:: sphx-glr-script-out .. code-block:: none [[ 1.47073968 -0.029986 -0.029986 ] [ 0.23485302 0.90111657 0.00477373] [ 0.23485302 0.00477373 0.90111657]] .. GENERATED FROM PYTHON SOURCE LINES 306-310 Export and plot of results ~~~~~~~~~~~~~~~~~~~~~~~~~~ Results are exported as VTK or XDMF files using `meshio `_. .. GENERATED FROM PYTHON SOURCE LINES 310-312 .. code-block:: Python fem.save(region, field, filename="result.vtk") .. GENERATED FROM PYTHON SOURCE LINES 313-319 Any tensor at quadrature points, shifted or projected to and averaged at mesh-points, is evaluated for ``quad`` and ``hexahedron`` cell types by :func:`~felupe.topoints` or :func:`~felupe.project`, respectively. For example, the calculation of the Cauchy stress involves the conversion from the first Piola-Kirchhoff stress to the Cauchy stress followed by the shift or the projection. The stress results at mesh points are passed as a dictionary to the ``point_data`` argument. .. GENERATED FROM PYTHON SOURCE LINES 319-334 .. code-block:: Python from felupe.math import det, dot, tovoigt, transpose s = dot(P(F)[0], transpose(F[0])) / det(F[0]) # s = solid.evaluate.cauchy_stress() fem.save( region, field, filename="result_with_cauchy.vtk", point_data={ "CauchyStressShifted": fem.topoints(s, region), "CauchyStressProjected": fem.project(tovoigt(s), region), }, ) .. GENERATED FROM PYTHON SOURCE LINES 335-341 Instead of writing the simulation results to a VTK-file, point- and cell-data related to a solid body may also be viewed with PyVista. Once again, additional point-data like our stress results at mesh points are passed as a dictionary to the ``point_data`` argument of the :meth:`~felupe.SolidBody.view`-method. Cauchy stresses on solid bodies are also pre-defined in FElupe and may be plotted as cell- or (shifted / projected) point-data. .. GENERATED FROM PYTHON SOURCE LINES 341-352 .. code-block:: Python view = solid.view( point_data={ "Cauchy Stress (Shifted)": fem.topoints(s, region), "Cauchy Stress (Projected)": fem.project(tovoigt(s), region), } ) view.plot("Cauchy Stress (Shifted)", component=0).show() solid.plot("Cauchy Stress", component=0).show() # solid.plot("Cauchy Stress", component=0, project=fem.topoints).show() # solid.plot("Cauchy Stress", component=0, project=fem.project).show() .. tab-set:: .. tab-item:: Static Scene .. image-sg:: /tutorial/images/sphx_glr_extut04_building_blocks_004.png :alt: extut04 building blocks :srcset: /tutorial/images/sphx_glr_extut04_building_blocks_004.png :class: sphx-glr-single-img .. tab-item:: Interactive Scene .. offlineviewer:: /home/docs/checkouts/readthedocs.org/user_builds/felupe/checkouts/stable/docs/tutorial/images/sphx_glr_extut04_building_blocks_004.vtksz .. tab-set:: .. tab-item:: Static Scene .. image-sg:: /tutorial/images/sphx_glr_extut04_building_blocks_005.png :alt: extut04 building blocks :srcset: /tutorial/images/sphx_glr_extut04_building_blocks_005.png :class: sphx-glr-single-img .. tab-item:: Interactive Scene .. offlineviewer:: /home/docs/checkouts/readthedocs.org/user_builds/felupe/checkouts/stable/docs/tutorial/images/sphx_glr_extut04_building_blocks_005.vtksz .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 1.586 seconds) .. _sphx_glr_download_tutorial_extut04_building_blocks.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: extut04_building_blocks.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: extut04_building_blocks.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: extut04_building_blocks.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_