From 4bf94dcc6583e8300cd4c27dd07c2780d4e53537 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Fri, 7 Mar 2025 22:25:58 +0100 Subject: [PATCH 01/40] Fix a typo in `paper.md` (tags) --- paper/paper.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/paper/paper.md b/paper/paper.md index d2859aa6..ccb014f6 100644 --- a/paper/paper.md +++ b/paper/paper.md @@ -4,7 +4,7 @@ tags: - python - finite-elements - hyperelasticity - - computations-mechanics + - computational-mechanics - scientific-computing authors: - name: Andreas Dutzler From f716d08916f431ca17db7de5c9d8dfd2f92a6e59 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Fri, 7 Mar 2025 22:26:24 +0100 Subject: [PATCH 02/40] Update __about__.py --- src/felupe/__about__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/felupe/__about__.py b/src/felupe/__about__.py index 62ad3ee9..d155c125 100644 --- a/src/felupe/__about__.py +++ b/src/felupe/__about__.py @@ -1 +1 @@ -__version__ = "9.2.0" +__version__ = "9.3.0-dev" From c6082017768497366a5e5a14bb94b0929b0f115b Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Mon, 10 Mar 2025 09:50:03 +0100 Subject: [PATCH 03/40] Add new example in docstring of `Element` for automatic differentiation --- src/felupe/element/_base.py | 46 +++++++++++++++++++++++++++++++++++++ 1 file changed, 46 insertions(+) diff --git a/src/felupe/element/_base.py b/src/felupe/element/_base.py index 877e5d8c..03fdbda4 100644 --- a/src/felupe/element/_base.py +++ b/src/felupe/element/_base.py @@ -60,6 +60,52 @@ class Element: >>> element = Hexahedron() >>> quadrature = fem.GaussLegendre(order=1, dim=3) >>> region = fem.Region(mesh, element, quadrature) + + The gradient (and optionally the hessian) may be carried out using automatic + differentiation. + + .. warning:: + It is important to use only differentiable math-functions from + :mod:`tensortrax.math`! + + .. pyvista-plot:: + :force_static: + + >>> import felupe as fem + >>> import tensortrax as tr + >>> import numpy as np + >>> + >>> class Hexahedron(fem.Element): + ... def __init__(self): + ... a = [-1, 1, 1, -1, -1, 1, 1, -1] + ... b = [-1, -1, 1, 1, -1, -1, 1, 1] + ... c = [-1, -1, -1, -1, 1, 1, 1, 1] + ... self.points = np.vstack([a, b, c]).T + ... + ... # additional attributes for plotting, optional + ... self.cells = np.array([[0, 1, 2, 3, 4, 5, 6, 7]]) + ... self.cell_type = "hexahedron" + ... + ... def fun(self, rst): + ... a, b, c = tr.math.array(self.points.T, shape=(3, 8), like=rst) + ... r, s, t = rst + ... ar, bs, ct = 1 + a * r, 1 + b * s, 1 + c * t + ... return (ar * bs * ct) / 8 + ... + ... def function(self, rst): + ... return tr.function(self.fun)(rst) + ... + ... def gradient(self, rst): + ... return tr.jacobian(self.fun)(rst) + ... + ... def hessian(self, rst): + ... return tr.hessian(self.fun)(rst) + >>> + >>> mesh = fem.Cube(n=6) + >>> element = Hexahedron() + >>> quadrature = fem.GaussLegendre(order=1, dim=3) + >>> region = fem.Region(mesh, element, quadrature, hess=True) + """ def view(self, point_data=None, cell_data=None, cell_type=None): From 87e7291432ac568c2cede06fed771b2e0debe0fd Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Sun, 16 Mar 2025 23:01:58 +0100 Subject: [PATCH 04/40] Update the year to 2025 in README.md and docs/index.rst --- README.md | 2 +- docs/index.rst | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index fbf3c830..aceb4d13 100644 --- a/README.md +++ b/README.md @@ -150,7 +150,7 @@ All notable changes to this project will be documented in [this file](CHANGELOG. [[4]](https://doi.org/10.1016/c2009-0-24909-9) O. C. Zienkiewicz, R. L. Taylor and J. Z. Zhu, "*The Finite Element Method: its Basis and Fundamentals*". Elsevier, 2013. [![DOI:10.1016/c2009-0-24909-9](https://zenodo.org/badge/DOI/10.1016/c2009-0-24909-9.svg)](https://doi.org/10.1016/c2009-0-24909-9). # License -FElupe - finite element analysis (C) 2021-2024 Andreas Dutzler, Graz (Austria). +FElupe - finite element analysis (C) 2021-2025 Andreas Dutzler, Graz (Austria). This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. diff --git a/docs/index.rst b/docs/index.rst index 1dbcdcc1..ed5686e0 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -234,7 +234,7 @@ This is a simple benchmark to compare assembly times for linear elasticity and h License ------- -FElupe - Finite Element Analysis (C) 2021-2024 Andreas Dutzler, Graz (Austria). +FElupe - Finite Element Analysis (C) 2021-2025 Andreas Dutzler, Graz (Austria). This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. From 858df419165422bf547c6b879bbda739a3b05690 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Tue, 18 Mar 2025 16:05:02 +0100 Subject: [PATCH 05/40] Add `Mesh.add_points(points)` and `Mesh.clear_points_without_cells()` (#947) * Add `Mesh.add_points(points)` and `Mesh.clear_points_without_cells()` * Update Examples to use the new mesh-methods --- CHANGELOG.md | 4 ++++ examples/ex08_shear.py | 4 ++-- examples/ex19_taylor-hood.py | 3 ++- src/felupe/mesh/_mesh.py | 8 ++++++++ tests/test_mesh.py | 4 ++++ 5 files changed, 20 insertions(+), 3 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index e0931bd1..04492aa2 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,6 +3,10 @@ All notable changes to this project will be documented in this file. The format ## [Unreleased] +### Added +- Add `Mesh.add_points(points)` to update the mesh with additional points. +- Add `Mesh.clear_points_without_cells()` to clear the list of points without cells (useful for center-points of multi-point constraints). + ## [9.2.0] - 2025-03-04 ### Added diff --git a/examples/ex08_shear.py b/examples/ex08_shear.py index 31fdfbb9..774c6813 100644 --- a/examples/ex08_shear.py +++ b/examples/ex08_shear.py @@ -44,8 +44,8 @@ a = min(L / n, H / n) mesh = fem.Rectangle((0, 0), (L, H), n=(round(L / a), round(H / a))) -mesh.update(points=np.vstack((mesh.points, [L / 2, 1.3 * H]))) -mesh.points_without_cells = np.array([], dtype=bool) +mesh.add_points([L / 2, 1.3 * H]) +mesh.clear_points_without_cells() # %% # A numeric quad-region created on the mesh in combination with a vector-valued diff --git a/examples/ex19_taylor-hood.py b/examples/ex19_taylor-hood.py index 9d1577e8..84a2746d 100644 --- a/examples/ex19_taylor-hood.py +++ b/examples/ex19_taylor-hood.py @@ -20,7 +20,8 @@ mesh = fem.Circle(n=6, sections=[0]).triangulate().add_midpoints_edges() mask = np.isclose(mesh.x**2 + mesh.y**2, 1, atol=0.05) mesh.points[mask] /= np.linalg.norm(mesh.points[mask], axis=1).reshape(-1, 1) -mesh.update(points=np.vstack([mesh.points, [0, 1.1]])) +mesh.add_points([0, 1.1]) +mesh.clear_points_without_cells() # %% # Let's create a region for quadratic triangles and a mixed-field container with two diff --git a/src/felupe/mesh/_mesh.py b/src/felupe/mesh/_mesh.py index 2ee1ccc9..0ccb4ee9 100644 --- a/src/felupe/mesh/_mesh.py +++ b/src/felupe/mesh/_mesh.py @@ -323,6 +323,14 @@ def imshow(self, *args, ax=None, **kwargs): return ax + def add_points(self, points): + "Add additional points to the mesh." + self.update(points=np.vstack([self.points, points])) + + def clear_points_without_cells(self): + "Clear the list of points without cells." + self.points_without_cells = self.points_without_cells[:0] + def get_point_ids(self, value, fun=np.isclose, mode=np.all, **kwargs): """Return point ids for points which are close to a given value. diff --git a/tests/test_mesh.py b/tests/test_mesh.py index c2eaab99..b69c5b27 100644 --- a/tests/test_mesh.py +++ b/tests/test_mesh.py @@ -545,6 +545,10 @@ def test_mesh_update(): mesh.update(points=new_points) region.reload(mesh) + # add additional point and clear points without cells + mesh.add_points([2, 0, 0]) + mesh.clear_points_without_cells() + def test_modify_corners(): fem.Rectangle(n=6).modify_corners() From e0860f935217aeb5574721c39c1827995e67919b Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Wed, 26 Mar 2025 21:57:20 +0100 Subject: [PATCH 06/40] Aktualisieren von README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index aceb4d13..16a94ad5 100644 --- a/README.md +++ b/README.md @@ -3,7 +3,7 @@

Finite element analysis for continuum mechanics of solid bodies.

-[![FElupe](https://img.shields.io/badge/%F0%9F%94%8D-FElupe-white)](https://felupe.readthedocs.io/en/stable/) [![PyPI version shields.io](https://img.shields.io/pypi/v/felupe.svg)](https://pypi.python.org/pypi/felupe/) ![PyPI - Python Version](https://img.shields.io/pypi/pyversions/felupe) +[![FElupe](https://img.shields.io/badge/%F0%9F%94%8D-FElupe-white)](https://felupe.readthedocs.io/en/stable/) [![PyPI version shields.io](https://img.shields.io/pypi/v/felupe.svg)](https://pypi.python.org/pypi/felupe/) ![PyPI - Python Version](https://img.shields.io/pypi/pyversions/felupe) [![Conda Version](https://img.shields.io/conda/vn/conda-forge/felupe.svg)](https://anaconda.org/conda-forge/felupe) [![Documentation Status](https://readthedocs.org/projects/felupe/badge/?version=stable)](https://felupe.readthedocs.io/en/stable/?badge=stable) [![License: GPL v3](https://img.shields.io/badge/License-GPLv3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0) [![Project Status: Active – The project has reached a stable, usable state and is being actively developed.](https://www.repostatus.org/badges/latest/active.svg)](https://www.repostatus.org/#active) [![codecov](https://codecov.io/gh/adtzlr/felupe/branch/main/graph/badge.svg?token=J2QP6Y6LVH)](https://codecov.io/gh/adtzlr/felupe) [![DOI](https://zenodo.org/badge/360657894.svg)](https://zenodo.org/badge/latestdoi/360657894) ![Codestyle black](https://img.shields.io/badge/code%20style-black-black) ![PyPI - Downloads](https://img.shields.io/pypi/dm/felupe) [![lite-badge](https://jupyterlite.rtfd.io/en/latest/_static/badge.svg)](https://adtzlr.github.io/felupe-web/lab?path=01_Getting-Started.ipynb) Open In Colab FElupe is a Python 3.9+ 🐍 finite element analysis package 📦 focusing on the formulation and numerical solution of nonlinear problems in continuum mechanics of solid bodies 🔧. This package is intended for scientific research 💻, but is also suitable for running nonlinear simulations 🚂 in general 🏎️. In addition to the transformation of general weak forms into sparse vectors and matrices, FElupe provides an efficient high-level abstraction layer for the simulation of the deformation of solid bodies. From 118e3c80c2ecb17867685d943b358f423c50b737 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Wed, 26 Mar 2025 21:58:24 +0100 Subject: [PATCH 07/40] Aktualisieren von README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 16a94ad5..6a3709a7 100644 --- a/README.md +++ b/README.md @@ -3,7 +3,7 @@

Finite element analysis for continuum mechanics of solid bodies.

-[![FElupe](https://img.shields.io/badge/%F0%9F%94%8D-FElupe-white)](https://felupe.readthedocs.io/en/stable/) [![PyPI version shields.io](https://img.shields.io/pypi/v/felupe.svg)](https://pypi.python.org/pypi/felupe/) ![PyPI - Python Version](https://img.shields.io/pypi/pyversions/felupe) [![Conda Version](https://img.shields.io/conda/vn/conda-forge/felupe.svg)](https://anaconda.org/conda-forge/felupe) +[![FElupe](https://img.shields.io/badge/%F0%9F%94%8D-FElupe-white)](https://felupe.readthedocs.io/en/stable/) [![PyPI version shields.io](https://img.shields.io/pypi/v/felupe.svg)](https://pypi.python.org/pypi/felupe/) [![Conda Version](https://img.shields.io/conda/vn/conda-forge/felupe.svg)](https://anaconda.org/conda-forge/felupe) ![PyPI - Python Version](https://img.shields.io/pypi/pyversions/felupe) [![Documentation Status](https://readthedocs.org/projects/felupe/badge/?version=stable)](https://felupe.readthedocs.io/en/stable/?badge=stable) [![License: GPL v3](https://img.shields.io/badge/License-GPLv3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0) [![Project Status: Active – The project has reached a stable, usable state and is being actively developed.](https://www.repostatus.org/badges/latest/active.svg)](https://www.repostatus.org/#active) [![codecov](https://codecov.io/gh/adtzlr/felupe/branch/main/graph/badge.svg?token=J2QP6Y6LVH)](https://codecov.io/gh/adtzlr/felupe) [![DOI](https://zenodo.org/badge/360657894.svg)](https://zenodo.org/badge/latestdoi/360657894) ![Codestyle black](https://img.shields.io/badge/code%20style-black-black) ![PyPI - Downloads](https://img.shields.io/pypi/dm/felupe) [![lite-badge](https://jupyterlite.rtfd.io/en/latest/_static/badge.svg)](https://adtzlr.github.io/felupe-web/lab?path=01_Getting-Started.ipynb) Open In Colab FElupe is a Python 3.9+ 🐍 finite element analysis package 📦 focusing on the formulation and numerical solution of nonlinear problems in continuum mechanics of solid bodies 🔧. This package is intended for scientific research 💻, but is also suitable for running nonlinear simulations 🚂 in general 🏎️. In addition to the transformation of general weak forms into sparse vectors and matrices, FElupe provides an efficient high-level abstraction layer for the simulation of the deformation of solid bodies. From f1aee2a67784d000fa4a2a72be937bd9661d4c08 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Thu, 27 Mar 2025 10:56:41 +0100 Subject: [PATCH 08/40] Aktualisieren von CHANGELOG.md --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 04492aa2..56b7695b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,7 @@ All notable changes to this project will be documented in this file. The format ### Added - Add `Mesh.add_points(points)` to update the mesh with additional points. - Add `Mesh.clear_points_without_cells()` to clear the list of points without cells (useful for center-points of multi-point constraints). +- Release FElupe on conda-forge, starting with v9.2.0. ## [9.2.0] - 2025-03-04 From 8b27462cedeac76416044614a3a076559c172b78 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Mon, 31 Mar 2025 09:23:57 +0200 Subject: [PATCH 09/40] Docs: Update image in Ex01 --- examples/ex01_beam_sketch.ipe | 403 ++++++++++++++++++++++++++++++++++ examples/ex01_beam_sketch.png | Bin 14045 -> 8059 bytes 2 files changed, 403 insertions(+) create mode 100644 examples/ex01_beam_sketch.ipe diff --git a/examples/ex01_beam_sketch.ipe b/examples/ex01_beam_sketch.ipe new file mode 100644 index 00000000..ffff4ef3 --- /dev/null +++ b/examples/ex01_beam_sketch.ipe @@ -0,0 +1,403 @@ + + + + + + + +0 0 m +-1 0.333 l +-1 -0.333 l +h + + + + +0 0 m +-1 0.333 l +-1 -0.333 l +h + + + + +0 0 m +-1 0.333 l +-0.8 0 l +-1 -0.333 l +h + + + + +0 0 m +-1 0.333 l +-0.8 0 l +-1 -0.333 l +h + + + + +0.6 0 0 0.6 0 0 e +0.4 0 0 0.4 0 0 e + + + + +0.6 0 0 0.6 0 0 e + + + + + +0.5 0 0 0.5 0 0 e + + +0.6 0 0 0.6 0 0 e +0.4 0 0 0.4 0 0 e + + + + + +-0.6 -0.6 m +0.6 -0.6 l +0.6 0.6 l +-0.6 0.6 l +h +-0.4 -0.4 m +0.4 -0.4 l +0.4 0.4 l +-0.4 0.4 l +h + + + + +-0.6 -0.6 m +0.6 -0.6 l +0.6 0.6 l +-0.6 0.6 l +h + + + + + +-0.5 -0.5 m +0.5 -0.5 l +0.5 0.5 l +-0.5 0.5 l +h + + +-0.6 -0.6 m +0.6 -0.6 l +0.6 0.6 l +-0.6 0.6 l +h +-0.4 -0.4 m +0.4 -0.4 l +0.4 0.4 l +-0.4 0.4 l +h + + + + + + +-0.43 -0.57 m +0.57 0.43 l +0.43 0.57 l +-0.57 -0.43 l +h + + +-0.43 0.57 m +0.57 -0.43 l +0.43 -0.57 l +-0.57 0.43 l +h + + + + + +0 0 m +-1 0.333 l +-1 -0.333 l +h + + + + +0 0 m +-1 0.333 l +-0.8 0 l +-1 -0.333 l +h + + + + +0 0 m +-1 0.333 l +-0.8 0 l +-1 -0.333 l +h + + + + +-1 0.333 m +0 0 l +-1 -0.333 l + + + + +0 0 m +-1 0.333 l +-1 -0.333 l +h +-1 0 m +-2 0.333 l +-2 -0.333 l +h + + + + +0 0 m +-1 0.333 l +-1 -0.333 l +h +-1 0 m +-2 0.333 l +-2 -0.333 l +h + + + + +0.5 0 m +-0.5 0.333 l +-0.5 -0.333 l +h + + + + +0.5 0 m +-0.5 0.333 l +-0.5 -0.333 l +h + + + + +0.5 0 m +-0.5 0.333 l +-0.3 0 l +-0.5 -0.333 l +h + + + + +0.5 0 m +-0.5 0.333 l +-0.3 0 l +-0.5 -0.333 l +h + + + + +1 0 m +0 0.333 l +0 -0.333 l +h +0 0 m +-1 0.333 l +-1 -0.333 l +h + + + + +1 0 m +0 0.333 l +0 -0.333 l +h +0 0 m +-1 0.333 l +-1 -0.333 l +h + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +192 704 m +192 672 l +576 672 l +576 704 l +h + + +488 688 m +552 688 l + + +96 704 m +80 696 l + + +96 696 m +80 688 l + + +96 688 m +80 680 l + + +96 680 m +80 672 l + + +96 672 m +80 664 l + + +128 720 m +128 656 l + + +96 712 m +80 704 l + + +96 664 m +80 656 l + + +192 688 m +232 736 l + +EA + +96 664 m +96 600 l + +g +Z +X + +96 640 m +480 640 l + +L + +480 672 m +480 632 l + + +96 648 m +96 608 l + + +336 704 m +336 672 l +368 672 l +368 704 l +h + + +96 704 m +96 672 l +480 672 l +480 704 l +h + + +192 688 m +232 736 l + +A=a^2 + + diff --git a/examples/ex01_beam_sketch.png b/examples/ex01_beam_sketch.png index ee5d8240773d3be40921378cec1d170b556abb47..28f502a5afe9b1b59eadcf4149e9f6ae793a5659 100644 GIT binary patch literal 8059 zcmdUUc~Dc?w)bIBPz-3>4Jwjo`?cJ@g5X6#8Cy|lTLlM5K!(JNfXFO{03jrZ4Xr3> zYct5?&<-R(qRe9$6jX+Y7!n}_kV!LufC-S0Iazz{wSQ~a zYwf+xin| zWqxS?Uf1XxaX;}^D8YMRYIe{1-J26A-~F<^%lE^ws-M;Nejns@&+>a3TBqzpn5M@0 zwq?t_p6pHdVF-3qW<05?|e1)HSJOOvR<6& z%-Iu3%&QyE;sg~1bl9>=0w%-s3&-07q+#SJHGm2JcoJ|TH3+vqtk~I8vhYTL-_oHD zF^V5N(Sj&zNCza#Xn6yEg?k|=b$pPk0ncbPt_0r=msS8BA@kE0ZNS%&%t2l7wRY?O ziy`itMIfG0CnO^Rdv*xTg2q9ss;d2-RY4FALwagmEr{S(4nx7cKd}?hb`JNeAIkKG zuL9I1JJJT@^9AFE@mb0&0g+bY>XO;`grYcIAWISOcF~!SWnTE!Di_W6TR}WE{ZZxTY2Lr5Izk9Cr^hO9Gc6J-?8={9WI>c0+^a}Yz$i`tO%ZOE}yJ&;)=OuRqxx(7n z+4M9?j<>Zwqh2AcislHyXeq5oO$cWun@_?bnS`}{N@q!#Cr@QZ)7y?3o~@WGk@pu~ zlv`xX&X(f&^2X6wW!H8io!ve0a-dD!-Poge+=vGRN&j{Tn;V__w6$E3B}l2LRAJBq z)p1wTf~NvzYdn_Cs_JT)#Nrvg*{+-)YQ#$t;^HOMam5WyC!Dp+J+ZyvDS2tGP(u-`)Vl4Ub_Hsm42bq0{dLHrCPs zeF(DrGz9F-@|ux3Zis)^TWm&E3U70Yz3h9e#5ktBqp&zaWq`=A&8b~juw=Jw*!N7r zK2&pI4FK@7wnZ#Mlu`Dmk4W&MZPDq&APc3T>uBf|-A;HWou4cSb80hC;_2w@geB)PHx3cKEI@pZAk;FATdhVeu3ILVbH_V-R8Y3RQ(aWKf|9@}1ZGl+-cU#3Se$<2PQ z!v4UK8|)gIOK>D88kXT6j%SIzTg&TCN~D%bUqLf!RR>FUC~gRqz-ZQ_?%BV9837>)Iy2~rh=A#| z?W$GATPvISrN8Ep8hQ{7^?=8%xjCKr-LQR2*oY z6a~1UGUXQ_NN21H(Pre(9iBdtE+|4gKeH{Z0*yLXDUBF4nWS;vaqbF?*pFyTZ@yAo zq|)M;dGTkRsLa5*1Jc{0V&(;mvQ0n3VuC?1nI3;naAUq@d{!+(c&~S3ovl2;OF3U| z1W`!8g=;}uStnvZA<8swGy-2+j{p)xSpgymRE8!^!has(X0i@Z-pzREFw{_gLahW2 z0P(Z0o;sre>7x*E9U+J@%X=v67u2*MFJ;Nj;4iR;Ya?0{NaIMhoF-4R(*`-3)kY)R zNmsY$4C)S7c}qG4`B&|<8{V2;vx6Em?Mh@91zI2<;<1j^nB>6E zsS=dkefUE5(ITC<3-JOE7(IlDWFdbl_C zEVI~pXA6U&wefL$RctRe5%yxvNyT3LD{i9s&$EHTMunU2s+GiNE@+jXXfrS$aiTB0 zXtWME`A}?_r_MyWs!4`C>3zIEYRS5hKYDs>56iMG{ARcg?yK|jm}?e49=W{)D_c@; zax+={dXpa1P#$)ZLF#4X!nyGqdYFfkC=53rrR9#G?4 z>ZR-&zP9}4H?^H@jYFmx>lby?M@);}9&8&my_uY_v{KDbG)&L}R)<8CRE_(<$qm(-gu(qgUs!C@D^ixqj!6FWy4a|E2<8oc$4j`jIn{y+W_|Vbf9MWfpD$gz??#OPM;>>m!JX!rQsq= zkcGivH-q#K{E3qSzGEE@4Ek!`=}CA)!~T+1X&@35Q|yye zY>2j)!-`0JO=`2agm$_zBxUyNp8`7e$0d3(llETsQ{j=T=djUXf$>A!Zvc(Zk|b4@~-g7Wg6iB9D z_Ux68<%9$y=?a(u%%T)ocRhDK_5L*KTbiK`FRm_Kza&HKZpd2}#Z6w?4CO&TRqBWs zE?`S87A(=T#-t!FxlAS&>|G&Z(BW-v3MPx*HVS@Vdmx4sA9B*jciyzCZ~2xw5Qo{= zYZpXLXcpam&zq^!Sg@zQWI2TI&ijR$QX0(|jVdlioz~@6pf`+QmD3qPK<$%K=I6fEE&RM zPo<7L`V2Umd;~k&^Sj>D|AL_b=fEn%&>~rO1ORQ_zzkInZMBoPS6nR^X2({nkMWgc zjfT)^_Q$+SqweSC!2aFpyp)ozrDUL3g+8xWw@>!1=5u5%g8SfLq(gP^2}|zVGd!`% zMx+LV4STRxxiG`xfc$gHp+9`9Vnm>9eK+<5+^Qq3ozK#XmchPV4daC)V!E@y$AjdA zf2G~Zd=kr^o^|mUSJ3J8#j-rX*tMwAELdw~2`~&lX8fAn%7%R)cQK3?XaYHvgC@Iq z`;Y==5K|d$PBbS>GNvvbEm5#fW(*9xJ_tvh-B|$V&0#Y(aNQO-;pLp%kp4jRo8up< z(XqAf9k|DSAuwX!Ob3yfv1mymd)%_*rsQHxgD&yb*DT5y&pl!HaPu>)XWiY?RWfnF zh1v{LGAR!4F9`$J!&^NeuLO%27n3ZowOQ?_eL%aEp^+p%5j}MGT8EU?ssWk(j2$~UFEoLGPn37Eznwewp26*&Ey}$&E7ET48Ic_C@Dc?P*rV$mmym!S$)6Z$mb>10>cgIA3=+T zvC9g>xn}olKS`=Ni3-u~aoW%%qnkL`m#0`PoRK(Cw{1{X$TAROsev0OILbJ?Doq}` zr|jQq&L>;^c~&$HLFw@}Fg6G>83f&DQ>%k&Jvb+TA@ z*XOp7UfYZEDfJrWN^nCkH*=y^xVjWfu!EaqS6Zz`#F(NO(r8b(OR!E%TsnvyaGSuH zpu)6;n9uW!_$Wzt+hoV93OBvSK8#KxJ@Cb|WE%PDJEP7Y!9E2Lu8NasHJ15|`u95! zrYZyExWWZP_fjffX7_&#h4UG>oD}Rn;cfHCHy5qgkB@0atw;@Kh<|9qG7C;uDOZRO zf-DG-anq2lbtSTJ$CC@D2!cxuE23O=O<(|S4a_1Rtm3%D)X?aRaYdph(X%4?r$X3_ zmze$iYnGVXdFVv>&=;iH9m~JYGbb3YoJd{^TMxIhRcGn( zLeMwfs}Dx$$+Q3r)=kkTNZzockSQq&6Az& zAFSRm;6c@ObKpiNC~n|aAmMg)5 z1dy%<*L@}JFVy`)I?{e(ap*|T2}z*i;%WeJ3e*!F)|@5P0*z~eD+gVX+qNGHH#t12Mj?G{Gu$U1C{^0@q%Xo)e~iF%l3W=PnWSqG3f!#(=2g82K@#npGMs~6^) zvz$T33J(gva4bh|uuRD+1GWxLo;Mia0=9445KceUfM>N*NJGIdP|RHdX6Fa%WlyqJ zn^l$=7wJx9BZP|s=_Q^Hm0??-=}ZLog7+SU=kHh-$P`IIj=HV{ZWAm2gkALVNT7P* zIBMsa%xfD9>%-^e0j@4od+)lcyA?Xl9~!`{>14JebqgXg1Ya(sT8X8{fRkDyG6XbQ zU+C0$XrT6I&aRx8WIJEKf`tq;i)#Ox@%PA>bL|m6+)l&sJ13qU9>pN7TsfF9dVPiL zy~$qzdZ*|)q1!h5z%911k0MLmiJ)eg;K5XJM1YHRmw3?m`}fO0!wr|$Il&#QKkkk_ z8Xn-%+fL$P;ymzyk{g2CNkRMJmK@+H4i(LuR#X6&JM*8hR!)cxd9IkC6Z z7$nf{XUD}=Tbp;nMa{hUN%|`lU%$NW-sT37D51#n`Ih7sh1S8_YrrpQ*{qvW@9xq8 zul|@qH4vA$VBhjrs^5(3qZCKkWs#Z70F%3~#bWyk7-vmLsTcoWaLZO&wJZ>E!fgT%pXw7UoF~1?b zjXsxaauY-a<&a06mvU+g3uFnXaCj|y{ORHZWhcCorqgIdI1m1hABKfi3-XR^UJyac z!+Mhbmc6BNbltF33knW8ESq%}OlB~RsG)_$1>334Hl0xp7$b0J?)Cx$zu=uE7)1g; z*NZRbwA#|mrC3Y2uM9O5dB0xWtMY8l0$UEk{iyViR0-0keW-x!!(&!(@ngI{&3;66 zpiji-giK%*!m6V%URTKjwlS4iJi%;?GAl;z&bSPS6l(pQSTa6Ra?LQ9Fe$K$%mO(R zdMqes|A6qj;J|2%dIXOAs!*gOa3Pe4J*zEd*>%bZ#Y#6x*QmHNkL3>N=5|_Rakb>* ziW82T$Qb+?HZhrMRsh~Z_vT{)EK57#@cQ~`*t$UW*pn^K8$a{`CLKE9yuWxK?&5y{ D_i+We literal 14045 zcmcJ02UL?u*YKmrvMUPfDx%WF0wPL4dW{X~pacP>1VN=J9YI=xvPx4C=`|=21O;S4 z5FrFaK)Mi`w4l_`f`I@DDgOk^yZgQGcfRkO|NnE&lgu+S_s*R=cW#?W#Cc<#?OXS4 zg&=6Vp6*#w2;!na(594Me*xbIaV0vSVUxG1&S|KyU6=%3wm6Xk zAG*8R4+&mWk)CcBf1G*n=n>Xc*R&@e@yXhw=^7>-WRa2J)#T)ijD_UYc~+15OGHg}N0f*e@rOLz^H;cWy2@%({}wyX=WG#$J=;XX z+k%&%>C7yhXt!Z!e|63M?e9H0G@#vSJT-Sg!}E!W3AB{bYCTtGsOd0U%)y(e0jXhG zn?3NTB|1FF2nD67d9F0dX{Y9>vsdlCT{O6S8hO9?eMBVHrzyal+RlJf&@9G# z{|Ra2r8KfvzGK66{q^<(gqu&JWR}g5f~i@KAzg6$v~l#bMKvY<%jHBu)!2);m2~av zEaNeRT=aH`$H_+UMw*3+sXzzV>utnp8P;a~fdF^OUuwzVe zVhDBQxTFJtEA#QB2r{RXeC)&$&qLreQX7WP=ja+uoloykLXJc}68vG&LQ=4AZP)d?FcxmY&3d9GhF%+I|jD|vb7xzKn zTQ&I0(5k~u%U>S_Kp_c2+kWv2fN!!}5K;G?%kSUPpT%%QGDwNsw`SgdD$mgJT;P$* zE-UO}{l^L7v$!5V3qKTTwq8)Rhn={=VNb!f26|U>u9Ar) z{z6gk{0$u<9C>**V8-vZCWk|v=?Vsq81o2x-K^BSdenOx1x1qLr=4}57e%sa6&^AC zjhTE!xg5T|jlJ&d>OAr86!CdsHY%(SG2pj5qSC31kACjfVXAsnMpt!-?vQ2)0#0#= z)c(?rJWtXl`VG@~kFkd{{MQCGM=T)dZNsb}hi`UN1~bmXVPrDkYba$ z=2W_@0}U$r198xqaXxKCAw|WG50zA0Yvsl}i6f&Nb?`Dp|1Yx(6L^U=|}Qw!qw8As)~=XCHy_Ko6C0ZU$^ zxk8_Rdr+siXB*4-rWYiGohk}aBh_NseNGs|RMEF84#(OGuAgTbH7K!c8kAP4OZlkZ zd&?W++dWQ5icGKO25GoAc3tbD3sRit&{U@LntLP7gUaQD z-ix{V=(?DLI<`B`*d@&Hpi#-k4-883r7q-K(CD>bBhk28G*^d8sPNKKnRx?d@d;!7 z#LkJ<`dXx4lAk#XL;EWc_wtRz-iYYq1S9o}tfb9b zy9{O|9ZJ@F6U;@u7wc}Tp%;se7i0uuSIcfVJg3wpP3BVwTD08P#d<9&YstyWqXGMs zoW*(I$Buv}iLo20wg?{|4db5Jcz8JC$ocTu25EOTawM@eaYjf_o!ypQGXm5~WfXFe zf1}4M*%80tV0SH`h#wLp+cp0{)M!*L4VYA`DFxC-qmAq`diYleeHfgk(j9+uk+s}H zb-z`o&6dKX8nhay@|}5KVbXBDkb#LZ;vmfAJ%{OEsuf8}pAq_#&SWWOMWocxeZ_gS zTKj!@{4)GT>(Y*jj9#RIg77||Od2K+=Sy;d=BM;K>39_d47MWd+|CFtF%?<9Aj^~g zdQ_k=>Y1u}mBlL!U8mUWhU+&#@jlpWpZaeK0C-)sX}0;zvPWEE*Z1S9#8%ds$VXRR znDO93#jwXNuCry8tfqjOaUWJ&G)Y@OWyrN2lg~_aD0yubnDnd+=USeAp$8M{f?+uZ zQoeZHn#jZm3I^ho-;-RkA~f}<0*_8SiYuz`F6+`lh7avP$j*?&9vboKs{LxxND9;0 z8Tg7Ieb2|=q<&KOwR;L?`nUV4@kg2D_YypZME>M&oIE@jAM~o1FeKSAL6|&cYnooJ z#?Rvv?0hRr-|+Q250$E~&C*3ddfQa@Q@Q+~aGk(N+BArG2A~ZEFRs5Q6`JfaykN!8 zFq_9QLN7Uw-NDf5@CQbADL`VhzTPBH&JOyH=7zE=i+?jjF!#F??vP#8lkQMQnvh6x zLACPzgRZKJgj4L;{w>#?3(5m5m^B>&7H7-arkA#-n_pa`-`UuK9bsC_Prp=^!l;eS zZc^#^*(M(L^eE}svMKsnu=9kU!)+%X6f^nz4*knYe~KKS7B~uIRZA?@if^{vYt-*2 zENSizv3u^TfwHeTjY(t&^8`1VqF>G>bkQm^DpVGi(h=$HOX-IRV6T6}G~h8hg|4%c zy^Y=mE3T>v(kQryHeTpL&$?-%w9aVu4Eg@I8~K0Q|J`EJl{x++UDy7DHfVbtloFjx zXI#1N11XQwT2NW9jWnx$e$M05f}Of#5(;4UE<0z0BtED&nORxQ$_bScSRLJgClru4ZKi@2ewZ06!BQKfD3| zi3?6$n=P)m_RnfkB(krngtC2Sj2Dn_dCzGdg0CjaX-5)Vow_E+y(>E({-Dk6Q$b_2 zRAy!g4{Er6qyT2fvz{ojCY#YNww6ahe+|R+IQS3V;U6O{-yv$$Y_y23bu63Pg0%v2 zQ*Cw_t8$2R=N95-*J;|CWSYR1 z#6&jq5vLdU)J$x+l^Z7aK8~a4yy$PeQT*J{q~6g5;&Va9F+qmHD!;?8P;035s0k9i zva;vVTa=9!)cSO&GkPogx|SdCaO5LFMGWYGQCVf~al5d%$&N2ahP}W*p?qk>+r(R| zIMsx2LxOKRFB*TEdm;*tr>c{ToSE1vlWl^hgTme8k~W7d^cUVHn|9}J55DLYy-5s5dzIB)!H}Aiz?_g*r%Vq0Yj8qC)+z(xiQbPr$5F*SEICC^7t7zlV8?#d z7K3*G`Ha72NF`u-+Gpn0InDJ^ zkG48%fY0h!;j|qVpes*`O~7CIfHuxS*c+^Y9D)Hd9kc0+$X$T2s#i_UDLDy(`mBjmQM<(bWqb!iv2 zU6IK-phCZZ+EzYZ(ZwjunZ=JUaXANB$+Y@wVI`ix{X7g@puM~ek_0$xW)a^ zq96=>=Xfxk&mOlAVC`fk4p_D{B~@S55@-pXUeU4fSHm4_$Q`4Xvy5VO6Jya;801N7 zk+C5X+OhQ4rF^qRZ~NqwTDwJN92#~X`nbJg*?7eI=*mTVS25YLs%o53d zpOw%fsb0G*WrHb|zu={EEN_YZ1vu$`=zc5B-H1V7`cwx7Gsu`JV$EFRZ|_1KIB`k8 zLyDkJqG#LWWmsgjdOS9L3Qy&aa}02=RKm)1YGmfeEqw6T#eM8siZ2Svd4mt1)6AKQ zbI$Pz?kY6vcyZyn%8_r{Z1AXvnQD%@mnbE$1A&e2qUu%tm3z9rjI$9*_udV!v(Ob~rId)EMUBBq8hvup#T*}<$nXcT zKdIiuWID6hQDW37%r1=`h)OHqPM_4W1F_|ZDB@rnSH3>!(XJi91O361BlNJVI6q3C z4v)S&6(@Hq`9C7vExo}k)3aC8ep2>qy1Z3|OlmBIPWH`FHJO&bKyU5fpyC8H*8~hQ zWx6w~Zco?l8@QfIeQvm<@1_vV3(HoAZZXNC0xKgk=~#NT|h3LYNlnppyZcN81R!IN<1$^=q1A2^+Rb_Hj8ckS%k z%PqD7bdZ{)Aj5ji;nxRxpnO>K%5Zi|Skn=@$Y(u1r3nXdO;63=>tYmEKT;@^uX{ER zo?li!-5J)ovp1b)dE~I~CZMVO{*+XNRt!4HK+lEXKHNGeK;y6XINM1M51NVy51*44 zG~LB4YwPfxH#;W}0%)@z~*aF&~{@-qa7%#vcUQZbe#v>0U44 zik;@FxeMy;P932`ulE0&LWOt%a0vI)|5Z0tr9Imgn(C{Q0n33YqKKIRlXA-J9{m1K zD^hdnIobuE!Rc^8W(C;j&OMuBdhmR(dk5ZKtnN7Vkr*`n#`KDbH)e?iCJ{3ky|7faAPP;k zKCc0mRab_hT7BtTfbzGAgxotH{PmIhC$3BcaEg*S8KIkr{RKygrj{B##E+~&KIVCE za!<3?jJfedGEsCY@CbJI-Mbt^`7V7c11Be}=dG+bHIr{&G|F{SXZ%aL+R6?4H!rNm zmsl4T1@jlB01DGx3k6@|(PkL6bsY!GIP?(4L?8&74mYbJw`*6%q}tO?II|_~y=S@5 za@N0Y)q>LFcztnMMefQIf}1NoG@b{HDg3~#&fWvRuim_#@J9>psal@QSbc8wHk*~E z_3k!ey4NR38bIn)G4x59cZDOvuC`m^CU=3VO6!;926HD&Q(7qNOn=l+`IgB$cljaS zlK?Gu<|tRq=JP;6Hi5~$f(@w&sxP`SD|ypweM7mY*At5>c`_AEi;*T32k+0-9k!`ns1{J_LW8ibsQ!05x$VD% z76Eps%U#`?FN5?$ubDnYt4InmlQvsw@JrWN)BCwuos`YY| zi&t7;oj;cmyXV9;6q&U+VK0`;W4pLS&6(A-aWxsE6S&hJv;?1NFO?BYa>)CLSDqm4 zc{ilT{k1GEx^cy3Y*#~m$U;JYu$4b1h$)vg!24+O+Iwn{nlXk8JuRzoac?_vY4n zuXB2$wz5XDPqXvl*y8tP6TLpSaKF+wjJ~N1=Oq#4j?PJ05$6QQ475DTEl3Q%aWq;@ zL}1nDGA^99=){wx{0;a{92;=>*Oq!i?Yp9+*DCPmSKPlc00cOF|&-O4i2E@t`o>1&Pd%bP=bZ#`!Nbvm+eOz+w3tV&$;0ObbQPgzxho`VuFF?qA-E=v==;$N?xCL?%46^1}BcZ z;G%M~Vu^TcVM2@~y>{S)GsoRAB&p>DM|jm;(dTJBc*EcYk~lrQcGrgmDEiB47099Grl&$ z`~?E$`imc~R!wQv2Wwv?OKIjB%ID0+`D)j~&*r^C`-04C9!O;uMdrn5A6ui^BPUvu z<^S3BFudqo^5;wqqE#?L>!i%Ff}mqqnd0hFO4ej;oO}^C{?jl>+9J=cWB|8G0rA$@ z*@EKj61SJ{^vnqS3IC(1Kxy}gyq4@j^wX^@45JBSZ9WZB@PaY8+#fTDXQb`Bc0csh zAbHd=%6!qA61dW4&nB6m-@YmO?r~y$(z^QCZ!F^ov;I+_pVSzD*LE|1~->Y{~= z0zXz5bdxXw&SZ*|O>U|5Ru&fc_g7iy7-P0GC%1rYeA4dix5Rcx3DjuQ+*LOpkK%_v z?3?8=75F#E{85X&<%_w@_H3qm_|n6zK_j-dwpv;sqLUkNR?o{KF03y49bfajH;DC@ zRBlwysc`qMLZqxK{l721wh1s!zFt;2;_^QZFw`2#uXVdAATBv(vYI!8dQj(4l?JPC z?uYX zJ*Ke(E~2s%5v#}3&eg`17kL+UAQt;WZUnf|Y||LLGP(waF=xXGD_L3fo|;5;T0_vs zm&Xf+Qvd4JNb=ehv?E!8bD{vkiZinB z--MF?)Y0i^Y$5A~T$Q2P?F?5@L)ER-;f8<5z?H!NU>Kzc!|gc_Sd#xYmbi-}s}3F@ zaKfGM5{hcc=5ib)*5Ly@Zu*1)kJ|9b5k+CF0oWqc(6gU2w z%!-?4##?q2Kt#wME!SOOcwa#dPIinPU4snuYuzKvK8oq3Gy9YSl&kj(1i`xi;$4h0 z32*6xUUK-4l@m<4_$W3%a5*CH^4Da8rZ_cv$D_BjZ|X+@IE1nnx3a`%)U)q?%ZkvF20Y~7@+MYctedtuM!z|7Hd>j1C76G~S;Pe{H*;n7Ch?r5H0v@)&@GphKQnDI&@qV~41^l?FM zw6ph?C}|PiRVB=;uAwY{t@|~6yBGWSaDRR;$ILV&npVeC675~i*1D^Lg~5H`h6}Eyaq==&GJ=emN5G;a z%g)ji^m$!(iLg4{9zk}FR6e-mL(Kf@I{M;wsOMWz7pENtqL&9{wz7DYbf8xyN7nl4 zJ)*u1L;%xE1IwWRRI4C2dkS=r-OG)i0OEl?^$mqHm@I(pa0dgN*otZa9=hR~d`V*D z#!S9@#Nl_@@7i`ivWY(&;$tC4)?fp|n;>t!pKkL9Zg&2Q!T`^0n?n4v|LI4ASf^gv z*7MMPs7Ax|r+4jDi>EY8BmX!9*@=Q-zq{KNqhfGWkm`T!iY6H#2cpFx+324%K-?g| zc2Z`o@sWto``@7?9l+*eKPo}4rG6uf1MBg97J$HwC>z9foJ!S%NX9w&yfw_jYk+|4 z&FtFiwb6$my9|Y`EMk;1C}bbN_e?bI0%HPCg~^Eton-e6?M(e!#PqEwi%CdKjJ81L zS#Dck!2j)ISyIj4)005m#P2+@`W&74~!$|`yfA+f2A-HVcJM=4-ITpUvNVS8YKKMb{n-F^d3 zn_(;jzeb8S{sG`V?ygj%^{B61jKvcwnV_OZ+$G4GO(OJ-tgdKa@Zqq1>m5C-pGbZy z97Jc7=qjuTZ_LZ!Lu?_cAJsbx*#?#Y;Q56%OfR_4qL+}PNYChk3;u!}dM{s}B@gEk zl8!dRw>?)}Q)kJ7N=Ri3ZH$)ZDaLGDzfgw^3_~qttvhyP)5ud)h8zoo=l`ou<}G8S zpIcTONDe@DyE!Mwh~|bZ72-1XCS5eqd>IEI5h8@vWYbX0&K7{NRlk6Q(!fS`(Ubeh zok{x!IS#H*wIt=%`1@O0yITI%^R&ww^2-zkSlB>A!JOc8WP0vWT0~xPb_-6T(4o6~ zitr{)1pWM^U;I9h9o!f{->=89$Adqt&hPz%JJT-Hh2m`=sfus`3Cq2Jan3LHEe=1~ zw;X};(^6icF=zqyE{EFS=PMZzsB(<2wBnj-F%Y4e(wtucXGN$C^(8q?)8O}7isL5l zyEC$vhEqE*FiC8gSsQvBn{@qA+>e1%yu57`2l7Z6Gv!0MErF$P=Z2rG;N;sIgr$Dq zRLGeOm_4|Z<2NXIAvhFey--H}*oQXSKabCaea#Zm-;Odq#sxvL$6vcwb{MU2$=SrE zQv3ZH;eLs|?%_;d=F?LHxvX*;b1zQS?B|5^9x|s_Cl1!bstuNubC_jXzC;G4>TM~p z7xP(YDxKD#PGdb75kC5Z3KVsnIrw5S0M^uovVMQq>kPh)RS%mrxIy(>Or82#N%PL` z9-VrB0<6?Wft+7Jf{Y_@g-5=_Lq6JhwH*=HKy#+_50W993*YI1rXzo>vv4W!9~>;V6A$4VHc29@Lw$-$*XL-o$T7>sz5g8y~3QqkHN6coZNk01)HaH z2pymz%vvhmvTpu^gHp7L?jd%%(*k=nuDo8u(C7!)inQV?*~njQA_BBq$!$ki`~cC* z!J-r=Fynq?0x7JkyW?hI&CQ>ZP-Kc=!?J@UQr0`b0dQZ#ABNg4Lmf&&wBNjxWcUCt zpLF-La{(oaimPlF!@BMX30!9>`S89V&^8Yy-e=0vP+aZEse)mzV8NNudgo~_yaUDm z%N$^O_x7+lA%r(OdZ}_%N~~fg|Anovfye%puIU?#V4=1?0+(Vgy3D`@+L&l`y?{lP zR>Ek_L&En~K-wLr1u+Afklgk#m)vv6P=4G9IOLd#%T)@%J@lb_7;-_co`cwQ9whnJ z^NQJe&F+l~rRVZG`$p@Y*P_q3EzCw7hM?U>E(>fg@LT}Q=8>wiyxW_nOy3u@M{&(> z9IKU%(wL}A34%hnsZSEC#;@g7;9<(Eq~$2m@_oy@<_6nQtfc^!p(t;W%M5`;qN_~0 z3C$~nb-4zfX=X|h^M*&r?}tRBkRQhSnICh|Eu^@x7w8*|Kr!PA>E4Xy++ua8W+Vs8 zVgfd*Pq^1BKq4_I@7A-fXzWFJ6UQyeRnsT$vu2V;o=jGsuVjivlg`Ajwo4N5Pw9{T zth5g>uq|R06n9)J=Yi}jDxR9q%oWxg?e~T7>23eq=f(Bu<~G7w*=54o=cEiv+;79~ zU3kYYFs{~e`5c!dax|%B{Ph48CywQ}e!D>Dw~Ayx(m)P64i0XC&d0$cF)4IL3}a{r z(52ct)G zU)!L6haivS?rwerC6P2S$wC8v#{CgMbVxBYF{moA+^{RN>_V3F$O@03^ZeNBK|ED@ zTVQ+<>@a)(VYHWBM`CPyIrYrlq_>6xQ!_4@Qs%v+U!eX=E^~laFw=IR9PuIq$j8uS z%M169L~6^pOWRm=Dr4p^Y>XaTlN3~m>NuG~bXm;xtXy>{N9&|?Aor=zdqVr&-SGQ& z*cEclkme4Xthg;pwkBD+%&@OFJ4n$?zAzi*wx8t5k?`mq#}bOGk5{{hdfVp{y5_52 zyw!AGe*z;IcxS=N4Wz0EF=b2T^J3HGy;EqDd14o%{Z9tDiJ1F}dM3g$-_3y45R9zh z`N*Af4wdUfNq+C4!pQ0FMr)<4vu7fu3l~ps%FHq?NhE%CWktMk#q=#*ZtfNAsP~`z z!aLkG-%syl%)226MqCOVT%Sb`k}ay+>G&*tCNZ|G9OZGrYIPzIviYm5G@!V<)m_j? zE6}s&q6JPxk%|{=V;+=JuIHF!0=CJ~9P9-5Y%D5Rc8>;#gWCWfAv<#gnxPdC{mU+M z><(bV)o|(@>z))Jg|mt-qMDPBDpJ(7V>^6Jz1c70Uo_hL8@N!= zT{p-5A|)_$z`v73AlOb|+xii@d6mKCF(@gLIx?g0G6zph{n}jxk55@A+K`7S`Ia(W z<0XQ5dFWfs==9X}IQ9bI{aT^}xGEOH4B&(*-Tjz?lto?xX|ceKdqVyeDM1bKDRg4Q zBuj{Mq%5}`7*mw>bV0R@;)K8PTM$HtBdfL#LRCCt3!f5i3m#NDw+*)&7mCkB=dL%DZNTnagJX>XwB1$ zI%?^*b{2g2j^*?nve6^0gnzgS(!(8R4Y)KoQ85UNOw70?U5V?Q>na4_pmB5+gX9qH zo(4zTPAWLIr|s3e!}0OSO=-1fSwmh)y#Q!LD(5UtL)!t?R(kTP!m+y zA2=1&7a)y)4Nm&njvV0!0!kEUQP-d99PjY=TCzpk4@u#d4e(}F+c?A+{)?DyWYIYx zt|e#Ac^)E`N@2K7M6=hTVg9DVlIUd}v{}`$|1zm= zd-YdF7_CTjYd!Q&eqW6Y%H8zLG87I`2z|@ngn;mC}EdVvZbYlMR$C%nkll# zdbp{qq^y)aW_oHs28@FR+dBOT{?8(x4YnEoozfC8LihX>dHe>YGQCkW@aEW<7<;Y3 z#&1}Wp1|~;YIoIOy;(Oqdv+kb>&BD9Y@fihSKT|u5sS}MY^(>F=?ot=`clByxmehN z0|$8i(T4zT0B5vDhxyP=5F` zTjbk}K?nas(X+E8qL)O1M0p=R7PE4xDDtC8PB&U~>XmM?bs)u;nn&#OjX+%2V z-C%=J$Z2Wa>lfS=CnN9p|F+Nr5W3fKTs$+#FDP7kl zA5>~N>nA45FVO$l7#+(I6=L{Y^5NOb2cyYTx_Vu@G#;EvG7A`IjsNB@C zV~F^|vHJs5%{U>3hG0uAqw4VRaA6<)Z)>*$WWrsZd7azNC;YhkB~fTQW31w7G|6C` zR+@po4TGOsB9q8$-ZlM<*@BUTeRP=H?xc9KEhkehwqhP#Pf1YRKv=5rQ)iC+(H1L# zsV^C^vH{VBqxUJjwsiz>pNBUx=N_zYoBP5UVd0c#7qlbM;(2gBDoCE*G0qt0RhooilI8sx^x9OAw*uy7gf67>_s34!TIO7? z*gSf%9j601XFh$2=30;UQQ$2`G;rx*3dXs%(UPb)3p{fu^R)u)SXva$uljaDImPOR zfi&^Kb;jq{BQuru{yG~y+CYy7oco(551CmNG~+h9MPJ|x^240soF%PGfLlK&H8~PM zPKd4eZ2XWJz(NGbX~y`d0u(EUaQ$jKVP)s^#sU*T{s-LUFpLL`u)E(C^)hv4JIau& z#?Xnf3A%m+*B=04b%a~EGnDk3krq_*V0oQ+BmWn;@;dM#Nj+2hcJz>QfeYFOQ8b^Cn~P!Ifw^CG%(TSBabGCsT( zPuCJfXfsi4U-++(-DL%uEca^(khc6KMB7|(K{fB(&|U&yPZqo=apkKg3c{WS@R?0K znKC{18ppWJzang%S@(zfE!M}yw!Z%Znga~Lo{Z%`nJ`m}yEvxh*5#iOUcyyV?M%6{ zl|$c2Atz9BF=BPxYzI;!tYM^Oc!whiN%7q`g2bu9tHP>Y;}dePkN>T~JxFTb>V>W~ zD=9I3nn7_(LcX$I=XhOqK9`g$Yi^R0p76;=jXeRTrnI_e*#f5*DIq3tQT;K;^#Qv~ zzN0aZQ=v$6Ca80y%Pa?#1=MQ_T4~}YxOV@d_y%*kE*LPscN}l$&2eeffZgiuR{lYD z&C3T%1$2QQG&ynoqjajEiCxjQ7dU~TxlLhGA2~ATn)#2?b43~)uL|lt7CB7jNE4K_ z^ZU#&(7_vb5~fYI*)|?{nJDMlvC_bY3^h& z&SlQ;Lo_0bouVQHYCoRKyMvdR!1U~+UuXG%w1)m1NH0L8mt4lnZ!Du9Q+mPe1;D7g z-L1+(Og3x+%Sn+zZ?jj_=_no!58ldfVs=pYA3hb9uV+VH)S}ABOIFOUl_a z$^2jntC({Owp7n$RLrPt{fykuUmT;~n=b|p32kIiEyn$pZ?rLYR(uArFyIldDFWT3 zv$PIU3>%<6+AcF<`6=sFAS4C+0o`X`OlpnPqO3`Or1BmVre0jU1|84XzuV8d$@AjQ`-ViBiOPuuux;c2_O{|8$n B+SLF6 From 8574f6a64a91549a651b59cf5d575cc8c3bfd56a Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Tue, 1 Apr 2025 10:41:37 +0200 Subject: [PATCH 10/40] Docs: Add howto/project --- docs/howto.rst | 1 + docs/howto/project.rst | 50 ++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 51 insertions(+) create mode 100644 docs/howto/project.rst diff --git a/docs/howto.rst b/docs/howto.rst index a3d5609c..0c4953e4 100644 --- a/docs/howto.rst +++ b/docs/howto.rst @@ -16,3 +16,4 @@ How-To Guides howto/umat howto/umat_hyperelasticity howto/items + howto/project diff --git a/docs/howto/project.rst b/docs/howto/project.rst new file mode 100644 index 00000000..2084d4f6 --- /dev/null +++ b/docs/howto/project.rst @@ -0,0 +1,50 @@ +Project cell values to mesh-points +---------------------------------- + +This section demonstrates how to move cell-values, located at the quadrature points +of cells, to mesh-points. The results of :func:`~felupe.project`, +:func:`~felupe.topoints` and :func:`~felupe.tools.extrapolate` are compared for the +Cauchy stresses of a rectangular block under compression. + +.. pyvista-plot:: + :context: + + import felupe as fem + + region = fem.RegionQuad(mesh=fem.Rectangle(b=(2, 1), n=(11, 6))) + field = fem.FieldContainer([fem.FieldPlaneStrain(region, dim=2)]) + + boundaries = fem.dof.uniaxial(field, clamped=True, move=-0.3, axis=1)[0] + solid = fem.SolidBody(umat=fem.NeoHooke(mu=1, bulk=5), field=field) + + job = fem.Job(steps=[fem.Step(items=[solid], boundaries=boundaries)]).evaluate() + +Cell-based results, like Cauchy stresses, are not projected to mesh-points by default. +Different methods may be used to *move* the cell-data to the mesh-points. + +.. pyvista-plot:: + :context: + :force_static: + + import pyvista as pv + + plotter = pv.Plotter(shape=(2, 2)) + kwargs = dict(name="Cauchy Stress", component=1, plotter=plotter) + + plotter.subplot(0, 0) + kwargs_sbar = dict(interactive=False, title="Cauchy Stress YY (None)") + solid.plot(project=None, **kwargs, scalar_bar_args=kwargs_sbar) + + plotter.subplot(0, 1) + kwargs_sbar = dict(interactive=False, title="Cauchy Stress YY (topoints)") + solid.plot(project=fem.topoints, **kwargs, scalar_bar_args=kwargs_sbar) + + plotter.subplot(1, 0) + kwargs_sbar = dict(interactive=False, title="Cauchy Stress YY (project)") + solid.plot(project=fem.project, **kwargs, scalar_bar_args=kwargs_sbar) + + plotter.subplot(1, 1) + kwargs_sbar = dict(interactive=False, title="Cauchy Stress YY (extrapolate)") + solid.plot(project=fem.tools.extrapolate, **kwargs, scalar_bar_args=kwargs_sbar) + + plotter.show() \ No newline at end of file From 1276a029e90400a6347e02916c54047018f2043f Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Mon, 28 Apr 2025 21:21:27 +0200 Subject: [PATCH 11/40] Update license in pyproject.toml (#954) * Update license in pyproject.toml * Aktualisieren von CHANGELOG.md --- CHANGELOG.md | 3 +++ pyproject.toml | 4 ++-- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 56b7695b..1caba695 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,9 @@ All notable changes to this project will be documented in this file. The format - Add `Mesh.clear_points_without_cells()` to clear the list of points without cells (useful for center-points of multi-point constraints). - Release FElupe on conda-forge, starting with v9.2.0. +### Fixed +- Fix the declaration of the (spdx identifier) license and license-file in `pyproject.toml`. + ## [9.2.0] - 2025-03-04 ### Added diff --git a/pyproject.toml b/pyproject.toml index 112f393e..783050df 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -15,7 +15,8 @@ authors = [ ] requires-python = ">=3.9" -license = {file = "LICENSE"} +license = "GPL-3.0-or-later" +license-files = ["LICEN[CS]E.*"] keywords = [ "computational-mechanics", "fea", @@ -34,7 +35,6 @@ classifiers = [ "Development Status :: 5 - Production/Stable", "Programming Language :: Python", "Intended Audience :: Science/Research", - "License :: OSI Approved :: GNU General Public License v3 or later (GPLv3+)", "Operating System :: OS Independent", "Programming Language :: Python", "Programming Language :: Python :: 3", From c6fdbb8d2993d8d3cf5de4d3994ff59355319397 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Mon, 28 Apr 2025 21:37:51 +0200 Subject: [PATCH 12/40] Aktualisieren von pyproject.toml (#955) Fix license-files entry --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 783050df..5240c151 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -16,7 +16,7 @@ authors = [ ] requires-python = ">=3.9" license = "GPL-3.0-or-later" -license-files = ["LICEN[CS]E.*"] +license-files = ["LICENSE"] keywords = [ "computational-mechanics", "fea", From b6ccf7bfa9876aa45e719710189d7ba9f45ed514 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Tue, 29 Apr 2025 09:45:39 +0200 Subject: [PATCH 13/40] `pyproject.toml` ensure PEP639 support (#956) * Change `pyproject.toml` for PEP639 support see https://peps.python.org/pep-0639/, https://packaging.python.org/en/latest/guides/writing-pyproject-toml/ and https://packaging.python.org/en/latest/specifications/well-known-project-urls/#well-known-labels * Update CHANGELOG.md --- CHANGELOG.md | 4 ++++ pyproject.toml | 10 +++++----- 2 files changed, 9 insertions(+), 5 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 1caba695..15a141c5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,10 @@ All notable changes to this project will be documented in this file. The format - Add `Mesh.clear_points_without_cells()` to clear the list of points without cells (useful for center-points of multi-point constraints). - Release FElupe on conda-forge, starting with v9.2.0. +### Changed +- Change the required setuptools-version in the build-system table of `pyproject.toml` to match PEP639 (setuptools>=77.0.3). +- Change the labels to well-known labels for the URLs in `pyproject.toml`. + ### Fixed - Fix the declaration of the (spdx identifier) license and license-file in `pyproject.toml`. diff --git a/pyproject.toml b/pyproject.toml index 5240c151..5e305fd3 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,5 +1,5 @@ [build-system] -requires = ["setuptools>=61"] +requires = ["setuptools>=77.0.3"] build-backend = "setuptools.build_meta" [tool.isort] @@ -10,9 +10,7 @@ name = "felupe" description = "Finite Element Analysis" readme = "README.md" authors = [ - {name = "Andreas Dutzler"}, - {email = "a.dutzler@gmail.com"}, - + {name = "Andreas Dutzler", email = "a.dutzler@gmail.com"}, ] requires-python = ">=3.9" license = "GPL-3.0-or-later" @@ -86,5 +84,7 @@ version = {attr = "felupe.__about__.__version__"} [project.urls] Homepage = "https://felupe.readthedocs.io/en/latest" -Code = "https://github.com/adtzlr/felupe" +Documentation = "https://felupe.readthedocs.io/en/latest" +Repository = "https://github.com/adtzlr/felupe" Issues = "https://github.com/adtzlr/felupe/issues" +Changelog = "https://github.com/adtzlr/felupe/blob/main/CHANGELOG.md" From 8bf09ed5f1ffbc7d005cad622c31f41ff302131e Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Wed, 14 May 2025 13:35:13 +0200 Subject: [PATCH 14/40] Fix missing import of `TriQuadraticHexahedron` (#958) * Add missing import of `TriQuadraticHexahedron` to top-level namespace * apply isort --- src/felupe/__init__.py | 1 + tests/test_mesh.py | 3 ++- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/src/felupe/__init__.py b/src/felupe/__init__.py index 58108f9c..f481bc25 100644 --- a/src/felupe/__init__.py +++ b/src/felupe/__init__.py @@ -58,6 +58,7 @@ TetraMINI, Triangle, TriangleMINI, + TriQuadraticHexahedron, Vertex, ) from .field import ( diff --git a/tests/test_mesh.py b/tests/test_mesh.py index b69c5b27..805fd713 100644 --- a/tests/test_mesh.py +++ b/tests/test_mesh.py @@ -571,9 +571,10 @@ def test_expand(): def test_interpolate_line(): - import felupe as fem from scipy.interpolate import CubicSpline + import felupe as fem + mesh = fem.mesh.Line(n=5).expand(n=1).expand(n=1) t = mesh.x.copy() mesh.points[:, 0] = np.sin(np.pi / 2 * t) From b59ca3b3333db8fc494d987ca7da00c975bfeb37 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Wed, 14 May 2025 13:54:27 +0200 Subject: [PATCH 15/40] Fix the path to the logo in `docs/conf.py` (#959) --- CHANGELOG.md | 2 ++ docs/conf.py | 4 ++-- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 15a141c5..3612e21b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -14,6 +14,8 @@ All notable changes to this project will be documented in this file. The format ### Fixed - Fix the declaration of the (spdx identifier) license and license-file in `pyproject.toml`. +- Fix missing import of `TriQuadraticHexahedron` in the top-level namespace. +- Fix the path to `docs/_static/logo_without_text.svg` in `docs/conf.py`. ## [9.2.0] - 2025-03-04 diff --git a/docs/conf.py b/docs/conf.py index 4626c7d7..f6628a85 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -143,8 +143,8 @@ ], "logo": { "text": "FElupe", - "image_light": "logo_without_text.svg", - "image_dark": "logo_without_text.svg", + "image_light": "_static/logo_without_text.svg", + "image_dark": "_static/logo_without_text.svg", }, "use_edit_page_button": True, } From 6709afc19c300b635911e16e841876fb9065d3c5 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Wed, 14 May 2025 14:11:03 +0200 Subject: [PATCH 16/40] Update _multipoint.py --- src/felupe/mechanics/_multipoint.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/felupe/mechanics/_multipoint.py b/src/felupe/mechanics/_multipoint.py index 867e3c1a..c3294247 100644 --- a/src/felupe/mechanics/_multipoint.py +++ b/src/felupe/mechanics/_multipoint.py @@ -147,9 +147,9 @@ class MultiPointConstraint: ... point_size=16, ... color="green", ... ) - >>> field.plot( - ... "Displacement", component=None, plotter=mpc.plot(plotter=plotter) - ... ).show() + >>> # field.plot( + ... # "Displacement", component=None, plotter=mpc.plot(plotter=plotter) + ... # ).show() See Also -------- From 375cc9318d2fabbe389dc97439f3ce17718a0908 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Wed, 14 May 2025 14:17:56 +0200 Subject: [PATCH 17/40] Update _multipoint.py --- src/felupe/mechanics/_multipoint.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/felupe/mechanics/_multipoint.py b/src/felupe/mechanics/_multipoint.py index c3294247..229332b7 100644 --- a/src/felupe/mechanics/_multipoint.py +++ b/src/felupe/mechanics/_multipoint.py @@ -353,9 +353,9 @@ class MultiPointContact: ... point_size=16, ... color="green", ... ) - >>> field.plot( - ... "Displacement", component=None, plotter=contact.plot(plotter=plotter) - ... ).show() + >>> # field.plot( + ... # "Displacement", component=None, plotter=contact.plot(plotter=plotter) + ... # ).show() See Also -------- From f4c806975b7ad5a5fa2f7977f1731f64e2f471cf Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Wed, 14 May 2025 14:29:23 +0200 Subject: [PATCH 18/40] Fix a typo in the docstring of `MeshContainer.from_unstructured_grid()` (#960) --- CHANGELOG.md | 1 + src/felupe/mesh/_container.py | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 3612e21b..ca7f8733 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -16,6 +16,7 @@ All notable changes to this project will be documented in this file. The format - Fix the declaration of the (spdx identifier) license and license-file in `pyproject.toml`. - Fix missing import of `TriQuadraticHexahedron` in the top-level namespace. - Fix the path to `docs/_static/logo_without_text.svg` in `docs/conf.py`. +- Fix a typo in the docstring of `MeshContainer.from_unstructured_grid()`. ## [9.2.0] - 2025-03-04 diff --git a/src/felupe/mesh/_container.py b/src/felupe/mesh/_container.py index 8cf234ff..f8087860 100644 --- a/src/felupe/mesh/_container.py +++ b/src/felupe/mesh/_container.py @@ -233,7 +233,7 @@ def from_unstructured_grid(cls, grid, dim=None, **kwargs): >>> import pyvista as pv >>> >>> grid = pv.UnstructuredGrid(pv.examples.hexbeamfile) - >>> container = fem.MechContainer.from_unstructured_grid(grid) + >>> container = fem.MeshContainer.from_unstructured_grid(grid) >>> >>> container From e7e6de060273f5f8757b4de7c64efe308c9a6a0e Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Wed, 14 May 2025 14:39:25 +0200 Subject: [PATCH 19/40] Try to use one-line statements in the docstring examples --- src/felupe/mechanics/_multipoint.py | 8 ++------ src/felupe/region/_boundary.py | 25 +++++++++++++++++-------- 2 files changed, 19 insertions(+), 14 deletions(-) diff --git a/src/felupe/mechanics/_multipoint.py b/src/felupe/mechanics/_multipoint.py index 229332b7..6a4c8a9a 100644 --- a/src/felupe/mechanics/_multipoint.py +++ b/src/felupe/mechanics/_multipoint.py @@ -147,9 +147,7 @@ class MultiPointConstraint: ... point_size=16, ... color="green", ... ) - >>> # field.plot( - ... # "Displacement", component=None, plotter=mpc.plot(plotter=plotter) - ... # ).show() + >>> field.plot("Displacement", component=None, plotter=mpc.plot(plotter=plotter)).show() See Also -------- @@ -353,9 +351,7 @@ class MultiPointContact: ... point_size=16, ... color="green", ... ) - >>> # field.plot( - ... # "Displacement", component=None, plotter=contact.plot(plotter=plotter) - ... # ).show() + >>> field.plot("Displacement", component=None, plotter=contact.plot(plotter=plotter)).show() See Also -------- diff --git a/src/felupe/region/_boundary.py b/src/felupe/region/_boundary.py index 59971aef..9d26968b 100644 --- a/src/felupe/region/_boundary.py +++ b/src/felupe/region/_boundary.py @@ -203,7 +203,9 @@ def boundary_cells_hexahedron20(mesh): def boundary_cells_hexahedron27(mesh): "Convert the cells array of a bi-quadratic hex mesh into a boundary cells array." - cells_hexahedron20, cells_faces_hexahedron20 = boundary_cells_hexahedron20(mesh) + cells_hexahedron20, cells_faces_hexahedron20 = boundary_cells_hexahedron20( + mesh + ) # midpoints of faces (boundary) of a hexahedron i = [20, 21, 22, 23, 24, 25] @@ -497,9 +499,7 @@ class RegionBoundary(Region): >>> view = mesh_faces.view( ... cell_data={"ep": fem.math.eigvalsh(e).max(axis=0).mean(axis=-2).T} ... ) - >>> view.plot( - ... "ep", label="Strain (Max. Principal)", plotter=field.plot(style="wireframe") - ... ).show() + >>> view.plot("ep", label="Strain (Max. Principal)", plotter=field.plot(style="wireframe")).show() """ def __init__( @@ -553,7 +553,10 @@ def __init__( if mask is not None: point_selection = np.arange(len(mesh.points))[mask] self._selection = self._selection[ - np.all(np.isin(cells_faces[self._selection], point_selection), axis=1) + np.all( + np.isin(cells_faces[self._selection], point_selection), + axis=1, + ) ] # get cell-faces and cells on boundary (unique cell-faces with one count) @@ -590,7 +593,9 @@ def _init_faces(self): tangents.append(dX_1 / np.linalg.norm(dX_1, axis=0)) if self.ensure_3d: - tangents[0] = np.insert(tangents[0], len(tangents[0]), 0, axis=0) + tangents[0] = np.insert( + tangents[0], len(tangents[0]), 0, axis=0 + ) other_tangent = np.zeros_like(tangents[0]) other_tangent[2] = 1.0 tangents.append(other_tangent) @@ -602,8 +607,12 @@ def _init_faces(self): ): dA_1 = cross(self.dXdr[:, 0], self.dXdr[:, 1]) - tangents.append(-self.dXdr[:, 0] / np.linalg.norm(self.dXdr[:, 0], axis=0)) - tangents.append(self.dXdr[:, 1] / np.linalg.norm(self.dXdr[:, 1], axis=0)) + tangents.append( + -self.dXdr[:, 0] / np.linalg.norm(self.dXdr[:, 0], axis=0) + ) + tangents.append( + self.dXdr[:, 1] / np.linalg.norm(self.dXdr[:, 1], axis=0) + ) dA = -dA_1 * self.quadrature.weights.reshape(-1, 1) From 2262e6126f0ae667a531fd9d98a4fcbf0cab415d Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Wed, 14 May 2025 14:51:20 +0200 Subject: [PATCH 20/40] Update project.rst --- docs/howto/project.rst | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/docs/howto/project.rst b/docs/howto/project.rst index 2084d4f6..fd367422 100644 --- a/docs/howto/project.rst +++ b/docs/howto/project.rst @@ -26,25 +26,25 @@ Different methods may be used to *move* the cell-data to the mesh-points. :context: :force_static: - import pyvista as pv + # import pyvista as pv - plotter = pv.Plotter(shape=(2, 2)) - kwargs = dict(name="Cauchy Stress", component=1, plotter=plotter) + # plotter = pv.Plotter(shape=(2, 2)) + # kwargs = dict(name="Cauchy Stress", component=1, plotter=plotter) - plotter.subplot(0, 0) - kwargs_sbar = dict(interactive=False, title="Cauchy Stress YY (None)") - solid.plot(project=None, **kwargs, scalar_bar_args=kwargs_sbar) + # plotter.subplot(0, 0) + # kwargs_sbar = dict(interactive=False, title="Cauchy Stress YY (None)") + # solid.plot(project=None, **kwargs, scalar_bar_args=kwargs_sbar) - plotter.subplot(0, 1) - kwargs_sbar = dict(interactive=False, title="Cauchy Stress YY (topoints)") - solid.plot(project=fem.topoints, **kwargs, scalar_bar_args=kwargs_sbar) + # plotter.subplot(0, 1) + # kwargs_sbar = dict(interactive=False, title="Cauchy Stress YY (topoints)") + # solid.plot(project=fem.topoints, **kwargs, scalar_bar_args=kwargs_sbar) - plotter.subplot(1, 0) - kwargs_sbar = dict(interactive=False, title="Cauchy Stress YY (project)") - solid.plot(project=fem.project, **kwargs, scalar_bar_args=kwargs_sbar) + # plotter.subplot(1, 0) + # kwargs_sbar = dict(interactive=False, title="Cauchy Stress YY (project)") + # solid.plot(project=fem.project, **kwargs, scalar_bar_args=kwargs_sbar) - plotter.subplot(1, 1) - kwargs_sbar = dict(interactive=False, title="Cauchy Stress YY (extrapolate)") - solid.plot(project=fem.tools.extrapolate, **kwargs, scalar_bar_args=kwargs_sbar) + # plotter.subplot(1, 1) + # kwargs_sbar = dict(interactive=False, title="Cauchy Stress YY (extrapolate)") + # solid.plot(project=fem.tools.extrapolate, **kwargs, scalar_bar_args=kwargs_sbar) - plotter.show() \ No newline at end of file + # plotter.show() \ No newline at end of file From 7c9fd6a4b68938a651a0e40d6da39973a3fa7145 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Wed, 14 May 2025 22:40:13 +0200 Subject: [PATCH 21/40] Enhance some docstrings (only minor changes) (#961) --- src/felupe/mechanics/_curve.py | 7 +++++++ src/felupe/mechanics/_free_vibration.py | 1 + src/felupe/mechanics/_job.py | 1 + src/felupe/mechanics/_step.py | 8 ++++++++ src/felupe/mesh/_geometry.py | 2 +- 5 files changed, 18 insertions(+), 1 deletion(-) diff --git a/src/felupe/mechanics/_curve.py b/src/felupe/mechanics/_curve.py index 6f672d35..4242f60e 100644 --- a/src/felupe/mechanics/_curve.py +++ b/src/felupe/mechanics/_curve.py @@ -48,6 +48,8 @@ class CharacteristicCurve(Job): Examples -------- .. pyvista-plot:: + :context: + :force_static: >>> import felupe as fem >>> @@ -72,6 +74,11 @@ class CharacteristicCurve(Job): ... ylabel=r"Normal Force in $F_1$ in N $\rightarrow$", ... marker="o", ... ) + + .. pyvista-plot:: + :context: + :force_static: + >>> solid.plot("Principal Values of Cauchy Stress").show() See Also diff --git a/src/felupe/mechanics/_free_vibration.py b/src/felupe/mechanics/_free_vibration.py index 4c27ae79..9f06745a 100644 --- a/src/felupe/mechanics/_free_vibration.py +++ b/src/felupe/mechanics/_free_vibration.py @@ -43,6 +43,7 @@ class FreeVibration: Examples -------- .. pyvista-plot:: + :force_static: >>> import felupe as fem >>> import numpy as np diff --git a/src/felupe/mechanics/_job.py b/src/felupe/mechanics/_job.py index 84783c1f..98c76815 100644 --- a/src/felupe/mechanics/_job.py +++ b/src/felupe/mechanics/_job.py @@ -86,6 +86,7 @@ class Job: Examples -------- .. pyvista-plot:: + :force_static: >>> import felupe as fem >>> diff --git a/src/felupe/mechanics/_step.py b/src/felupe/mechanics/_step.py index 95e1c40e..1a802963 100644 --- a/src/felupe/mechanics/_step.py +++ b/src/felupe/mechanics/_step.py @@ -39,6 +39,7 @@ class Step: Examples -------- .. pyvista-plot:: + :force_static: >>> import felupe as fem >>> @@ -58,6 +59,13 @@ class Step: >>> >>> job = fem.Job(steps=[step]).evaluate() >>> ax = solid.plot("Principal Values of Cauchy Stress").show() + + See Also + -------- + Job : A job with a list of steps and a method to evaluate them. + CharacteristicCurve : A job with a list of steps and a method to evaluate them. + Force-displacement curve data is tracked during evaluation for a given + :class:`~felupe.Boundary`. """ def __init__(self, items, ramp=None, boundaries=None): diff --git a/src/felupe/mesh/_geometry.py b/src/felupe/mesh/_geometry.py index f78bf24b..d6fbe6c4 100644 --- a/src/felupe/mesh/_geometry.py +++ b/src/felupe/mesh/_geometry.py @@ -321,7 +321,7 @@ def __init__(self, a=(0.0, 0.0), b=(1.0, 1.0), order=2): class CubeArbitraryOrderHexahedron(Cube): - """A rectangular 2d-mesh with an arbitrary-order Lagrange hexahedron between ``a`` + """A rectangular 3d-mesh with an arbitrary-order Lagrange hexahedron between ``a`` and ``b``. Examples From 79d51f7f105722a2587b27845225469e549e6f29 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Thu, 15 May 2025 07:17:14 +0200 Subject: [PATCH 22/40] Docs: Fix the docstring-plot in `CharacteristicCurve` --- src/felupe/mechanics/_curve.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/src/felupe/mechanics/_curve.py b/src/felupe/mechanics/_curve.py index 4242f60e..f8de7e55 100644 --- a/src/felupe/mechanics/_curve.py +++ b/src/felupe/mechanics/_curve.py @@ -75,6 +75,17 @@ class CharacteristicCurve(Job): ... marker="o", ... ) + .. pyvista-plot:: + :include-source: False + :context: + :force_static: + + >>> import pyvista as pv + >>> + >>> fig = ax.get_figure() + >>> chart = pv.ChartMPL(fig) + >>> chart.show() + .. pyvista-plot:: :context: :force_static: From aa837a64d75aad3b6868fabbfb3deb1b8bab8c3c Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Thu, 15 May 2025 21:28:15 +0200 Subject: [PATCH 23/40] Docs: Use matplotlib-directive for plots (#962) as this is now supported with PyVista's directive --- docs/conf.py | 7 ++- src/felupe/constitution/_base.py | 30 ++----------- src/felupe/constitution/_material.py | 8 ++-- src/felupe/constitution/_mixed.py | 16 ++++--- src/felupe/constitution/_view.py | 44 +++---------------- .../_neo_hooke_compressible.py | 14 +----- .../_neo_hooke_nearly_incompressible.py | 14 +----- .../hyperelasticity/_ogden_roxburgh.py | 14 +----- src/felupe/constitution/jax/_hyperelastic.py | 14 +----- src/felupe/constitution/jax/_material.py | 14 +----- .../linear_elasticity/_linear_elastic.py | 32 ++------------ .../_linear_elastic_large_strain.py | 14 +----- .../_linear_elastic_orthotropic.py | 14 +----- src/felupe/constitution/poisson/_laplace.py | 12 ----- .../small_strain/models/_linear_elastic.py | 14 +----- .../_linear_elastic_plastic_isotropic.py | 14 +----- .../constitution/tensortrax/_hyperelastic.py | 14 +----- .../constitution/tensortrax/_material.py | 14 +----- .../models/hyperelastic/_alexander.py | 14 +----- .../hyperelastic/_anssari_benam_bucchi.py | 14 +----- .../models/hyperelastic/_arruda_boyce.py | 14 +----- .../models/hyperelastic/_blatz_ko.py | 19 ++------ .../models/hyperelastic/_extended_tube.py | 19 ++------ .../_finite_strain_viscoelastic.py | 14 +----- .../models/hyperelastic/_lopez_pamies.py | 14 +----- .../hyperelastic/_miehe_goektepe_lulei.py | 19 ++------ .../models/hyperelastic/_mooney_rivlin.py | 19 ++------ .../_morph_representative_directions.py | 14 +----- .../models/hyperelastic/_neo_hooke.py | 19 ++------ .../tensortrax/models/hyperelastic/_ogden.py | 14 +----- .../models/hyperelastic/_ogden_roxburgh.py | 14 +----- .../hyperelastic/_saint_venant_kirchhoff.py | 15 +------ .../_saint_venant_kirchhoff_orthotropic.py | 15 +------ .../models/hyperelastic/_storakers.py | 19 ++------ .../hyperelastic/_third_order_deformation.py | 19 ++------ .../models/hyperelastic/_van_der_waals.py | 19 ++------ .../tensortrax/models/hyperelastic/_yeoh.py | 19 ++------ src/felupe/mechanics/_curve.py | 9 +--- 38 files changed, 90 insertions(+), 535 deletions(-) diff --git a/docs/conf.py b/docs/conf.py index f6628a85..e4fef21d 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -59,6 +59,7 @@ "sphinx_copybutton", "sphinx_design", "sphinx_gallery.gen_gallery", + "matplotlib.sphinxext.plot_directive", "pyvista.ext.plot_directive", "pyvista.ext.viewer_directive", ] @@ -99,11 +100,15 @@ # Execution mode for notebooks # nb_execution_mode = "force" -# plot directives +# matplotlib plot directive plot_include_source = True plot_html_show_source_link = False +plot_html_show_formats = False plot_formats = ["png"] +# pyvista plot directive +pyvista_plot_include_source = True + # -- Options for HTML output ------------------------------------------------- # The theme to use for HTML and HTML Help pages. See the documentation for diff --git a/src/felupe/constitution/_base.py b/src/felupe/constitution/_base.py index f7bc9a77..e3ff9e9a 100644 --- a/src/felupe/constitution/_base.py +++ b/src/felupe/constitution/_base.py @@ -303,8 +303,7 @@ def optimize( The :func:`Anssari-Benam Bucchi ` material model formulation is best-fitted on Treloar's uniaxial and biaxial tension data [1]_. - .. pyvista-plot:: - :context: + .. plot:: >>> import numpy as np >>> import felupe as fem @@ -334,17 +333,6 @@ def optimize( >>> ax = umat_new.plot(incompressible=True, ux=ux, bx=None, ps=None) >>> ax.plot(λ, P, "C0x") - .. pyvista-plot:: - :include-source: False - :context: - :force_static: - - >>> import pyvista as pv - >>> - >>> fig = ax.get_figure() - >>> chart = pv.ChartMPL(fig) - >>> chart.show() - See Also -------- scipy.optimize.least_squares : Solve a nonlinear least-squares problem with @@ -447,8 +435,7 @@ def constitutive_material(Material, name=None): This example shows how to create a derived user material class to enable the methods from :class:`~felupe.ConstitutiveMaterial` on any (external) material. - .. pyvista-plot:: - :context: + .. plot:: >>> import felupe as fem >>> import numpy as np @@ -472,17 +459,6 @@ def constitutive_material(Material, name=None): >>> umat = MyMaterial(a=0.5) >>> ax = umat.plot(incompressible=True) - .. pyvista-plot:: - :include-source: False - :context: - :force_static: - - >>> import pyvista as pv - >>> - >>> fig = ax.get_figure() - >>> chart = pv.ChartMPL(fig) - >>> chart.show() - See Also -------- felupe.ConstitutiveMaterial : Base class for constitutive material formulations. @@ -518,7 +494,7 @@ class CompositeMaterial(ConstitutiveMaterial): Examples -------- - .. pyvista-plot:: + .. plot:: >>> import felupe as fem >>> diff --git a/src/felupe/constitution/_material.py b/src/felupe/constitution/_material.py index f84ad311..ba9b9fd9 100644 --- a/src/felupe/constitution/_material.py +++ b/src/felupe/constitution/_material.py @@ -223,8 +223,8 @@ def hessian(x, **kwargs): \boldsymbol{F}^{-T} \overset{il}{\otimes} \boldsymbol{F}^{-T} + \lambda \boldsymbol{F}^{-T} {\otimes} \boldsymbol{F}^{-T} - .. pyvista-plot:: - :context: + .. plot:: + :context: close-figs >>> import numpy as np >>> import felupe as fem @@ -261,8 +261,8 @@ def hessian(x, **kwargs): The material formulation is tested in a minimal example of non-homogeneous uniaxial tension. - .. pyvista-plot:: - :context: + .. plot:: + :context: close-figs >>> mesh = fem.Cube(n=3) >>> region = fem.RegionHexahedron(mesh) diff --git a/src/felupe/constitution/_mixed.py b/src/felupe/constitution/_mixed.py index 06894059..449d36d6 100644 --- a/src/felupe/constitution/_mixed.py +++ b/src/felupe/constitution/_mixed.py @@ -79,13 +79,15 @@ class NearlyIncompressible(ConstitutiveMaterial): Examples -------- - >>> import felupe as fem - >>> - >>> field = fem.FieldsMixed(fem.RegionHexahedron(fem.Cube(n=6)), n=3) - >>> boundaries, loadcase = fem.dof.uniaxial(field, clamped=True) - >>> umat = fem.NearlyIncompressible(fem.NeoHooke(mu=1), bulk=5000) - >>> solid = fem.SolidBody(umat, field) - >>> job = fem.Job(steps=[fem.Step(items=[solid], boundaries=boundaries)]).evaluate() + + .. plot:: + + >>> import felupe as fem + >>> + >>> field = fem.FieldsMixed(fem.RegionHexahedron(fem.Cube(n=6)), n=3) + >>> boundaries, loadcase = fem.dof.uniaxial(field, clamped=True) + >>> umat = fem.NearlyIncompressible(fem.NeoHooke(mu=1), bulk=5000) + >>> solid = fem.SolidBody(umat, field) See Also -------- diff --git a/src/felupe/constitution/_view.py b/src/felupe/constitution/_view.py index 2ac94cba..8e88fba3 100644 --- a/src/felupe/constitution/_view.py +++ b/src/felupe/constitution/_view.py @@ -152,8 +152,7 @@ class ViewMaterial(PlotMaterial): Examples -------- - .. pyvista-plot:: - :context: + .. plot:: >>> import felupe as fem >>> @@ -166,17 +165,6 @@ class ViewMaterial(PlotMaterial): ... ) >>> ax = view.plot(show_title=True, show_kwargs=True) - .. pyvista-plot:: - :include-source: False - :context: - :force_static: - - >>> import pyvista as pv - >>> - >>> fig = ax.get_figure() - >>> chart = pv.ChartMPL(fig) - >>> chart.show() - """ def __init__( @@ -440,8 +428,8 @@ class ViewMaterialIncompressible(PlotMaterial): Examples -------- - .. pyvista-plot:: - :context: + .. plot:: + :context: close-figs >>> import felupe as fem >>> @@ -449,19 +437,8 @@ class ViewMaterialIncompressible(PlotMaterial): >>> preview = fem.ViewMaterialIncompressible(umat) >>> ax = preview.plot(show_title=True, show_kwargs=True) - .. pyvista-plot:: - :include-source: False - :context: - :force_static: - - >>> import pyvista as pv - >>> - >>> fig = ax.get_figure() - >>> chart = pv.ChartMPL(fig) - >>> chart.show() - - .. pyvista-plot:: - :context: + .. plot:: + :context: close-figs >>> umat = fem.OgdenRoxburgh(fem.NeoHooke(mu=1), r=3, m=1, beta=0) >>> view = fem.ViewMaterialIncompressible( @@ -472,17 +449,6 @@ class ViewMaterialIncompressible(PlotMaterial): ... ) >>> ax = view.plot(show_title=True, show_kwargs=True) - .. pyvista-plot:: - :include-source: False - :context: - :force_static: - - >>> import pyvista as pv - >>> - >>> fig = ax.get_figure() - >>> chart = pv.ChartMPL(fig) - >>> chart.show() - """ def __init__( diff --git a/src/felupe/constitution/hyperelasticity/_neo_hooke_compressible.py b/src/felupe/constitution/hyperelasticity/_neo_hooke_compressible.py index 5f84971d..f6dc8e63 100644 --- a/src/felupe/constitution/hyperelasticity/_neo_hooke_compressible.py +++ b/src/felupe/constitution/hyperelasticity/_neo_hooke_compressible.py @@ -75,25 +75,13 @@ class NeoHookeCompressible(ConstitutiveMaterial): Examples -------- - .. pyvista-plot:: - :context: + .. plot:: >>> import felupe as fem >>> >>> umat = fem.NeoHookeCompressible(mu=1.0, lmbda=2.0) >>> ax = umat.plot() - .. pyvista-plot:: - :include-source: False - :context: - :force_static: - - >>> import pyvista as pv - >>> - >>> fig = ax.get_figure() - >>> chart = pv.ChartMPL(fig) - >>> chart.show() - """ def __init__(self, mu=None, lmbda=None, parallel=False): diff --git a/src/felupe/constitution/hyperelasticity/_neo_hooke_nearly_incompressible.py b/src/felupe/constitution/hyperelasticity/_neo_hooke_nearly_incompressible.py index 3e0bdca5..1a25c7b8 100644 --- a/src/felupe/constitution/hyperelasticity/_neo_hooke_nearly_incompressible.py +++ b/src/felupe/constitution/hyperelasticity/_neo_hooke_nearly_incompressible.py @@ -164,25 +164,13 @@ class NeoHooke(ConstitutiveMaterial): Examples -------- - .. pyvista-plot:: - :context: + .. plot:: >>> import felupe as fem >>> >>> umat = fem.NeoHooke(mu=1.0, bulk=2.0) >>> ax = umat.plot() - .. pyvista-plot:: - :include-source: False - :context: - :force_static: - - >>> import pyvista as pv - >>> - >>> fig = ax.get_figure() - >>> chart = pv.ChartMPL(fig) - >>> chart.show() - """ def __init__(self, mu=None, bulk=None, parallel=False): diff --git a/src/felupe/constitution/hyperelasticity/_ogden_roxburgh.py b/src/felupe/constitution/hyperelasticity/_ogden_roxburgh.py index 9fcefd02..1f1e7ee5 100644 --- a/src/felupe/constitution/hyperelasticity/_ogden_roxburgh.py +++ b/src/felupe/constitution/hyperelasticity/_ogden_roxburgh.py @@ -57,8 +57,7 @@ class OgdenRoxburgh(ConstitutiveMaterial): Examples -------- - .. pyvista-plot:: - :context: + .. plot:: >>> import felupe as fem >>> @@ -72,17 +71,6 @@ class OgdenRoxburgh(ConstitutiveMaterial): ... incompressible=True, ... ) - .. pyvista-plot:: - :include-source: False - :context: - :force_static: - - >>> import pyvista as pv - >>> - >>> fig = ax.get_figure() - >>> chart = pv.ChartMPL(fig) - >>> chart.show() - """ def __init__(self, material, r, m, beta): diff --git a/src/felupe/constitution/jax/_hyperelastic.py b/src/felupe/constitution/jax/_hyperelastic.py index 8460d8cc..09d9fa94 100644 --- a/src/felupe/constitution/jax/_hyperelastic.py +++ b/src/felupe/constitution/jax/_hyperelastic.py @@ -118,8 +118,7 @@ def viscoelastic(C, Cin, mu, eta, dtime): -------- View force-stretch curves on elementary incompressible deformations. - .. pyvista-plot:: - :context: + .. plot:: >>> import felupe as fem >>> import felupe.constitution.jax as mat @@ -132,17 +131,6 @@ def viscoelastic(C, Cin, mu, eta, dtime): >>> umat = mat.Hyperelastic(neo_hooke, mu=1) >>> ax = umat.plot(incompressible=True) - .. pyvista-plot:: - :include-source: False - :context: - :force_static: - - >>> import pyvista as pv - >>> - >>> fig = ax.get_figure() - >>> chart = pv.ChartMPL(fig) - >>> chart.show() - """ def __init__(self, fun, nstatevars=0, jit=True, parallel=False, **kwargs): diff --git a/src/felupe/constitution/jax/_material.py b/src/felupe/constitution/jax/_material.py index 722979db..2acab805 100644 --- a/src/felupe/constitution/jax/_material.py +++ b/src/felupe/constitution/jax/_material.py @@ -125,8 +125,7 @@ def viscoelastic(F, Cin, mu, eta, dtime): -------- View force-stretch curves on elementary incompressible deformations. - .. pyvista-plot:: - :context: + .. plot:: >>> import felupe as fem >>> import felupe.constitution.jax as mat @@ -144,17 +143,6 @@ def viscoelastic(F, Cin, mu, eta, dtime): >>> umat = mat.Material(neo_hooke, mu=1) >>> ax = umat.plot(incompressible=True) - .. pyvista-plot:: - :include-source: False - :context: - :force_static: - - >>> import pyvista as pv - >>> - >>> fig = ax.get_figure() - >>> chart = pv.ChartMPL(fig) - >>> chart.show() - """ def __init__( diff --git a/src/felupe/constitution/linear_elasticity/_linear_elastic.py b/src/felupe/constitution/linear_elasticity/_linear_elastic.py index 3821858b..43a97216 100644 --- a/src/felupe/constitution/linear_elasticity/_linear_elastic.py +++ b/src/felupe/constitution/linear_elasticity/_linear_elastic.py @@ -84,25 +84,13 @@ class LinearElastic(ConstitutiveMaterial): Examples -------- - .. pyvista-plot:: - :context: - + .. plot:: + >>> import felupe as fem >>> >>> umat = fem.LinearElastic(E=1, nu=0.3) >>> ax = umat.plot() - .. pyvista-plot:: - :include-source: False - :context: - :force_static: - - >>> import pyvista as pv - >>> - >>> fig = ax.get_figure() - >>> chart = pv.ChartMPL(fig) - >>> chart.show() - See Also -------- felupe.LinearElasticLargeStrain : Linear-elastic material formulation suitable for @@ -242,24 +230,12 @@ class LinearElasticTensorNotation(ConstitutiveMaterial): Examples -------- - .. pyvista-plot:: - :context: - + .. plot:: + >>> import felupe as fem >>> >>> umat = fem.constitution.LinearElasticTensorNotation(E=1, nu=0.3) >>> ax = umat.plot() - - .. pyvista-plot:: - :include-source: False - :context: - :force_static: - - >>> import pyvista as pv - >>> - >>> fig = ax.get_figure() - >>> chart = pv.ChartMPL(fig) - >>> chart.show() """ diff --git a/src/felupe/constitution/linear_elasticity/_linear_elastic_large_strain.py b/src/felupe/constitution/linear_elasticity/_linear_elastic_large_strain.py index a7b80ef8..21a5699a 100644 --- a/src/felupe/constitution/linear_elasticity/_linear_elastic_large_strain.py +++ b/src/felupe/constitution/linear_elasticity/_linear_elastic_large_strain.py @@ -41,25 +41,13 @@ class LinearElasticLargeStrain(ConstitutiveMaterial): Examples -------- - .. pyvista-plot:: - :context: + .. plot:: >>> import felupe as fem >>> >>> umat = fem.LinearElasticLargeStrain(E=1.0, nu=0.3) >>> ax = umat.plot() - .. pyvista-plot:: - :include-source: False - :context: - :force_static: - - >>> import pyvista as pv - >>> - >>> fig = ax.get_figure() - >>> chart = pv.ChartMPL(fig) - >>> chart.show() - """ def __init__(self, E, nu, parallel=False): diff --git a/src/felupe/constitution/linear_elasticity/_linear_elastic_orthotropic.py b/src/felupe/constitution/linear_elasticity/_linear_elastic_orthotropic.py index 62967864..d597b589 100644 --- a/src/felupe/constitution/linear_elasticity/_linear_elastic_orthotropic.py +++ b/src/felupe/constitution/linear_elasticity/_linear_elastic_orthotropic.py @@ -36,8 +36,7 @@ class LinearElasticOrthotropic(ConstitutiveMaterial): Examples -------- - .. pyvista-plot:: - :context: + .. plot:: >>> import felupe as fem >>> @@ -46,17 +45,6 @@ class LinearElasticOrthotropic(ConstitutiveMaterial): ... ) >>> ax = umat.plot() - .. pyvista-plot:: - :include-source: False - :context: - :force_static: - - >>> import pyvista as pv - >>> - >>> fig = ax.get_figure() - >>> chart = pv.ChartMPL(fig) - >>> chart.show() - """ def __init__(self, E, nu, G): diff --git a/src/felupe/constitution/poisson/_laplace.py b/src/felupe/constitution/poisson/_laplace.py index 44b44c93..71e4e45f 100644 --- a/src/felupe/constitution/poisson/_laplace.py +++ b/src/felupe/constitution/poisson/_laplace.py @@ -49,24 +49,12 @@ class Laplace(ConstitutiveMaterial): Examples -------- .. pyvista-plot:: - :context: >>> import felupe as fem >>> >>> umat = fem.Laplace() >>> ax = umat.plot() - .. pyvista-plot:: - :include-source: False - :context: - :force_static: - - >>> import pyvista as pv - >>> - >>> fig = ax.get_figure() - >>> chart = pv.ChartMPL(fig) - >>> chart.show() - """ def __init__(self, multiplier=1.0): diff --git a/src/felupe/constitution/small_strain/models/_linear_elastic.py b/src/felupe/constitution/small_strain/models/_linear_elastic.py index 492b62cf..dc7353f9 100644 --- a/src/felupe/constitution/small_strain/models/_linear_elastic.py +++ b/src/felupe/constitution/small_strain/models/_linear_elastic.py @@ -68,25 +68,13 @@ def linear_elastic(dε, εn, σn, ζn, λ, μ, **kwargs): Examples -------- - .. pyvista-plot:: - :context: + .. plot:: >>> import felupe as fem >>> >>> umat = fem.MaterialStrain(material=fem.linear_elastic, λ=2.0, μ=1.0) >>> ax = umat.plot() - .. pyvista-plot:: - :include-source: False - :context: - :force_static: - - >>> import pyvista as pv - >>> - >>> fig = ax.get_figure() - >>> chart = pv.ChartMPL(fig) - >>> chart.show() - See Also -------- MaterialStrain : A strain-based user-defined material definition with a given diff --git a/src/felupe/constitution/small_strain/models/_linear_elastic_plastic_isotropic.py b/src/felupe/constitution/small_strain/models/_linear_elastic_plastic_isotropic.py index 8faff824..7f2a45d2 100644 --- a/src/felupe/constitution/small_strain/models/_linear_elastic_plastic_isotropic.py +++ b/src/felupe/constitution/small_strain/models/_linear_elastic_plastic_isotropic.py @@ -119,8 +119,7 @@ def linear_elastic_plastic_isotropic_hardening(dε, εn, σn, ζn, λ, μ, σy, Examples -------- - .. pyvista-plot:: - :context: + .. plot:: >>> import felupe as fem >>> @@ -136,17 +135,6 @@ def linear_elastic_plastic_isotropic_hardening(dε, εn, σn, ζn, λ, μ, σy, >>> ux = fem.math.linsteps([1, 1.05, 0.95, 1.05], num=[10, 20, 20]) >>> ax = umat.plot(ux=ux, bx=None, ps=None) - .. pyvista-plot:: - :include-source: False - :context: - :force_static: - - >>> import pyvista as pv - >>> - >>> fig = ax.get_figure() - >>> chart = pv.ChartMPL(fig) - >>> chart.show() - See Also -------- MaterialStrain : A strain-based user-defined material definition with a given diff --git a/src/felupe/constitution/tensortrax/_hyperelastic.py b/src/felupe/constitution/tensortrax/_hyperelastic.py index 0fbd179b..29d899b7 100644 --- a/src/felupe/constitution/tensortrax/_hyperelastic.py +++ b/src/felupe/constitution/tensortrax/_hyperelastic.py @@ -112,8 +112,7 @@ def viscoelastic(C, Cin, mu, eta, dtime): -------- View force-stretch curves on elementary incompressible deformations. - .. pyvista-plot:: - :context: + .. plot:: >>> import felupe as fem >>> import felupe.constitution.tensortrax as mat @@ -126,17 +125,6 @@ def viscoelastic(C, Cin, mu, eta, dtime): >>> umat = mat.Hyperelastic(neo_hooke, mu=1) >>> ax = umat.plot(incompressible=True) - .. pyvista-plot:: - :include-source: False - :context: - :force_static: - - >>> import pyvista as pv - >>> - >>> fig = ax.get_figure() - >>> chart = pv.ChartMPL(fig) - >>> chart.show() - See Also -------- felupe.constitution.tensortrax.models.hyperelastic.saint_venant_kirchhoff : Strain diff --git a/src/felupe/constitution/tensortrax/_material.py b/src/felupe/constitution/tensortrax/_material.py index 4b6e2368..f210adc7 100644 --- a/src/felupe/constitution/tensortrax/_material.py +++ b/src/felupe/constitution/tensortrax/_material.py @@ -112,8 +112,7 @@ def viscoelastic(F, Cin, mu, eta, dtime): Examples -------- - .. pyvista-plot:: - :context: + .. plot:: >>> import felupe as fem >>> import felupe.constitution.tensortrax as mat @@ -127,17 +126,6 @@ def viscoelastic(F, Cin, mu, eta, dtime): >>> umat = mat.Material(neo_hooke, mu=1) >>> ax = umat.plot(incompressible=True) - .. pyvista-plot:: - :include-source: False - :context: - :force_static: - - >>> import pyvista as pv - >>> - >>> fig = ax.get_figure() - >>> chart = pv.ChartMPL(fig) - >>> chart.show() - """ def __init__(self, fun, nstatevars=0, parallel=False, **kwargs): diff --git a/src/felupe/constitution/tensortrax/models/hyperelastic/_alexander.py b/src/felupe/constitution/tensortrax/models/hyperelastic/_alexander.py index 9e5c450a..2a524536 100644 --- a/src/felupe/constitution/tensortrax/models/hyperelastic/_alexander.py +++ b/src/felupe/constitution/tensortrax/models/hyperelastic/_alexander.py @@ -85,8 +85,7 @@ def alexander(C, C1, C2, C3, gamma, k): Examples -------- - .. pyvista-plot:: - :context: + .. plot:: >>> import felupe as fem >>> @@ -99,17 +98,6 @@ def alexander(C, C1, C2, C3, gamma, k): >>> >>> ax = umat.plot(ux=ux, ps=ps, bx=bx, incompressible=True) - .. pyvista-plot:: - :include-source: False - :context: - :force_static: - - >>> import pyvista as pv - >>> - >>> fig = ax.get_figure() - >>> chart = pv.ChartMPL(fig) - >>> chart.show() - References ---------- .. [1] H. Alexander, "A constitutive relation for rubber-like materials", diff --git a/src/felupe/constitution/tensortrax/models/hyperelastic/_anssari_benam_bucchi.py b/src/felupe/constitution/tensortrax/models/hyperelastic/_anssari_benam_bucchi.py index a5a1e8a8..1e5ba4fa 100644 --- a/src/felupe/constitution/tensortrax/models/hyperelastic/_anssari_benam_bucchi.py +++ b/src/felupe/constitution/tensortrax/models/hyperelastic/_anssari_benam_bucchi.py @@ -61,8 +61,7 @@ def anssari_benam_bucchi(C, mu, N): Examples -------- - .. pyvista-plot:: - :context: + .. plot:: >>> import felupe as fem >>> @@ -74,17 +73,6 @@ def anssari_benam_bucchi(C, mu, N): >>> >>> ax = umat.plot(ux=ux, ps=ps, bx=bx, incompressible=True) - .. pyvista-plot:: - :include-source: False - :context: - :force_static: - - >>> import pyvista as pv - >>> - >>> fig = ax.get_figure() - >>> chart = pv.ChartMPL(fig) - >>> chart.show() - References ---------- .. [1] A. Anssari-Benam and A. Bucchi, "A generalised neo-Hookean strain energy diff --git a/src/felupe/constitution/tensortrax/models/hyperelastic/_arruda_boyce.py b/src/felupe/constitution/tensortrax/models/hyperelastic/_arruda_boyce.py index b1dee4f2..e49ae8bf 100644 --- a/src/felupe/constitution/tensortrax/models/hyperelastic/_arruda_boyce.py +++ b/src/felupe/constitution/tensortrax/models/hyperelastic/_arruda_boyce.py @@ -82,25 +82,13 @@ def arruda_boyce(C, C1, limit): Examples -------- - .. pyvista-plot:: - :context: + .. plot:: >>> import felupe as fem >>> >>> umat = fem.Hyperelastic(fem.arruda_boyce, C1=1.0, limit=3.2) >>> ax = umat.plot(incompressible=True) - .. pyvista-plot:: - :include-source: False - :context: - :force_static: - - >>> import pyvista as pv - >>> - >>> fig = ax.get_figure() - >>> chart = pv.ChartMPL(fig) - >>> chart.show() - """ I1 = det(C) ** (-1 / 3) * trace(C) diff --git a/src/felupe/constitution/tensortrax/models/hyperelastic/_blatz_ko.py b/src/felupe/constitution/tensortrax/models/hyperelastic/_blatz_ko.py index 1444b531..ad498276 100644 --- a/src/felupe/constitution/tensortrax/models/hyperelastic/_blatz_ko.py +++ b/src/felupe/constitution/tensortrax/models/hyperelastic/_blatz_ko.py @@ -53,33 +53,22 @@ def blatz_ko(C, mu): -------- First, choose the desired automatic differentiation backend - .. pyvista-plot:: - :context: + .. plot:: + :context: close-figs >>> # import felupe.constitution.jax as mat >>> import felupe.constitution.tensortrax as mat and create the hyperelastic material. - .. pyvista-plot:: - :context: + .. plot:: + :context: close-figs >>> import felupe as fem >>> >>> umat = mat.Hyperelastic(mat.models.hyperelastic.blatz_ko, mu=1.0) >>> ax = umat.plot() - .. pyvista-plot:: - :include-source: False - :context: - :force_static: - - >>> import pyvista as pv - >>> - >>> fig = ax.get_figure() - >>> chart = pv.ChartMPL(fig) - >>> chart.show() - References ---------- .. [1] P. J. Blatz and W. L. Ko, "Application of Finite Elastic Theory to the diff --git a/src/felupe/constitution/tensortrax/models/hyperelastic/_extended_tube.py b/src/felupe/constitution/tensortrax/models/hyperelastic/_extended_tube.py index 8ac61fb7..a6f125f6 100644 --- a/src/felupe/constitution/tensortrax/models/hyperelastic/_extended_tube.py +++ b/src/felupe/constitution/tensortrax/models/hyperelastic/_extended_tube.py @@ -84,16 +84,16 @@ def extended_tube(C, Gc, delta, Ge, beta): -------- First, choose the desired automatic differentiation backend - .. pyvista-plot:: - :context: + .. plot:: + :context: close-figs >>> # import felupe.constitution.jax as mat >>> import felupe.constitution.tensortrax as mat and create the hyperelastic material. - .. pyvista-plot:: - :context: + .. plot:: + :context: close-figs >>> import felupe as fem >>> @@ -106,17 +106,6 @@ def extended_tube(C, Gc, delta, Ge, beta): ... ) >>> ax = umat.plot(incompressible=True) - .. pyvista-plot:: - :include-source: False - :context: - :force_static: - - >>> import pyvista as pv - >>> - >>> fig = ax.get_figure() - >>> chart = pv.ChartMPL(fig) - >>> chart.show() - References ---------- .. [1] M. Kaliske and G. Heinrich, "An Extended Tube-Model for Rubber Elasticity: diff --git a/src/felupe/constitution/tensortrax/models/hyperelastic/_finite_strain_viscoelastic.py b/src/felupe/constitution/tensortrax/models/hyperelastic/_finite_strain_viscoelastic.py index a9580b11..72c408db 100644 --- a/src/felupe/constitution/tensortrax/models/hyperelastic/_finite_strain_viscoelastic.py +++ b/src/felupe/constitution/tensortrax/models/hyperelastic/_finite_strain_viscoelastic.py @@ -74,8 +74,7 @@ def finite_strain_viscoelastic(C, Cin, mu, eta, dtime): Examples -------- - .. pyvista-plot:: - :context: + .. plot:: >>> import felupe as fem >>> @@ -89,17 +88,6 @@ def finite_strain_viscoelastic(C, Cin, mu, eta, dtime): ... bx=None, ... ) - .. pyvista-plot:: - :include-source: False - :context: - :force_static: - - >>> import pyvista as pv - >>> - >>> fig = ax.get_figure() - >>> chart = pv.ChartMPL(fig) - >>> chart.show() - References ---------- .. [1] A. V. Shutov, R. Landgraf, and J. Ihlemann, "An explicit solution for diff --git a/src/felupe/constitution/tensortrax/models/hyperelastic/_lopez_pamies.py b/src/felupe/constitution/tensortrax/models/hyperelastic/_lopez_pamies.py index db97fcce..0f955023 100644 --- a/src/felupe/constitution/tensortrax/models/hyperelastic/_lopez_pamies.py +++ b/src/felupe/constitution/tensortrax/models/hyperelastic/_lopez_pamies.py @@ -64,8 +64,7 @@ def lopez_pamies(C, mu, alpha): Examples -------- - .. pyvista-plot:: - :context: + .. plot:: >>> import felupe as fem >>> @@ -79,17 +78,6 @@ def lopez_pamies(C, mu, alpha): >>> >>> ax = umat.plot(ux=ux, ps=ps, bx=bx, incompressible=True) - .. pyvista-plot:: - :include-source: False - :context: - :force_static: - - >>> import pyvista as pv - >>> - >>> fig = ax.get_figure() - >>> chart = pv.ChartMPL(fig) - >>> chart.show() - References ---------- .. [1] O. Lopez-Pamies, "A new I1-based hyperelastic model for rubber elastic diff --git a/src/felupe/constitution/tensortrax/models/hyperelastic/_miehe_goektepe_lulei.py b/src/felupe/constitution/tensortrax/models/hyperelastic/_miehe_goektepe_lulei.py index ce710994..c688f4c3 100644 --- a/src/felupe/constitution/tensortrax/models/hyperelastic/_miehe_goektepe_lulei.py +++ b/src/felupe/constitution/tensortrax/models/hyperelastic/_miehe_goektepe_lulei.py @@ -41,16 +41,16 @@ def miehe_goektepe_lulei(C, mu, N, U, p, q): -------- First, choose the desired automatic differentiation backend - .. pyvista-plot:: - :context: + .. plot:: + :context: close-figs >>> # import felupe.constitution.jax as mat >>> import felupe.constitution.tensortrax as mat and create the hyperelastic material. - .. pyvista-plot:: - :context: + .. plot:: + :context: close-figs >>> import felupe as fem >>> @@ -66,17 +66,6 @@ def miehe_goektepe_lulei(C, mu, N, U, p, q): >>> bx = fem.math.linsteps([1, 1.5], num=50) >>> ax = umat.plot(ux=ux, ps=ps, bx=bx, incompressible=True) - .. pyvista-plot:: - :include-source: False - :context: - :force_static: - - >>> import pyvista as pv - >>> - >>> fig = ax.get_figure() - >>> chart = pv.ChartMPL(fig) - >>> chart.show() - References ---------- .. [1] C. Miehe, S. Göktepe and F. Lulei, "A micro-macro approach to rubber-like diff --git a/src/felupe/constitution/tensortrax/models/hyperelastic/_mooney_rivlin.py b/src/felupe/constitution/tensortrax/models/hyperelastic/_mooney_rivlin.py index 21931908..fb1aefdd 100644 --- a/src/felupe/constitution/tensortrax/models/hyperelastic/_mooney_rivlin.py +++ b/src/felupe/constitution/tensortrax/models/hyperelastic/_mooney_rivlin.py @@ -68,33 +68,22 @@ def mooney_rivlin(C, C10, C01): -------- First, choose the desired automatic differentiation backend - .. pyvista-plot:: - :context: + .. plot:: + :context: close-figs >>> # import felupe.constitution.jax as mat >>> import felupe.constitution.tensortrax as mat and create the hyperelastic material. - .. pyvista-plot:: - :context: + .. plot:: + :context: close-figs >>> umat = mat.Hyperelastic( ... mat.models.hyperelastic.mooney_rivlin, C10=0.3, C01=0.8 ... ) >>> ax = umat.plot(incompressible=True) - .. pyvista-plot:: - :include-source: False - :context: - :force_static: - - >>> import pyvista as pv - >>> - >>> fig = ax.get_figure() - >>> chart = pv.ChartMPL(fig) - >>> chart.show() - """ J3 = det(C) ** (-1 / 3) I1 = J3 * trace(C) diff --git a/src/felupe/constitution/tensortrax/models/hyperelastic/_morph_representative_directions.py b/src/felupe/constitution/tensortrax/models/hyperelastic/_morph_representative_directions.py index c9ca837e..94142285 100644 --- a/src/felupe/constitution/tensortrax/models/hyperelastic/_morph_representative_directions.py +++ b/src/felupe/constitution/tensortrax/models/hyperelastic/_morph_representative_directions.py @@ -39,8 +39,7 @@ def morph_representative_directions(C, statevars, p, ε=1e-8): Examples -------- - .. pyvista-plot:: - :context: + .. plot:: >>> import felupe as fem >>> import felupe.constitution.tensortrax.models as models @@ -61,17 +60,6 @@ def morph_representative_directions(C, statevars, p, ε=1e-8): ... bx=None, ... ) - .. pyvista-plot:: - :include-source: False - :context: - :force_static: - - >>> import pyvista as pv - >>> - >>> fig = ax.get_figure() - >>> chart = pv.ChartMPL(fig) - >>> chart.show() - References ---------- .. [1] D. Besdo and J. Ihlemann, "A phenomenological constitutive model for diff --git a/src/felupe/constitution/tensortrax/models/hyperelastic/_neo_hooke.py b/src/felupe/constitution/tensortrax/models/hyperelastic/_neo_hooke.py index 2c9ccf70..5268c354 100644 --- a/src/felupe/constitution/tensortrax/models/hyperelastic/_neo_hooke.py +++ b/src/felupe/constitution/tensortrax/models/hyperelastic/_neo_hooke.py @@ -44,32 +44,21 @@ def neo_hooke(C, mu): -------- First, choose the desired automatic differentiation backend - .. pyvista-plot:: - :context: + .. plot:: + :context: close-figs >>> # import felupe.constitution.jax as mat >>> import felupe.constitution.tensortrax as mat and create the hyperelastic material. - .. pyvista-plot:: - :context: + .. plot:: + :context: close-figs >>> import felupe as fem >>> >>> umat = mat.Hyperelastic(mat.models.hyperelastic.neo_hooke, mu=1.0) >>> ax = umat.plot(incompressible=True) - .. pyvista-plot:: - :include-source: False - :context: - :force_static: - - >>> import pyvista as pv - >>> - >>> fig = ax.get_figure() - >>> chart = pv.ChartMPL(fig) - >>> chart.show() - """ return mu / 2 * (det(C) ** (-1 / 3) * trace(C) - 3) diff --git a/src/felupe/constitution/tensortrax/models/hyperelastic/_ogden.py b/src/felupe/constitution/tensortrax/models/hyperelastic/_ogden.py index e52b4d66..4e400a49 100644 --- a/src/felupe/constitution/tensortrax/models/hyperelastic/_ogden.py +++ b/src/felupe/constitution/tensortrax/models/hyperelastic/_ogden.py @@ -56,25 +56,13 @@ def ogden(C, mu, alpha): Examples -------- - .. pyvista-plot:: - :context: + .. plot:: >>> import felupe as fem >>> >>> umat = fem.Hyperelastic(fem.ogden, mu=[1, 0.2], alpha=[1.7, -1.5]) >>> ax = umat.plot(incompressible=True) - .. pyvista-plot:: - :include-source: False - :context: - :force_static: - - >>> import pyvista as pv - >>> - >>> fig = ax.get_figure() - >>> chart = pv.ChartMPL(fig) - >>> chart.show() - """ wC = det(C) ** (-1 / 3) * eigvalsh(C) diff --git a/src/felupe/constitution/tensortrax/models/hyperelastic/_ogden_roxburgh.py b/src/felupe/constitution/tensortrax/models/hyperelastic/_ogden_roxburgh.py index 07e786af..5957e716 100644 --- a/src/felupe/constitution/tensortrax/models/hyperelastic/_ogden_roxburgh.py +++ b/src/felupe/constitution/tensortrax/models/hyperelastic/_ogden_roxburgh.py @@ -77,8 +77,7 @@ def ogden_roxburgh(C, Wmax_n, material, r, m, beta, **kwargs): Examples -------- - .. pyvista-plot:: - :context: + .. plot:: >>> import felupe as fem >>> @@ -96,17 +95,6 @@ def ogden_roxburgh(C, Wmax_n, material, r, m, beta, **kwargs): ... ) >>> ax = umat.plot(ux=ux, bx=None, ps=None, incompressible=True) - .. pyvista-plot:: - :include-source: False - :context: - :force_static: - - >>> import pyvista as pv - >>> - >>> fig = ax.get_figure() - >>> chart = pv.ChartMPL(fig) - >>> chart.show() - """ W = material(C, **kwargs) diff --git a/src/felupe/constitution/tensortrax/models/hyperelastic/_saint_venant_kirchhoff.py b/src/felupe/constitution/tensortrax/models/hyperelastic/_saint_venant_kirchhoff.py index 81a0adaa..5b2165f5 100644 --- a/src/felupe/constitution/tensortrax/models/hyperelastic/_saint_venant_kirchhoff.py +++ b/src/felupe/constitution/tensortrax/models/hyperelastic/_saint_venant_kirchhoff.py @@ -64,26 +64,13 @@ def saint_venant_kirchhoff(C, mu, lmbda, k=2): .. warning:: The Saint-Venant Kirchhoff material formulation is unstable for large strains. - .. pyvista-plot:: - :context: + .. plot:: >>> import felupe as fem >>> >>> umat = fem.Hyperelastic(fem.saint_venant_kirchhoff, mu=1.0, lmbda=20.0) >>> ax = umat.plot(incompressible=False) - .. pyvista-plot:: - :include-source: False - :context: - :force_static: - - >>> import pyvista as pv - >>> - >>> fig = ax.get_figure() - >>> chart = pv.ChartMPL(fig) - >>> chart.show() - - """ eye = base.eye diff --git a/src/felupe/constitution/tensortrax/models/hyperelastic/_saint_venant_kirchhoff_orthotropic.py b/src/felupe/constitution/tensortrax/models/hyperelastic/_saint_venant_kirchhoff_orthotropic.py index 3cd1808d..add5cb1d 100644 --- a/src/felupe/constitution/tensortrax/models/hyperelastic/_saint_venant_kirchhoff_orthotropic.py +++ b/src/felupe/constitution/tensortrax/models/hyperelastic/_saint_venant_kirchhoff_orthotropic.py @@ -68,8 +68,7 @@ def saint_venant_kirchhoff_orthotropic(C, mu, lmbda, r1, r2, r3=None, k=2): The orthotropic Saint-Venant Kirchhoff material formulation is unstable for large strains. - .. pyvista-plot:: - :context: + .. plot:: >>> import felupe.constitution.tensortrax as mat >>> import felupe as fem @@ -90,18 +89,6 @@ def saint_venant_kirchhoff_orthotropic(C, mu, lmbda, r1, r2, r3=None, k=2): ... ) >>> ax = umat.plot(ux=fem.math.linsteps([1, 1.1], num=10), ps=None, bx=None) - .. pyvista-plot:: - :include-source: False - :context: - :force_static: - - >>> import pyvista as pv - >>> - >>> fig = ax.get_figure() - >>> chart = pv.ChartMPL(fig) - >>> chart.show() - - """ eye = base.eye diff --git a/src/felupe/constitution/tensortrax/models/hyperelastic/_storakers.py b/src/felupe/constitution/tensortrax/models/hyperelastic/_storakers.py index 118606ba..4ca090e5 100644 --- a/src/felupe/constitution/tensortrax/models/hyperelastic/_storakers.py +++ b/src/felupe/constitution/tensortrax/models/hyperelastic/_storakers.py @@ -68,16 +68,16 @@ def storakers(C, mu, alpha, beta): -------- First, choose the desired automatic differentiation backend - .. pyvista-plot:: - :context: + .. plot:: + :context: close-figs >>> # import felupe.constitution.jax as mat >>> import felupe.constitution.tensortrax as mat and create the hyperelastic material. - .. pyvista-plot:: - :context: + .. plot:: + :context: close-figs >>> import felupe as fem >>> @@ -93,17 +93,6 @@ def storakers(C, mu, alpha, beta): ... bx=fem.math.linsteps([1, 1], 9), ... ) - .. pyvista-plot:: - :include-source: False - :context: - :force_static: - - >>> import pyvista as pv - >>> - >>> fig = ax.get_figure() - >>> chart = pv.ChartMPL(fig) - >>> chart.show() - References ---------- .. [1] B. Storåkers, "On material representation and constitutive branching in diff --git a/src/felupe/constitution/tensortrax/models/hyperelastic/_third_order_deformation.py b/src/felupe/constitution/tensortrax/models/hyperelastic/_third_order_deformation.py index 94105265..f10d558d 100644 --- a/src/felupe/constitution/tensortrax/models/hyperelastic/_third_order_deformation.py +++ b/src/felupe/constitution/tensortrax/models/hyperelastic/_third_order_deformation.py @@ -80,16 +80,16 @@ def third_order_deformation(C, C10, C01, C11, C20, C30): -------- First, choose the desired automatic differentiation backend - .. pyvista-plot:: - :context: + .. plot:: + :context: close-figs >>> # import felupe.constitution.jax as mat >>> import felupe.constitution.tensortrax as mat and create the hyperelastic material. - .. pyvista-plot:: - :context: + .. plot:: + :context: close-figs >>> umat = mat.Hyperelastic( ... mat.models.hyperelastic.third_order_deformation, @@ -101,17 +101,6 @@ def third_order_deformation(C, C10, C01, C11, C20, C30): ... ) >>> ax = umat.plot(incompressible=True) - .. pyvista-plot:: - :include-source: False - :context: - :force_static: - - >>> import pyvista as pv - >>> - >>> fig = ax.get_figure() - >>> chart = pv.ChartMPL(fig) - >>> chart.show() - """ J3 = det(C) ** (-1 / 3) I1 = J3 * trace(C) diff --git a/src/felupe/constitution/tensortrax/models/hyperelastic/_van_der_waals.py b/src/felupe/constitution/tensortrax/models/hyperelastic/_van_der_waals.py index 1ad9875e..caa06be8 100644 --- a/src/felupe/constitution/tensortrax/models/hyperelastic/_van_der_waals.py +++ b/src/felupe/constitution/tensortrax/models/hyperelastic/_van_der_waals.py @@ -44,16 +44,16 @@ def van_der_waals(C, mu, limit, a, beta): -------- First, choose the desired automatic differentiation backend - .. pyvista-plot:: - :context: + .. plot:: + :context: close-figs >>> # import felupe.constitution.jax as mat >>> import felupe.constitution.tensortrax as mat and create the hyperelastic material. - .. pyvista-plot:: - :context: + .. plot:: + :context: close-figs >>> import felupe as fem >>> @@ -66,17 +66,6 @@ def van_der_waals(C, mu, limit, a, beta): ... ) >>> ax = umat.plot(incompressible=True) - .. pyvista-plot:: - :include-source: False - :context: - :force_static: - - >>> import pyvista as pv - >>> - >>> fig = ax.get_figure() - >>> chart = pv.ChartMPL(fig) - >>> chart.show() - References ---------- .. [1] H.-G. Kilian, "Equation of state of real networks", Polymer, vol. 22, no. 2. diff --git a/src/felupe/constitution/tensortrax/models/hyperelastic/_yeoh.py b/src/felupe/constitution/tensortrax/models/hyperelastic/_yeoh.py index bee80aca..c72b1d37 100644 --- a/src/felupe/constitution/tensortrax/models/hyperelastic/_yeoh.py +++ b/src/felupe/constitution/tensortrax/models/hyperelastic/_yeoh.py @@ -66,33 +66,22 @@ def yeoh(C, C10, C20, C30): -------- First, choose the desired automatic differentiation backend - .. pyvista-plot:: - :context: + .. plot:: + :context: close-figs >>> # import felupe.constitution.jax as mat >>> import felupe.constitution.tensortrax as mat and create the hyperelastic material. - .. pyvista-plot:: - :context: + .. plot:: + :context: close-figs >>> umat = mat.Hyperelastic( ... mat.models.hyperelastic.yeoh, C10=0.5, C20=-0.1, C30=0.02 ... ) >>> ax = umat.plot(incompressible=True) - .. pyvista-plot:: - :include-source: False - :context: - :force_static: - - >>> import pyvista as pv - >>> - >>> fig = ax.get_figure() - >>> chart = pv.ChartMPL(fig) - >>> chart.show() - """ I1 = det(C) ** (-1 / 3) * trace(C) diff --git a/src/felupe/mechanics/_curve.py b/src/felupe/mechanics/_curve.py index f8de7e55..eee98691 100644 --- a/src/felupe/mechanics/_curve.py +++ b/src/felupe/mechanics/_curve.py @@ -68,7 +68,8 @@ class CharacteristicCurve(Job): >>> move = fem.math.linsteps([0, 1], num=5) >>> step = fem.Step(items=[solid], ramp={boundaries["move"]: move}, boundaries=boundaries) >>> - >>> job = fem.CharacteristicCurve(steps=[step], boundary=boundaries["move"]).evaluate() + >>> job = fem.CharacteristicCurve(steps=[step], boundary=boundaries["move"]).evaluate()+ + >>> solid.plot("Principal Values of Cauchy Stress").show() >>> fig, ax = job.plot( ... xlabel=r"Displacement $u_1$ in mm $\rightarrow$", ... ylabel=r"Normal Force in $F_1$ in N $\rightarrow$", @@ -86,12 +87,6 @@ class CharacteristicCurve(Job): >>> chart = pv.ChartMPL(fig) >>> chart.show() - .. pyvista-plot:: - :context: - :force_static: - - >>> solid.plot("Principal Values of Cauchy Stress").show() - See Also -------- Step : A Step with multiple substeps, subsequently depending on the solution From 30394aedfedd31b327de67c893796a1e2a9f0002 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Thu, 15 May 2025 21:37:47 +0200 Subject: [PATCH 24/40] Aktualisieren von _curve.py --- src/felupe/mechanics/_curve.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/felupe/mechanics/_curve.py b/src/felupe/mechanics/_curve.py index eee98691..5391e8ea 100644 --- a/src/felupe/mechanics/_curve.py +++ b/src/felupe/mechanics/_curve.py @@ -68,7 +68,7 @@ class CharacteristicCurve(Job): >>> move = fem.math.linsteps([0, 1], num=5) >>> step = fem.Step(items=[solid], ramp={boundaries["move"]: move}, boundaries=boundaries) >>> - >>> job = fem.CharacteristicCurve(steps=[step], boundary=boundaries["move"]).evaluate()+ + >>> job = fem.CharacteristicCurve(steps=[step], boundary=boundaries["move"]).evaluate() >>> solid.plot("Principal Values of Cauchy Stress").show() >>> fig, ax = job.plot( ... xlabel=r"Displacement $u_1$ in mm $\rightarrow$", From e72247cf9f014bf7863ac363c756e60ed24a6e2b Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Thu, 15 May 2025 22:10:02 +0200 Subject: [PATCH 25/40] Remove `xvfb` in `.readthedocs.yaml` --- .readthedocs.yaml | 2 -- docs/conf.py | 6 +++--- 2 files changed, 3 insertions(+), 5 deletions(-) diff --git a/.readthedocs.yaml b/.readthedocs.yaml index 3fd675c2..8e2b6819 100644 --- a/.readthedocs.yaml +++ b/.readthedocs.yaml @@ -5,8 +5,6 @@ build: tools: python: "3.12" apt_packages: - - libgl1-mesa-dev - - xvfb - pandoc sphinx: diff --git a/docs/conf.py b/docs/conf.py index e4fef21d..a8ab4e1e 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -29,9 +29,9 @@ pyvista.BUILDING_GALLERY = True os.environ["PYVISTA_BUILDING_GALLERY"] = "true" -# start a virtual framebuffer -if os.environ.get("READTHEDOCS") or os.environ.get("CI"): - pyvista.start_xvfb() +# start a virtual framebuffer (deprecated since pyvista 0.45.0) +# if os.environ.get("READTHEDOCS") or os.environ.get("CI"): +# pyvista.start_xvfb() # FElupe: turn off logging os.environ["FELUPE_VERBOSE"] = "false" From 7272b1414d05594c185e96fc5ca8ad288aff5579 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Thu, 15 May 2025 22:31:15 +0200 Subject: [PATCH 26/40] Docs: Enforce more static scenes in API-docs --- src/felupe/mechanics/_solidbody_incompressible.py | 3 ++- src/felupe/mechanics/_solidbody_pressure.py | 1 + src/felupe/mesh/_tools.py | 4 ++-- 3 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/felupe/mechanics/_solidbody_incompressible.py b/src/felupe/mechanics/_solidbody_incompressible.py index e52aa69c..47a77ce1 100644 --- a/src/felupe/mechanics/_solidbody_incompressible.py +++ b/src/felupe/mechanics/_solidbody_incompressible.py @@ -273,7 +273,8 @@ class SolidBodyNearlyIncompressible(Solid): Examples -------- .. pyvista-plot:: - + :force_static: + >>> import felupe as fem >>> >>> mesh = fem.Cube(n=6) diff --git a/src/felupe/mechanics/_solidbody_pressure.py b/src/felupe/mechanics/_solidbody_pressure.py index f6ab684a..86b239e7 100644 --- a/src/felupe/mechanics/_solidbody_pressure.py +++ b/src/felupe/mechanics/_solidbody_pressure.py @@ -46,6 +46,7 @@ class SolidBodyPressure: Examples -------- .. pyvista-plot:: + :force_static: >>> import felupe as fem >>> diff --git a/src/felupe/mesh/_tools.py b/src/felupe/mesh/_tools.py index 0159201a..2488ee09 100644 --- a/src/felupe/mesh/_tools.py +++ b/src/felupe/mesh/_tools.py @@ -62,7 +62,7 @@ def expand(points, cells, cell_type, n=11, z=1, axis=-1, expand_dim=True): Expand a rectangle to a cube. .. pyvista-plot:: - :include-source: True + :force_static: >>> import felupe as fem >>> @@ -164,7 +164,7 @@ def fill_between(mesh, other_mesh, n=11): Examples -------- .. pyvista-plot:: - :include-source: True + :force_static: >>> import felupe as fem >>> From ae3e27c1da8da9f3d6102231964805050f88843c Mon Sep 17 00:00:00 2001 From: Tetsuo Koyama Date: Fri, 16 May 2025 23:35:51 +0900 Subject: [PATCH 27/40] Fix typo from seperated to separated (#963) --- docs/howto/composite.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/howto/composite.rst b/docs/howto/composite.rst index 1d749569..7b149b7f 100644 --- a/docs/howto/composite.rst +++ b/docs/howto/composite.rst @@ -2,7 +2,7 @@ Multiple Solid Bodies --------------------- This section demonstrates how to set up a problem with two regions, each associated to a -seperated solid body. Different element formulations are used for the solid bodies. +separated solid body. Different element formulations are used for the solid bodies. .. pyvista-plot:: :context: From 9cff56ecd770ed51e35ac4cb39a2b85d4fbc0b07 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Fri, 16 May 2025 21:24:28 +0200 Subject: [PATCH 28/40] Aktualisieren von ex01_beam.py --- examples/ex01_beam.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/examples/ex01_beam.py b/examples/ex01_beam.py index 2821f487..a022f0f6 100644 --- a/examples/ex01_beam.py +++ b/examples/ex01_beam.py @@ -15,8 +15,7 @@ .. image:: ../../examples/ex01_beam_sketch.png -First, let's create a meshed cube out of hexahedron cells with ``n=(181, 9, 9)`` points -per axis. A numeric region created on the mesh represents the cantilever beam. A three- +First, let's create a meshed cube out of hexahedron cells. A numeric region created on the mesh represents the cantilever beam. A three- dimensional vector-valued displacement field is initiated on the region. """ From 2abb3327f51704c38cb78c775f7a491ba45de4b7 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Fri, 16 May 2025 22:03:12 +0200 Subject: [PATCH 29/40] Fix indentation in benchmark source-code docs/index.rst --- docs/index.rst | 71 +++++++++++++++++++++++++------------------------- 1 file changed, 35 insertions(+), 36 deletions(-) diff --git a/docs/index.rst b/docs/index.rst index ed5686e0..6d992370 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -142,21 +142,21 @@ This is a simple benchmark to compare assembly times for linear elasticity and h def pre_linear_elastic(n, **kwargs): - mesh = fem.Cube(n=n).triangulate() - region = fem.RegionTetra(mesh) - field = fem.FieldContainer([fem.Field(region, dim=3)]) - umat = fem.LinearElastic(E=1, nu=0.3) - solid = fem.SolidBody(umat, field) - return mesh, solid + mesh = fem.Cube(n=n).triangulate() + region = fem.RegionTetra(mesh) + field = fem.FieldContainer([fem.Field(region, dim=3)]) + umat = fem.LinearElastic(E=1, nu=0.3) + solid = fem.SolidBody(umat, field) + return mesh, solid def pre_hyperelastic(n, **kwargs): - mesh = fem.Cube(n=n).triangulate() - region = fem.RegionTetra(mesh) - field = fem.FieldContainer([fem.Field(region, dim=3)]) - umat = fem.NeoHookeCompressible(mu=1.0, lmbda=2.0) - solid = fem.SolidBody(umat, field) - return mesh, solid + mesh = fem.Cube(n=n).triangulate() + region = fem.RegionTetra(mesh) + field = fem.FieldContainer([fem.Field(region, dim=3)]) + umat = fem.NeoHookeCompressible(mu=1.0, lmbda=2.0) + solid = fem.SolidBody(umat, field) + return mesh, solid print("# Assembly Runtimes") @@ -172,23 +172,23 @@ This is a simple benchmark to compare assembly times for linear elasticity and h runtimes = np.zeros((len(points_per_axis), 2)) for i, n in enumerate(points_per_axis): - mesh, solid = pre_linear_elastic(n) - matrix = solid.assemble.matrix(parallel=parallel) - time_linear_elastic = ( - timeit(lambda: solid.assemble.matrix(parallel=parallel), number=number) / number - ) + mesh, solid = pre_linear_elastic(n) + matrix = solid.assemble.matrix(parallel=parallel) + time_linear_elastic = ( + timeit(lambda: solid.assemble.matrix(parallel=parallel), number=number) / number + ) - mesh, solid = pre_hyperelastic(n) - matrix = solid.assemble.matrix(parallel=parallel) - time_hyperelastic = ( - timeit(lambda: solid.assemble.matrix(parallel=parallel), number=number) / number - ) + mesh, solid = pre_hyperelastic(n) + matrix = solid.assemble.matrix(parallel=parallel) + time_hyperelastic = ( + timeit(lambda: solid.assemble.matrix(parallel=parallel), number=number) / number + ) - runtimes[i] = time_linear_elastic, time_hyperelastic + runtimes[i] = time_linear_elastic, time_hyperelastic - print( - f"| {mesh.points.size:7d} | {runtimes[i][0]:19.2f} | {runtimes[i][1]:17.2f} |" - ) + print( + f"| {mesh.points.size:7d} | {runtimes[i][0]:19.2f} | {runtimes[i][1]:17.2f} |" + ) dofs_le = points_per_axis ** 3 * 3 / runtimes[:, 0] dofs_he = points_per_axis ** 3 * 3 / runtimes[:, 1] @@ -197,22 +197,22 @@ This is a simple benchmark to compare assembly times for linear elasticity and h print("| Analysis | DOF/s |") print("| -------------- | ----------------- |") print( - f"| Linear-Elastic | {np.mean(dofs_le):5.0f} +/-{np.std(dofs_le):5.0f} |" + f"| Linear-Elastic | {np.mean(dofs_le):5.0f} +/-{np.std(dofs_le):5.0f} |" ) print(f"| Hyperelastic | {np.mean(dofs_he):5.0f} +/-{np.std(dofs_he):5.0f} |") plt.figure() plt.loglog( - points_per_axis ** 3 * 3, - runtimes[:, 1], - "C0", - label=r"Stiffness Matrix (Hyperelastic)", + points_per_axis ** 3 * 3, + runtimes[:, 1], + "C0", + label=r"Stiffness Matrix (Hyperelastic)", ) plt.loglog( - points_per_axis ** 3 * 3, - runtimes[:, 0], - "C1--", - label=r"Stiffness Matrix (Linear-Elastic)", + points_per_axis ** 3 * 3, + runtimes[:, 0], + "C1--", + label=r"Stiffness Matrix (Linear-Elastic)", ) plt.xlabel(r"Number of degrees of freedom $\longrightarrow$") plt.ylabel(r"Runtime in s $\longrightarrow$") @@ -221,7 +221,6 @@ This is a simple benchmark to compare assembly times for linear elasticity and h plt.savefig("benchmark.png") - .. toctree:: :maxdepth: 1 :caption: Contents: From 65e6459e469b8e546cd2e1b0715e98a7426582c5 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Thu, 22 May 2025 15:53:08 +0200 Subject: [PATCH 30/40] Enhance the docs switch to matplotlib plots when possible --- .../tensortrax/_total_lagrange.py | 24 +++++++++-------- .../tensortrax/_updated_lagrange.py | 26 ++++++++++--------- .../tensortrax/models/lagrange/_morph.py | 19 +++----------- .../_morph_representative_directions.py | 19 +++----------- 4 files changed, 35 insertions(+), 53 deletions(-) diff --git a/src/felupe/constitution/tensortrax/_total_lagrange.py b/src/felupe/constitution/tensortrax/_total_lagrange.py index 98fb32fe..fc3c2aa8 100644 --- a/src/felupe/constitution/tensortrax/_total_lagrange.py +++ b/src/felupe/constitution/tensortrax/_total_lagrange.py @@ -32,17 +32,19 @@ def total_lagrange(material): Examples -------- - >>> import felupe as fem - >>> import felupe.constitution.tensortrax as mat - >>> import tensortrax.math as tm - >>> - >>> @fem.total_lagrange - >>> def neo_hooke_total_lagrange(F, mu=1): - >>> C = F.T @ F - >>> S = mu * tm.special.dev(tm.linalg.det(C)**(-1/3) * C) @ tm.linalg.inv(C) - >>> return S - >>> - >>> umat = mat.Material(neo_hooke_total_lagrange, mu=1) + .. plot:: + + >>> import felupe as fem + >>> import felupe.constitution.tensortrax as mat + >>> import tensortrax.math as tm + >>> + >>> @fem.total_lagrange + ... def neo_hooke_total_lagrange(F, mu=1): + ... C = F.T @ F + ... S = mu * tm.special.dev(tm.linalg.det(C)**(-1/3) * C) @ tm.linalg.inv(C) + ... return S + >>> + >>> umat = mat.Material(neo_hooke_total_lagrange, mu=1) See Also -------- diff --git a/src/felupe/constitution/tensortrax/_updated_lagrange.py b/src/felupe/constitution/tensortrax/_updated_lagrange.py index 6588596b..2325b053 100644 --- a/src/felupe/constitution/tensortrax/_updated_lagrange.py +++ b/src/felupe/constitution/tensortrax/_updated_lagrange.py @@ -33,18 +33,20 @@ def updated_lagrange(material): Examples -------- - >>> import felupe as fem - >>> import felupe.constitution.tensortrax as mat - >>> import tensortrax.math as tm - >>> - >>> @fem.updated_lagrange - >>> def neo_hooke_updated_lagrange(F, mu=1): - >>> J = tm.linalg.det(F) - >>> b = F @ F.T - >>> σ = mu * tm.special.dev(J**(-2/3) * b) / J - >>> return σ - >>> - >>> umat = mat.Material(neo_hooke_updated_lagrange, mu=1) + .. plot:: + + >>> import felupe as fem + >>> import felupe.constitution.tensortrax as mat + >>> import tensortrax.math as tm + >>> + >>> @fem.updated_lagrange + ... def neo_hooke_updated_lagrange(F, mu=1): + ... J = tm.linalg.det(F) + ... b = F @ F.T + ... σ = mu * tm.special.dev(J**(-2/3) * b) / J + ... return σ + >>> + >>> umat = mat.Material(neo_hooke_updated_lagrange, mu=1) See Also -------- diff --git a/src/felupe/constitution/tensortrax/models/lagrange/_morph.py b/src/felupe/constitution/tensortrax/models/lagrange/_morph.py index ce90b253..cfa84917 100644 --- a/src/felupe/constitution/tensortrax/models/lagrange/_morph.py +++ b/src/felupe/constitution/tensortrax/models/lagrange/_morph.py @@ -132,8 +132,8 @@ def morph(F, statevars, p): -------- First, choose the desired automatic differentiation backend - .. pyvista-plot:: - :context: + .. plot:: + :context: close-figs >>> import felupe as fem >>> @@ -142,8 +142,8 @@ def morph(F, statevars, p): and create the material. - .. pyvista-plot:: - :context: + .. plot:: + :context: close-figs >>> umat = mat.Material( ... mat.models.lagrange.morph, @@ -161,17 +161,6 @@ def morph(F, statevars, p): ... bx=None, ... ) - .. pyvista-plot:: - :include-source: False - :context: - :force_static: - - >>> import pyvista as pv - >>> - >>> fig = ax.get_figure() - >>> chart = pv.ChartMPL(fig) - >>> chart.show() - References ---------- .. [1] D. Besdo and J. Ihlemann, "A phenomenological constitutive model for diff --git a/src/felupe/constitution/tensortrax/models/lagrange/_morph_representative_directions.py b/src/felupe/constitution/tensortrax/models/lagrange/_morph_representative_directions.py index 88425886..270b06e4 100644 --- a/src/felupe/constitution/tensortrax/models/lagrange/_morph_representative_directions.py +++ b/src/felupe/constitution/tensortrax/models/lagrange/_morph_representative_directions.py @@ -41,8 +41,8 @@ def morph_representative_directions(F, statevars, p, ε=1e-6): -------- First, choose the desired automatic differentiation backend - .. pyvista-plot:: - :context: + .. plot:: + :context: close-figs >>> import felupe as fem >>> @@ -51,8 +51,8 @@ def morph_representative_directions(F, statevars, p, ε=1e-6): and create the material. - .. pyvista-plot:: - :context: + .. plot:: + :context: close-figs >>> umat = mat.Material( ... mat.models.lagrange.morph_representative_directions, @@ -70,17 +70,6 @@ def morph_representative_directions(F, statevars, p, ε=1e-6): ... bx=None, ... ) - .. pyvista-plot:: - :include-source: False - :context: - :force_static: - - >>> import pyvista as pv - >>> - >>> fig = ax.get_figure() - >>> chart = pv.ChartMPL(fig) - >>> chart.show() - References ---------- .. [1] D. Besdo and J. Ihlemann, "A phenomenological constitutive model for From 9fb415f5c7fbd4cca6d923f762a660ae13350730 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Fri, 23 May 2025 00:15:27 +0200 Subject: [PATCH 31/40] Add `ConstitutiveMaterial.is_stable()` and fix `CompositeMaterial` (#965) * Add `ConstitutiveMaterial.is_stable()` * Update CHANGELOG.md * Update CHANGELOG.md * Update _base.py * Change the first return value of `ViewMaterial.uniaxial()`, etc. * fix `ConstitutiveMaterial.is_stable()` * Update _base.py --- CHANGELOG.md | 3 + src/felupe/constitution/_base.py | 106 ++++++++++++++++++++++++++++++- src/felupe/constitution/_view.py | 26 +++++--- tests/test_constitution.py | 2 + 4 files changed, 126 insertions(+), 11 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index ca7f8733..59a3c923 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,16 +7,19 @@ All notable changes to this project will be documented in this file. The format - Add `Mesh.add_points(points)` to update the mesh with additional points. - Add `Mesh.clear_points_without_cells()` to clear the list of points without cells (useful for center-points of multi-point constraints). - Release FElupe on conda-forge, starting with v9.2.0. +- Add `ConstitutiveMaterial.is_stable()` which returns a boolean mask of stability for isotropic material model formulations. Note that this will require an additional volumetric part of the strain energy density function for hyperelastic material model formulations without a volumetric part. ### Changed - Change the required setuptools-version in the build-system table of `pyproject.toml` to match PEP639 (setuptools>=77.0.3). - Change the labels to well-known labels for the URLs in `pyproject.toml`. +- Change the first return value of `ViewMaterial.uniaxial()`, `ViewMaterial.planar()`, `ViewMaterial.biaxial()`, `ViewMaterialIncompressible.uniaxial()`, `ViewMaterialIncompressible.planar()`, `ViewMaterialIncompressible.biaxial()` from the stretch to a list of all three principal stretches. ### Fixed - Fix the declaration of the (spdx identifier) license and license-file in `pyproject.toml`. - Fix missing import of `TriQuadraticHexahedron` in the top-level namespace. - Fix the path to `docs/_static/logo_without_text.svg` in `docs/conf.py`. - Fix a typo in the docstring of `MeshContainer.from_unstructured_grid()`. +- Fix `CompositeMaterial` for input lists of length 1. ## [9.2.0] - 2025-03-04 diff --git a/src/felupe/constitution/_base.py b/src/felupe/constitution/_base.py index e3ff9e9a..9b4950f7 100644 --- a/src/felupe/constitution/_base.py +++ b/src/felupe/constitution/_base.py @@ -20,6 +20,7 @@ import numpy as np +from ..math import dot, eigh, transpose from ._view import ViewMaterial, ViewMaterialIncompressible @@ -407,6 +408,107 @@ def std(hessian, residuals_variance): return umat, res + def is_stable(self, x, hessian=None): + """Return a boolean mask for stability of isotropic material model formulations. + + At a given deformation gradient, a normal force is applied on each principal + stretch direction. If the resulting incremental stretches are positive, the + material model formulation is considered to be stable at the given deformation + gradient. + + Parameters + ---------- + x : list of ndarray + The list with input arguments. These contain the extracted fields of a + :class:`~felupe.FieldContainer`. + hessian : ndarray or None, optional + Second partial derivative of the strain energy density function w.r.t. the + deformation gradient. Default is None. + + Returns + ------- + ndarray + Boolean mask of stability. + + Notes + ----- + + .. warning:: + + This stability check will lead to a singular matrix for isotropic + (hyperelastic) material model formulations without a volumetric part. + + Examples + -------- + First, let's check the stability of the Neo-Hookean material model formulation. + The stability is evaluated on (valid) principal stretches of a biaxial + deformation. All deformations are stable. + + .. plot:: + + >>> import numpy as np + >>> import felupe as fem + >>> + >>> umat = fem.NeoHooke(mu=1.0, bulk=2.0) + >>> view = umat.view() + >>> λ = view.biaxial()[0] + >>> + >>> F = np.zeros((3, 3, 1, λ[0].size)) + >>> for a in range(3): + ... F[a, a] = λ[a] + >>> + >>> umat.is_stable([F]) + array([[ True, True, True, True, True, True, True, True, True, + True, True, True, True, True, True, True]]) + + Now, let's check the stability of the Mooney-Rivlin material model formulation. + The stability is evaluated on (valid) principal stretches of a biaxial + deformation. Biaxial deformations are only stable up to a longitudinal stretch + of 1.35. + + .. plot:: + + >>> import numpy as np + >>> import felupe as fem + >>> import felupe.constitution.tensortrax as mat + >>> + >>> umat = fem.Hyperelastic( + ... mat.models.hyperelastic.mooney_rivlin, + ... C10=0.25, + ... C01=0.25, + ... ) & fem.Volumetric(bulk=5000) + >>> view = umat.view() + >>> λ = view.biaxial()[0] + >>> + >>> F = np.zeros((3, 3, 1, λ[0].size)) + >>> for a in range(3): + ... F[a, a] = λ[a] + >>> + >>> umat.is_stable([F]) + array([[ True, True, True, True, True, True, True, True, False, + False, False, False, False, False, False, False]]) + + """ + N = eigh(dot(transpose(x[0]), x[0]))[1] + n = eigh(dot(x[0], transpose(x[0])))[1] + + if hessian is None: + hessian = self.hessian(x)[0] + + A = np.einsum( + "iJkL...,ia...,Ja...,kb...,Lb...->ab...", + hessian, + n, + N, + n, + N, + optimize=True, + ) + + stretch_change = np.diagonal(np.linalg.inv(A.T).T, axis1=0, axis2=1) + + return np.all(stretch_change > 0, axis=-1) + def __and__(self, other_material): return CompositeMaterial(self, other_material) @@ -512,12 +614,12 @@ def __init__(self, material, other_material): def gradient(self, x, **kwargs): gradients = [material.gradient(x, **kwargs) for material in self.materials] - nfields = len(x) - 1 + nfields = max(1, len(x) - 1) P = [np.sum([grad[i] for grad in gradients], axis=0) for i in range(nfields)] statevars_new = gradients[0][-1] return [*P, statevars_new] def hessian(self, x, **kwargs): hessians = [material.hessian(x, **kwargs) for material in self.materials] - nfields = len(x) - 1 + nfields = max(1, len(x) - 1) return [np.sum([hess[i] for hess in hessians], axis=0) for i in range(nfields)] diff --git a/src/felupe/constitution/_view.py b/src/felupe/constitution/_view.py index 8e88fba3..69f53390 100644 --- a/src/felupe/constitution/_view.py +++ b/src/felupe/constitution/_view.py @@ -72,14 +72,14 @@ def plot(self, ax=None, show_title=True, show_kwargs=True, **kwargs): shear and biaxial tension.""" import matplotlib.pyplot as plt - + data = self.evaluate() if ax is None: fig, ax = plt.subplots() - for stretch, force, label in data: - ax.plot(stretch, force, label=label, **kwargs) + for stretches, force, label in data: + ax.plot(stretches[0], force, label=label, **kwargs) ax.set_xlabel(r"Stretch $l/L$ $\rightarrow$") ax.set_ylabel("Normal force per undeformed area" + r" $N/A$ $\rightarrow$") @@ -254,7 +254,7 @@ def fun(λ3): warnings.warn("Uniaxial data with volume ratio det(F) <= 0 included.") P[0, 0][~valid] = np.nan - return λ1, P[0, 0].ravel(), "Uniaxial" + return (λ1, λ2, λ3), P[0, 0].ravel(), "Uniaxial" def planar(self, stretches=None, **kwargs): """Normal force per undeformed area vs stretch curve for a planar shear @@ -327,7 +327,7 @@ def fun(λ3): warnings.warn("Planar Shear data with volume ratio det(F) <= 0 included.") P[0, 0][~valid] = np.nan - return λ1, P[0, 0].ravel(), "Planar Shear" + return (λ1, λ2, λ3), P[0, 0].ravel(), "Planar Shear" def biaxial(self, stretches=None, **kwargs): """Normal force per undeformed area vs stretch curve for a equi-biaxial @@ -399,7 +399,7 @@ def fun(λ3): warnings.warn("Biaxial data with volume ratio det(F) <= 0 included.") P[0, 0][~valid] = np.nan - return λ1, P[0, 0].ravel(), "Biaxial" + return (λ1, λ2, λ3), P[0, 0].ravel(), "Biaxial" class ViewMaterialIncompressible(PlotMaterial): @@ -503,7 +503,11 @@ def uniaxial(self, stretches=None): else: P, self.statevars = self.umat.gradient([F, None]) - return λ1, (P[0, 0] - λ3 / λ1 * P[2, 2]).ravel(), "Uniaxial (Incompressible)" + return ( + (λ1, λ2, λ3), + (P[0, 0] - λ3 / λ1 * P[2, 2]).ravel(), + "Uniaxial (Incompressible)", + ) def planar(self, stretches=None): """Normal force per undeformed area vs stretch curve for a planar shear @@ -544,7 +548,7 @@ def planar(self, stretches=None): P, self.statevars = self.umat.gradient([F, None]) return ( - λ1, + (λ1, λ2, λ3), (P[0, 0] - λ3 / λ1 * P[2, 2]).ravel(), "Planar Shear (Incompressible)", ) @@ -586,4 +590,8 @@ def biaxial(self, stretches=None): else: P, self.statevars = self.umat.gradient([F, None]) - return λ1, (P[0, 0] - λ3 / λ1 * P[2, 2]).ravel(), "Biaxial (Incompressible)" + return ( + (λ1, λ2, λ3), + (P[0, 0] - λ3 / λ1 * P[2, 2]).ravel(), + "Biaxial (Incompressible)", + ) diff --git a/tests/test_constitution.py b/tests/test_constitution.py index 5bdf470e..1be92871 100644 --- a/tests/test_constitution.py +++ b/tests/test_constitution.py @@ -101,6 +101,8 @@ def test_nh(): assert P[0].shape == (3, 3, *F[0].shape[-2:]) assert A[0].shape == (3, 3, 3, 3, *F[0].shape[-2:]) + assert np.all(nh.is_stable(F)) + def test_linear(): r, F = pre(sym=False, add_identity=True) From c609486e5ede62a290df1cbe50dbdbf1a04c822c Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Wed, 28 May 2025 16:37:19 +0200 Subject: [PATCH 32/40] Add `mechanics.Truss` and `constitution.LinearElastic1D` (#966) * Add `Truss` and `LinearElastic1D` * Update CHANGELOG.md * Update the API docs * format black * `Truss`: Enforce `E` and `A` as arrays * Update test_mechanics.py * Change `Truss` to `TrussBody` --- CHANGELOG.md | 1 + docs/felupe/constitution/core.rst | 6 + docs/felupe/mechanics.rst | 6 + docs/tutorial/examples/extut02_job.py | 2 +- .../examples/extut03_building_blocks.py | 12 +- examples/ex01_beam.py | 2 +- examples/ex02_plate-with-hole.py | 6 +- examples/ex03_plasticity.py | 2 +- examples/ex04_balloon.py | 10 +- examples/ex05_rubber-metal-bushing.py | 10 +- examples/ex06_rubber-metal-spring.py | 10 +- examples/ex07_engine-mount.py | 10 +- examples/ex08_shear.py | 14 +- examples/ex09_numeric-continuation.py | 50 ++++- examples/ex10_poisson-equation.py | 5 +- examples/ex11_notch-stress.py | 10 +- examples/ex12_foot-bone.py | 4 +- examples/ex13_morph-rubber-wheel.py | 38 ++-- examples/ex14_hyperelasticity.py | 4 +- examples/ex15_hexmesh-metacone.py | 2 +- examples/ex16_deeplearning-torch.py | 8 +- examples/ex17_torsion-gif.py | 4 +- .../ex18_nonlinear-viscoelasticity-newton.py | 4 +- examples/ex19_taylor-hood.py | 2 +- examples/ex20_third-medium-contact.py | 4 +- src/felupe/__init__.py | 4 + src/felupe/constitution/__init__.py | 2 + .../linear_elasticity/__init__.py | 2 + .../linear_elasticity/_linear_elastic_1d.py | 117 ++++++++++ .../tensortrax/_updated_lagrange.py | 2 +- src/felupe/mechanics/__init__.py | 2 + src/felupe/mechanics/_pointload.py | 2 +- src/felupe/mechanics/_truss.py | 205 ++++++++++++++++++ src/felupe/region/_boundary.py | 16 +- tests/test_basis.py | 6 +- tests/test_bilinearform.py | 6 +- tests/test_constitution.py | 6 +- tests/test_constitution_jax.py | 6 +- tests/test_constitution_newton.py | 6 +- tests/test_dof.py | 6 +- tests/test_dtype.py | 6 +- tests/test_element.py | 6 +- tests/test_field.py | 6 +- tests/test_form.py | 6 +- tests/test_free_vibration.py | 6 +- tests/test_job.py | 2 +- tests/test_math.py | 6 +- tests/test_mechanics.py | 27 ++- tests/test_mpc.py | 6 +- tests/test_planestrain.py | 6 +- tests/test_plot.py | 6 +- tests/test_region.py | 6 +- tests/test_solve.py | 6 +- tests/test_tools.py | 6 +- 54 files changed, 551 insertions(+), 164 deletions(-) create mode 100644 src/felupe/constitution/linear_elasticity/_linear_elastic_1d.py create mode 100644 src/felupe/mechanics/_truss.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 59a3c923..91f7972d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,7 @@ All notable changes to this project will be documented in this file. The format - Add `Mesh.clear_points_without_cells()` to clear the list of points without cells (useful for center-points of multi-point constraints). - Release FElupe on conda-forge, starting with v9.2.0. - Add `ConstitutiveMaterial.is_stable()` which returns a boolean mask of stability for isotropic material model formulations. Note that this will require an additional volumetric part of the strain energy density function for hyperelastic material model formulations without a volumetric part. +- Add the linear-elastic material formulation `constitution.LinearElastic1D()` and a truss-body `mechanics.TrussBody()`. ### Changed - Change the required setuptools-version in the build-system table of `pyproject.toml` to match PEP639 (setuptools>=77.0.3). diff --git a/docs/felupe/constitution/core.rst b/docs/felupe/constitution/core.rst index e3246c4b..5c963230 100644 --- a/docs/felupe/constitution/core.rst +++ b/docs/felupe/constitution/core.rst @@ -18,6 +18,7 @@ This page contains the core (hard-coded) constitutive material model formulation .. autosummary:: LinearElastic + LinearElastic1D LinearElasticPlaneStress constitution.LinearElasticPlaneStrain constitution.LinearElasticTensorNotation @@ -65,6 +66,11 @@ This page contains the core (hard-coded) constitutive material model formulation :undoc-members: :inherited-members: +.. autoclass:: felupe.LinearElastic1D + :members: + :undoc-members: + :inherited-members: + .. autofunction:: felupe.linear_elastic .. autoclass:: felupe.LinearElasticLargeStrain diff --git a/docs/felupe/mechanics.rst b/docs/felupe/mechanics.rst index 2f9d4d8e..c8314278 100644 --- a/docs/felupe/mechanics.rst +++ b/docs/felupe/mechanics.rst @@ -14,6 +14,7 @@ Mechanics SolidBodyForce SolidBodyPressure SolidBodyCauchyStress + TrussBody **Steps and Jobs** @@ -73,6 +74,11 @@ Mechanics :undoc-members: :show-inheritance: +.. autoclass:: felupe.TrussBody + :members: + :undoc-members: + :show-inheritance: + .. autoclass:: felupe.Step :members: :undoc-members: diff --git a/docs/tutorial/examples/extut02_job.py b/docs/tutorial/examples/extut02_job.py index dc7f3993..7ce20c38 100644 --- a/docs/tutorial/examples/extut02_job.py +++ b/docs/tutorial/examples/extut02_job.py @@ -6,7 +6,7 @@ and run a **Job**. * create a **Step** with ramped boundary conditions - + * run a **Job** and export a XDMF time-series file This tutorial once again covers the essential high-level parts of creating and solving diff --git a/docs/tutorial/examples/extut03_building_blocks.py b/docs/tutorial/examples/extut03_building_blocks.py index fb5c06da..c242d23f 100644 --- a/docs/tutorial/examples/extut03_building_blocks.py +++ b/docs/tutorial/examples/extut03_building_blocks.py @@ -5,17 +5,17 @@ .. 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`. diff --git a/examples/ex01_beam.py b/examples/ex01_beam.py index a022f0f6..a9e439e4 100644 --- a/examples/ex01_beam.py +++ b/examples/ex01_beam.py @@ -5,7 +5,7 @@ .. topic:: Apply a gravity load on a solid body. * create a solid body and apply the gravity load - + * linear-elastic analysis The displacement due to gravity of a cantilever beam with young's modulus diff --git a/examples/ex02_plate-with-hole.py b/examples/ex02_plate-with-hole.py index 5dd16265..9a45e15b 100644 --- a/examples/ex02_plate-with-hole.py +++ b/examples/ex02_plate-with-hole.py @@ -5,11 +5,11 @@ .. topic:: Plane stress linear analysis. * create and mesh a plate with a hole - + * define a solid body with a linear-elastic plane stress material - + * create an external pressure load - + * export and plot stress results diff --git a/examples/ex03_plasticity.py b/examples/ex03_plasticity.py index c40e0392..1a6149df 100644 --- a/examples/ex03_plasticity.py +++ b/examples/ex03_plasticity.py @@ -12,7 +12,7 @@ The normal plastic strain due to a body force applied on a solid with a linear-elastic plastic material formulation with isotropic hardening with young's modulus -:math:`E=210000`, poisson ratio :math:`\nu=0.3`, isotropic hardening modulus +:math:`E=210000`, poisson ratio :math:`\nu=0.3`, isotropic hardening modulus :math:`K=1000`, yield stress :math:`\sigma_y=355`, length :math:`L=3` and cross section area :math:`A=1` is to be evaluated. diff --git a/examples/ex04_balloon.py b/examples/ex04_balloon.py index f22911fc..cafccaa1 100644 --- a/examples/ex04_balloon.py +++ b/examples/ex04_balloon.py @@ -7,14 +7,14 @@ * use FElupe with contique * plot pressure-displacement curve - + * view the deformed balloon .. admonition:: This example requires external packages. :class: hint - + .. code-block:: - + pip install contique With the help of `contique `_ it is possible to @@ -35,8 +35,8 @@ .. math:: :label: neo-hookean-strain-energy - \psi = \frac{\mu}{2} \left( - \text{tr}(\boldsymbol{C}) - \ln(\det(\boldsymbol{C})) + \psi = \frac{\mu}{2} \left( + \text{tr}(\boldsymbol{C}) - \ln(\det(\boldsymbol{C})) \right) """ diff --git a/examples/ex05_rubber-metal-bushing.py b/examples/ex05_rubber-metal-bushing.py index 100640b3..949c53e0 100644 --- a/examples/ex05_rubber-metal-bushing.py +++ b/examples/ex05_rubber-metal-bushing.py @@ -7,18 +7,18 @@ * create and stack meshes * define a boundary condition with torsional loading - + * work with multiple solid bodies - + * create a step and add it to a job - + * plot strains and stresses .. admonition:: This example requires external packages. :class: hint - + .. code-block:: - + pip install pypardiso An elastic bearing is subjected to combined multiaxial radial-torsional-cardanic diff --git a/examples/ex06_rubber-metal-spring.py b/examples/ex06_rubber-metal-spring.py index 2c759276..415b37b6 100644 --- a/examples/ex06_rubber-metal-spring.py +++ b/examples/ex06_rubber-metal-spring.py @@ -5,18 +5,18 @@ .. topic:: A hyperelastic spring with a simplified frictionless contact. * read a mesh file - + * define an isotropic hyperelastic solid body - + * setup a simplified frictionless elastic-to-rigid contact interaction - + * export and plot the log. strain .. admonition:: This example requires external packages. :class: hint - + .. code-block:: - + pip install pypardiso A `meshed three-dimensional geometry <../_static/ex06_rubber-metal-spring_mesh.vtk>`_ of diff --git a/examples/ex07_engine-mount.py b/examples/ex07_engine-mount.py index 93e97aa6..0a2c4e5b 100644 --- a/examples/ex07_engine-mount.py +++ b/examples/ex07_engine-mount.py @@ -5,11 +5,11 @@ .. topic:: A rubberlike-metal component used as an engine-mount. * read and combine mesh files - + * define an isotropic hyperelastic solid body - + * create consecutive steps and add them to a job - + * export and plot characteristic curves @@ -19,9 +19,9 @@ The air inside the structure is meshed as a hyperelastic solid with no volumetric part of the strain energy function for a simplified treatment of the rubber contact. The metal parts are simplified as rigid bodies. Three mesh files are provided for this -example: +example: -* a `mesh for the metal parts <../_static/ex07_engine-mount_mesh-metal.vtk>`_, +* a `mesh for the metal parts <../_static/ex07_engine-mount_mesh-metal.vtk>`_, * a `mesh for the rubber blocks <../_static/ex07_engine-mount_mesh-rubber.vtk>`_ as well as diff --git a/examples/ex08_shear.py b/examples/ex08_shear.py index 774c6813..f7453e22 100644 --- a/examples/ex08_shear.py +++ b/examples/ex08_shear.py @@ -5,20 +5,20 @@ .. topic:: Plane strain hyperelastic non-homogeneous shear loadcase * define a non-homogeneous shear loadcase - + * use a mixed hyperelastic formulation in plane strain - + * assign a micro-sphere material formulation * define a step and a job along with a callback-function - + * export and visualize principal stretches - + * plot force - displacement curves -Two rubber blocks of height :math:`H` and length :math:`L`, both glued to a -rigid plate on their top and bottom faces, are subjected to a displacement -controlled non-homogeneous shear deformation by :math:`u_{ext}` in combination +Two rubber blocks of height :math:`H` and length :math:`L`, both glued to a +rigid plate on their top and bottom faces, are subjected to a displacement +controlled non-homogeneous shear deformation by :math:`u_{ext}` in combination with a compressive normal force :math:`F`. .. image:: ../../examples/ex08_shear_sketch.svg diff --git a/examples/ex09_numeric-continuation.py b/examples/ex09_numeric-continuation.py index c0491e3c..9d6b36b2 100644 --- a/examples/ex09_numeric-continuation.py +++ b/examples/ex09_numeric-continuation.py @@ -12,9 +12,9 @@ .. topic:: Numeric continuation of a hyperelastic cube. * use FElupe with contique - + * on-the-fly XDMF-file export - + * plot a force-displacement curve """ @@ -25,7 +25,6 @@ import numpy as np import felupe as fem -import felupe.constitution.tensortrax as mat # %% # First, setup a problem as usual (mesh, region, field, boundaries and umat). The @@ -42,10 +41,24 @@ # partition degrees of freedom dof0, dof1 = fem.dof.partition(field, bounds) +import tensortrax.math as tm + # constitutive isotropic hyperelastic material formulation -yeoh = mat.Hyperelastic(mat.models.hyperelastic.yeoh, C10=0.5, C20=-0.25, C30=0.025) +import felupe.constitution.tensortrax as mat + + +def fun(C): + I3 = tm.linalg.det(C) + + return 0.5 * (I3 ** (-1 / 3) * tm.trace(C) - 3) + 10 * (tm.sqrt(I3) - 1) ** 2 / 2 + return 0.5 * (I3 ** (-1 / 3) * tm.trace(C) - 3) + 10 * (tm.sqrt(I3) - 1) ** 2 / 2 + # return 0.5 * (I3**(-1/3) * tm.trace(C) - 3) + 10 * tm.log(tm.sqrt(I3))**2 / 2 + + +yeoh = umat = fem.Hyperelastic(fun) + ax = yeoh.plot(incompressible=True, ux=np.linspace(1, 2.76), bx=None, ps=None) -body = fem.SolidBodyNearlyIncompressible(yeoh, field, bulk=5000) +body = fem.SolidBody(yeoh, field) # %% # An external normal force is applied at :math:`x=1` on a quarter model of a cube with @@ -56,7 +69,7 @@ # 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 +values_load = -np.vstack([v, np.ones_like(v) * 0, np.ones_like(v) * 0]).T load = fem.PointLoad(field, right, values_load) @@ -116,20 +129,21 @@ def step_to_xdmf(step, res): fun=fun, jac=[dfundx, dfundl], x0=field[0][dof1], + control0=(-1, 1), lpf0=0, - dxmax=0.06, - dlpfmax=0.02, - maxsteps=35, + dxmax=0.1, + dlpfmax=0.1, + maxsteps=100, rebalance=True, overshoot=1.05, callback=step_to_xdmf, - tol=1e-2, + tol=1e-8, ) X = np.array([res.x for res in Res]) # check the final lpf value -assert np.isclose(X[-1, 1], -0.3982995) +# assert np.isclose(X[-1, 1], -0.3982995) # %% # Finally, the force-displacement curve is plotted. It can be seen that the resulting @@ -142,3 +156,17 @@ def step_to_xdmf(step, res): plt.ylabel(r"load-proportionality-factor $\lambda$ $\longrightarrow$") field.plot("Displacement", component=0).show() + +umat = yeoh +stability = umat.is_stable(field.extract()) +# stability = umat.is_stable(field.extract(), hessian=solid.evaluate.hessian(field)[0]) +field.view(cell_data={"Stability": stability.sum(axis=0) / 8}).plot( + "Stability", cmap=plt.get_cmap("coolwarm_r", 8), clim=[0, 1] +).show() + +import numpy as np + +stretches = fem.math.eigvalsh(field.extract()[0])[:, 1, 0] + +F = np.diag(stretches) +s = umat.is_stable([F.reshape(3, 3, 1, 1)]) diff --git a/examples/ex10_poisson-equation.py b/examples/ex10_poisson-equation.py index c95db314..1725442c 100644 --- a/examples/ex10_poisson-equation.py +++ b/examples/ex10_poisson-equation.py @@ -7,12 +7,12 @@ .. math:: :label: poisson - + \text{div}(\boldsymbol{\nabla} u) + f = 0 \quad \text{in} \quad \Omega .. math:: :label: poisson-boundaries - + u &= 0 \quad \text{on} \quad \Gamma_u f &= 1 \quad \text{in} \quad \Omega @@ -28,6 +28,7 @@ \ d\Omega = \int_\Omega \delta u \cdot f \ d\Omega """ + # sphinx_gallery_thumbnail_number = -1 import felupe as fem diff --git a/examples/ex11_notch-stress.py b/examples/ex11_notch-stress.py index 31b1cfa9..2f0d7e09 100644 --- a/examples/ex11_notch-stress.py +++ b/examples/ex11_notch-stress.py @@ -5,20 +5,20 @@ .. topic:: Three-dimensional linear-elastic analysis. * create a hexahedron mesh - + * define a linear-elastic solid body - + * project the linear-elastic stress tensor to the mesh-points - + * plot the max. principal stress component * evaluate the fatigue life .. admonition:: This example requires external packages. :class: hint - + .. code-block:: - + pip install pypardiso A linear-elastic notched plate is subjected to uniaxial tension. The cell-based mean of diff --git a/examples/ex12_foot-bone.py b/examples/ex12_foot-bone.py index d58f3572..c4b6034e 100644 --- a/examples/ex12_foot-bone.py +++ b/examples/ex12_foot-bone.py @@ -6,9 +6,9 @@ .. admonition:: This example requires external packages. :class: hint - + .. code-block:: - + pip install pypardiso """ diff --git a/examples/ex13_morph-rubber-wheel.py b/examples/ex13_morph-rubber-wheel.py index af8106b0..9f47c1b2 100644 --- a/examples/ex13_morph-rubber-wheel.py +++ b/examples/ex13_morph-rubber-wheel.py @@ -20,17 +20,17 @@ .. math:: :label: morph-state-ex - + \boldsymbol{C} &= \boldsymbol{F}^T \boldsymbol{F} - + I_3 &= \det (\boldsymbol{C}) - + \hat{\boldsymbol{C}} &= I_3^{-1/3} \boldsymbol{C} - + \hat{\lambda}^2_\alpha &= \text{eigvals}(\hat{\boldsymbol{C}}) - + \hat{C}_T &= \max \left( \hat{\lambda}^2_\alpha - \hat{\lambda}^2_\beta \right) - + \hat{C}_T^S &= \max \left( \hat{C}_T, \hat{C}_{T,n}^S \right) A sigmoid-function is used inside the deformation-dependent variables :math:`\alpha`, @@ -38,13 +38,13 @@ .. math:: :label: morph-sigmoid-ex - + f(x) &= \frac{1}{\sqrt{1 + x^2}} - + \alpha &= p_1 + p_2 \ f(p_3\ C_T^S) - + \beta &= p_4\ f(p_3\ C_T^S) - + \gamma &= p_5\ C_T^S\ \left( 1 - f\left(\frac{C_T^S}{p_6}\right) \right) The rate of deformation is described by the Lagrangian tensor and its Tresca-invariant, @@ -57,17 +57,17 @@ .. math:: :label: morph-rate-of-deformation-ex - - \hat{\boldsymbol{L}} &= \text{sym}\left( - \text{dev}(\boldsymbol{C}^{-1} \Delta\boldsymbol{C}) + + \hat{\boldsymbol{L}} &= \text{sym}\left( + \text{dev}(\boldsymbol{C}^{-1} \Delta\boldsymbol{C}) \right) \hat{\boldsymbol{C}} - + \lambda_{\hat{\boldsymbol{L}}, \alpha} &= \text{eigvals}(\hat{\boldsymbol{L}}) - + \hat{L}_T &= \max \left( \lambda_{\hat{\boldsymbol{L}}, \alpha} - \lambda_{\hat{\boldsymbol{L}}, \beta} \right) - + \Delta\boldsymbol{C} &= \boldsymbol{C} - \boldsymbol{C}_n The additional stresses evolve between the limiting stresses, see Eq. @@ -76,17 +76,17 @@ .. math:: :label: morph-stresses-ex - + \boldsymbol{S}_L &= \left( \gamma \exp \left(p_7 \frac{\hat{\boldsymbol{L}}}{\hat{L}_T} \frac{\hat{C}_T}{\hat{C}_T^S} \right) + p8 \frac{\hat{\boldsymbol{L}}}{\hat{L}_T} \right) \boldsymbol{C}^{-1} - + \boldsymbol{S}_A &= \frac{ \boldsymbol{S}_{A,n} + \beta\ \hat{L}_T\ \boldsymbol{S}_L }{1 + \beta\ \hat{L}_T} - + \boldsymbol{S} &= 2 \alpha\ \text{dev}( \hat{\boldsymbol{C}} ) \boldsymbol{C}^{-1} + \text{dev}\left( \boldsymbol{S}_A\ \boldsymbol{C} \right) \boldsymbol{C}^{-1} diff --git a/examples/ex14_hyperelasticity.py b/examples/ex14_hyperelasticity.py index 5b50990f..8c88fd65 100644 --- a/examples/ex14_hyperelasticity.py +++ b/examples/ex14_hyperelasticity.py @@ -5,9 +5,9 @@ .. topic:: Best-fit of hyperelastic material parameters on given experimental data. * define a strain-energy function for a hyperelastic material formulation - + * best-fit the parameters to uniaxial and biaxial tension experiments - + * plot the standard deviations of the material parameters The :func:`Extended Tube ` material model formulation [1]_ is diff --git a/examples/ex15_hexmesh-metacone.py b/examples/ex15_hexmesh-metacone.py index de1c2709..a3cd337a 100644 --- a/examples/ex15_hexmesh-metacone.py +++ b/examples/ex15_hexmesh-metacone.py @@ -5,7 +5,7 @@ .. topic:: Create a 3d dynamic mesh for a metacone component out of hexahedrons. * apply :ref:`mesh-tools ` - + * create a :class:`~felupe.MeshContainer` for meshes associated to two materials """ diff --git a/examples/ex16_deeplearning-torch.py b/examples/ex16_deeplearning-torch.py index 6fdeacf5..90942884 100644 --- a/examples/ex16_deeplearning-torch.py +++ b/examples/ex16_deeplearning-torch.py @@ -7,9 +7,9 @@ * compression of a hyperelastic axisymmetric cylinder * evaluate displacements located at mesh-points - + * train a neural network on the displacement data - + * obtain displacements from the PyTorch model and plot the log. strains First, an axisymmetric model is created. The displacements are saved after each @@ -17,9 +17,9 @@ .. admonition:: This example requires external packages. :class: hint - + .. code-block:: - + pip install torch """ diff --git a/examples/ex17_torsion-gif.py b/examples/ex17_torsion-gif.py index 47b4b818..de3a5602 100644 --- a/examples/ex17_torsion-gif.py +++ b/examples/ex17_torsion-gif.py @@ -7,9 +7,9 @@ * assign non-homogeneous displacements to a boundary condition * create an animation - + * export a GIF-movie - + * evaluate the reaction moment on a boundary condition This example demonstrates how to create a :class:`~felupe.Boundary` for torsional diff --git a/examples/ex18_nonlinear-viscoelasticity-newton.py b/examples/ex18_nonlinear-viscoelasticity-newton.py index 7b271cd9..f23db54e 100644 --- a/examples/ex18_nonlinear-viscoelasticity-newton.py +++ b/examples/ex18_nonlinear-viscoelasticity-newton.py @@ -2,9 +2,9 @@ Uniaxial loading/unloading of a viscoelastic material (VHB 4910) ---------------------------------------------------------------- This example shows how to implement a constitutive material model for -rubber viscoelastic materials using a strain energy density function coupled +rubber viscoelastic materials using a strain energy density function coupled with an ODE for an internal (state) variable [1]_. The ODE is discretized using a -backward-Euler scheme and the resulting nonlinear algebraic equations for the +backward-Euler scheme and the resulting nonlinear algebraic equations for the internal variable are solved using Newton's method [2]_. """ diff --git a/examples/ex19_taylor-hood.py b/examples/ex19_taylor-hood.py index 84a2746d..b913a86e 100644 --- a/examples/ex19_taylor-hood.py +++ b/examples/ex19_taylor-hood.py @@ -2,7 +2,7 @@ Mixed-field hyperelasticity with quadratic triangles ---------------------------------------------------- A 90° section of a plane-strain circle is subjected to frictionless uniaxial compression -by a vertically moved rigid top plate. A mixed-field formulation is used with quadratic +by a vertically moved rigid top plate. A mixed-field formulation is used with quadratic triangles. """ diff --git a/examples/ex20_third-medium-contact.py b/examples/ex20_third-medium-contact.py index 1d4658ce..505d8f6f 100644 --- a/examples/ex20_third-medium-contact.py +++ b/examples/ex20_third-medium-contact.py @@ -5,9 +5,9 @@ .. topic:: Frictionless contact method simulated by a third medium [1]_. * create a mesh container with multiple meshes - + * define multiple solid bodies and create a top-level field - + * create an animation of the deformed solid This contact method uses a third medium for two solid contact bodies with a a Hessian- diff --git a/src/felupe/__init__.py b/src/felupe/__init__.py index f481bc25..9a609023 100644 --- a/src/felupe/__init__.py +++ b/src/felupe/__init__.py @@ -21,6 +21,7 @@ ConstitutiveMaterial, Laplace, LinearElastic, + LinearElastic1D, LinearElasticLargeStrain, LinearElasticOrthotropic, LinearElasticPlaneStress, @@ -85,6 +86,7 @@ SolidBodyPressure, StateNearlyIncompressible, Step, + TrussBody, ) from .mesh import Circle, Cube, Grid, Mesh, MeshContainer, Point, Rectangle from .quadrature import ( @@ -155,6 +157,7 @@ "AreaChange", "Laplace", "LinearElastic", + "LinearElastic1D", "LinearElasticLargeStrain", "LinearElasticOrthotropic", "LinearElasticPlaneStress", @@ -214,6 +217,7 @@ "SolidBodyPressure", "StateNearlyIncompressible", "Step", + "TrussBody", "MultiPointConstraint", "MultiPointContact", "GaussLegendre", diff --git a/src/felupe/constitution/__init__.py b/src/felupe/constitution/__init__.py index 1907f408..93dc02cd 100644 --- a/src/felupe/constitution/__init__.py +++ b/src/felupe/constitution/__init__.py @@ -6,6 +6,7 @@ from .hyperelasticity import NeoHooke, NeoHookeCompressible, OgdenRoxburgh, Volumetric from .linear_elasticity import ( LinearElastic, + LinearElastic1D, LinearElasticLargeStrain, LinearElasticOrthotropic, LinearElasticPlaneStrain, @@ -27,6 +28,7 @@ "NeoHookeCompressible", "Laplace", "LinearElastic", + "LinearElastic1D", "LinearElasticLargeStrain", "LinearElasticOrthotropic", "LinearElasticPlaneStrain", diff --git a/src/felupe/constitution/linear_elasticity/__init__.py b/src/felupe/constitution/linear_elasticity/__init__.py index 3e0f8ef9..69e4b0b1 100644 --- a/src/felupe/constitution/linear_elasticity/__init__.py +++ b/src/felupe/constitution/linear_elasticity/__init__.py @@ -5,6 +5,7 @@ LinearElasticPlaneStress, LinearElasticTensorNotation, ) +from ._linear_elastic_1d import LinearElastic1D from ._linear_elastic_large_strain import LinearElasticLargeStrain from ._linear_elastic_orthotropic import LinearElasticOrthotropic @@ -12,6 +13,7 @@ "lame_converter", "lame_converter_orthotropic", "LinearElastic", + "LinearElastic1D", "LinearElasticLargeStrain", "LinearElasticOrthotropic", "LinearElasticPlaneStrain", diff --git a/src/felupe/constitution/linear_elasticity/_linear_elastic_1d.py b/src/felupe/constitution/linear_elasticity/_linear_elastic_1d.py new file mode 100644 index 00000000..9592b45e --- /dev/null +++ b/src/felupe/constitution/linear_elasticity/_linear_elastic_1d.py @@ -0,0 +1,117 @@ +# -*- coding: utf-8 -*- +""" +This file is part of FElupe. + +FElupe is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +FElupe is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with FElupe. If not, see . +""" + +import numpy as np + +from .._base import ConstitutiveMaterial + + +class LinearElastic1D(ConstitutiveMaterial): + r"""Isotropic one-dimensional linear-elastic material formulation. + + Parameters + ---------- + E : float + Young's modulus. + + Notes + ----- + The stress-stretch relation of the linear-elastic material formulation is given in + Eq. :eq:`linear-elastic-1d` + + .. math:: + :label: linear-elastic-1d + + \sigma = E\ \left(\lambda - 1 \right) + + with the stretch from Eq. :eq:`linear-elastic-strain-1d`. + + .. math:: + :label: linear-elastic-strain-1d + + \lambda = \frac{l}{L} + + + Examples + -------- + .. plot:: + + >>> import felupe as fem + >>> + >>> umat = fem.LinearElastic1D(E=1) + >>> ax = umat.plot() + + """ + + def __init__(self, E): + self.E = np.array(E) + + self.kwargs = {"E": self.E} + + # aliases for gradient and hessian + self.stress = self.gradient + self.elasticity = self.hessian + + # initial variables for calling + # ``self.gradient(self.x)`` and ``self.hessian(self.x)`` + self.x = [np.ones(1), np.zeros(0)] + + def gradient(self, x, out=None): + r"""Evaluate the stress (as a function of the stretch). + + Parameters + ---------- + x : list of ndarray + List with the stretch :math:`\lambda` as first item. + out : ndarray or None, optional + A location into which the result is stored (default is None). + Not implemented. + + Returns + ------- + ndarray + Stress + + """ + + λ, statevars = x[0], x[-1] + + return [self.E * (λ - 1), statevars] + + def hessian(self, x=None, shape=(1,), dtype=None, out=None): + r"""Evaluate the elasticity. The stretch is only used for the shape of the + trailing axes. + + Parameters + ---------- + x : list of ndarray, optional + List with stretch :math:`\lambda` as first item (default is None). + shape : tuple of int, optional + Tuple with shape of the trailing axes (default is (1,)). + out : ndarray or None, optional + A location into which the result is stored (default is None). + Not implemented. + + Returns + ------- + ndarray + Elasticity + + """ + + return [self.E] diff --git a/src/felupe/constitution/tensortrax/_updated_lagrange.py b/src/felupe/constitution/tensortrax/_updated_lagrange.py index 2325b053..1a44ce3e 100644 --- a/src/felupe/constitution/tensortrax/_updated_lagrange.py +++ b/src/felupe/constitution/tensortrax/_updated_lagrange.py @@ -34,7 +34,7 @@ def updated_lagrange(material): Examples -------- .. plot:: - + >>> import felupe as fem >>> import felupe.constitution.tensortrax as mat >>> import tensortrax.math as tm diff --git a/src/felupe/mechanics/__init__.py b/src/felupe/mechanics/__init__.py index 69266ada..46c7ab4b 100644 --- a/src/felupe/mechanics/__init__.py +++ b/src/felupe/mechanics/__init__.py @@ -12,6 +12,7 @@ from ._solidbody_incompressible import SolidBodyNearlyIncompressible from ._solidbody_pressure import SolidBodyPressure from ._step import Step +from ._truss import TrussBody __all__ = [ "Assemble", @@ -32,4 +33,5 @@ "Step", "MultiPointConstraint", "MultiPointContact", + "TrussBody", ] diff --git a/src/felupe/mechanics/_pointload.py b/src/felupe/mechanics/_pointload.py index 389a60df..23c5e7d0 100644 --- a/src/felupe/mechanics/_pointload.py +++ b/src/felupe/mechanics/_pointload.py @@ -134,7 +134,7 @@ def plot( skip[values.shape[1] :] = True if values.shape[1] > 1: - skip[:values.shape[1]][np.isclose(values, 0).all(axis=0)] = True + skip[: values.shape[1]][np.isclose(values, 0).all(axis=0)] = True for a, (skip_axis, direction) in enumerate(zip(skip, np.eye(3))): d = np.broadcast_to(direction.reshape(1, 3), points[self.points].shape) diff --git a/src/felupe/mechanics/_truss.py b/src/felupe/mechanics/_truss.py new file mode 100644 index 00000000..9c0ae507 --- /dev/null +++ b/src/felupe/mechanics/_truss.py @@ -0,0 +1,205 @@ +# -*- coding: utf-8 -*- +""" +This file is part of FElupe. + +FElupe is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +FElupe is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with FElupe. If not, see . +""" + + +import inspect + +import numpy as np + +from ..assembly import IntegralForm +from ..math import dot, dya, identity, sqrt +from ._helpers import Assemble, Evaluate, Results +from ._solidbody import Solid + + +class TrussBody(Solid): + r"""A Truss body with methods for the assembly of sparse vectors/matrices. + + Parameters + ---------- + umat : class + A class which provides methods for evaluating the gradient and the hessian of + the strain energy density function per unit undeformed volume. The function + signatures must be ``dψdλ, ζ_new = umat.gradient([λ, ζ])`` for + the gradient and ``d2ψdλdλ = umat.hessian([λ, ζ])`` for the hessian of + the strain energy density function :math:`\psi(\lambda)`, where + :math:`\lambda` is the stretch and :math:`\zeta` + holds the array of internal state variables. + field : FieldContainer + A field container with one or more fields. + area : float or ndarray + The cross-sectional areas of the trusses. + statevars : ndarray or None, optional + Array of initial internal state variables (default is None). + + Examples + -------- + .. pyvista-plot:: + :force_static: + + >>> import felupe as fem + >>> + >>> mesh = fem.Mesh( + ... points=[[0, 0], [1, 1], [2.0, 0]], cells=[[0, 1], [1, 2]], cell_type="line" + ... ) + >>> region = fem.Region(mesh, fem.Line(), fem.GaussLobatto(order=0, dim=1), grad=False) + >>> field = fem.Field(region, dim=2).as_container() + >>> boundaries = fem.BoundaryDict(fixed=fem.Boundary(field[0], fy=0)) + >>> + >>> umat = fem.LinearElastic1D(E=[1, 1]) + >>> truss = fem.TrussBody(umat, field, area=[1, 1]) + >>> load = fem.PointLoad(field, [1]) + >>> + >>> mesh.plot(plotter=load.plot(plotter=boundaries.plot()), line_width=5).show() + >>> + >>> move = fem.math.linsteps([0, -0.1], num=5, axis=1, axes=2) + >>> step = fem.Step(items=[truss, load], ramp={load: move}, boundaries=boundaries) + >>> job = fem.Job(steps=[step]).evaluate() + + See Also + -------- + felupe.SolidBodyNearlyIncompressible : A (nearly) incompressible solid body with + methods for the assembly of sparse vectors/matrices. + """ + + def __init__(self, umat, field, area, statevars=None): + self.field = field + self.assemble = Assemble(vector=self._vector, matrix=self._matrix) + self.evaluate = Evaluate(gradient=self._gradient, hessian=self._hessian) + self.results = Results(stress=True, elasticity=True) + + self.umat = umat + self.area = np.array(area) + + self.form_vector = IntegralForm( + [None], + v=self.field, + dV=None, + grad_v=[False], + ) + self.form_matrix = IntegralForm( + [None], + v=self.field, + dV=None, + u=self.field, + grad_v=[False], + grad_u=[False], + ) + + if statevars is not None: + self.results.statevars = statevars + else: + statevars_shape = (0,) + if hasattr(umat, "x"): + statevars_shape = umat.x[-1].shape + self.results.statevars = np.zeros( + ( + *statevars_shape, + field.region.mesh.ncells, + ) + ) + + cells = self.field.region.mesh.cells + self.X = self.field.region.mesh.points[cells].T + dX = self.X[:, 1] - self.X[:, 0] + self.length_undeformed = sqrt(dot(dX, dX, mode=(1, 1))) + + def _kinematics(self, field): + + cells = self.field.region.mesh.cells + u = field[0].values[cells].T + x = self.X + u + + dx = x[:, 1] - x[:, 0] + + length_deformed = sqrt(dot(dx, dx, mode=(1, 1))) + + stretch = length_deformed / self.length_undeformed + normal_deformed = dx / length_deformed + + return stretch, length_deformed, normal_deformed + + def _gradient(self, field=None, args=(), kwargs=None): + if kwargs is None: + kwargs = {} + + if field is not None: + self.field = field + self.results.kinematics = self._kinematics(field) + + if "out" in inspect.signature(self.umat.gradient).parameters: + kwargs["out"] = self.results.gradient + + gradient = self.umat.gradient( + [*self.results.kinematics, self.results.statevars], *args, **kwargs + ) + self.results.gradient = self.results.stress = gradient[0] + self.results._statevars = gradient[-1] + + return self.results.gradient + + def _hessian(self, field=None, args=(), kwargs=None): + if kwargs is None: + kwargs = {} + + if field is not None: + self.field = field + self.results.kinematics = self._kinematics(field) + + if "out" in inspect.signature(self.umat.hessian).parameters: + kwargs["out"] = self.results.hessian + + hessian = self.umat.hessian( + [*self.results.kinematics, self.results.statevars], *args, **kwargs + ) + self.results.hessian = self.results.elasticity = hessian[0] + + return self.results.hessian + + def _vector(self, field=None, parallel=None, args=(), kwargs=None): + + self.results.stress = self._gradient(field, args=args, kwargs=kwargs) + normal_deformed = self.results.kinematics[2] + + force = self.results.stress * self.area + r = force * np.array([-normal_deformed, normal_deformed]) # (a, i, cell) + + self.results.force = self.form_vector.assemble(values=[r], parallel=parallel) + + return self.results.force + + def _matrix(self, field=None, parallel=None, args=(), kwargs=None): + + self.results.hessian = self._hessian(field, args=args, kwargs=kwargs) + L = self.length_undeformed + stretch, l, n = self.results.kinematics + + S = self.results.stress = self._gradient(field, args=args, kwargs=kwargs) + dSdE = self.results.elasticity = self._hessian(field, args=args, kwargs=kwargs) + + m = dya(n, n, mode=1) + eye = identity(dim=len(n), shape=(1,)) + K_EE = dSdE / L * self.area * m + S / l * self.area * (eye - m) + + K = np.stack([[K_EE, -K_EE], [-K_EE, K_EE]]) # (a, b, i, j, cell) + + self.results.stiffness = self.form_matrix.assemble( + values=[K.transpose([0, 2, 1, 3, 4])], parallel=parallel + ) + + return self.results.stiffness diff --git a/src/felupe/region/_boundary.py b/src/felupe/region/_boundary.py index 9d26968b..77f7da60 100644 --- a/src/felupe/region/_boundary.py +++ b/src/felupe/region/_boundary.py @@ -203,9 +203,7 @@ def boundary_cells_hexahedron20(mesh): def boundary_cells_hexahedron27(mesh): "Convert the cells array of a bi-quadratic hex mesh into a boundary cells array." - cells_hexahedron20, cells_faces_hexahedron20 = boundary_cells_hexahedron20( - mesh - ) + cells_hexahedron20, cells_faces_hexahedron20 = boundary_cells_hexahedron20(mesh) # midpoints of faces (boundary) of a hexahedron i = [20, 21, 22, 23, 24, 25] @@ -593,9 +591,7 @@ def _init_faces(self): tangents.append(dX_1 / np.linalg.norm(dX_1, axis=0)) if self.ensure_3d: - tangents[0] = np.insert( - tangents[0], len(tangents[0]), 0, axis=0 - ) + tangents[0] = np.insert(tangents[0], len(tangents[0]), 0, axis=0) other_tangent = np.zeros_like(tangents[0]) other_tangent[2] = 1.0 tangents.append(other_tangent) @@ -607,12 +603,8 @@ def _init_faces(self): ): dA_1 = cross(self.dXdr[:, 0], self.dXdr[:, 1]) - tangents.append( - -self.dXdr[:, 0] / np.linalg.norm(self.dXdr[:, 0], axis=0) - ) - tangents.append( - self.dXdr[:, 1] / np.linalg.norm(self.dXdr[:, 1], axis=0) - ) + tangents.append(-self.dXdr[:, 0] / np.linalg.norm(self.dXdr[:, 0], axis=0)) + tangents.append(self.dXdr[:, 1] / np.linalg.norm(self.dXdr[:, 1], axis=0)) dA = -dA_1 * self.quadrature.weights.reshape(-1, 1) diff --git a/tests/test_basis.py b/tests/test_basis.py index 22f93c78..73bf5446 100644 --- a/tests/test_basis.py +++ b/tests/test_basis.py @@ -1,11 +1,11 @@ # -*- coding: utf-8 -*- """ - _______ _______ ___ __ __ _______ _______ + _______ _______ ___ __ __ _______ _______ | || || | | | | || || | | ___|| ___|| | | | | || _ || ___| -| |___ | |___ | | | |_| || |_| || |___ +| |___ | |___ | | | |_| || |_| || |___ | ___|| ___|| |___ | || ___|| ___| -| | | |___ | || || | | |___ +| | | |___ | || || | | |___ |___| |_______||_______||_______||___| |_______| This file is part of felupe. diff --git a/tests/test_bilinearform.py b/tests/test_bilinearform.py index a3173861..15138f5f 100644 --- a/tests/test_bilinearform.py +++ b/tests/test_bilinearform.py @@ -1,11 +1,11 @@ # -*- coding: utf-8 -*- """ - _______ _______ ___ __ __ _______ _______ + _______ _______ ___ __ __ _______ _______ | || || | | | | || || | | ___|| ___|| | | | | || _ || ___| -| |___ | |___ | | | |_| || |_| || |___ +| |___ | |___ | | | |_| || |_| || |___ | ___|| ___|| |___ | || ___|| ___| -| | | |___ | || || | | |___ +| | | |___ | || || | | |___ |___| |_______||_______||_______||___| |_______| This file is part of felupe. diff --git a/tests/test_constitution.py b/tests/test_constitution.py index 1be92871..41036975 100644 --- a/tests/test_constitution.py +++ b/tests/test_constitution.py @@ -1,11 +1,11 @@ # -*- coding: utf-8 -*- """ - _______ _______ ___ __ __ _______ _______ + _______ _______ ___ __ __ _______ _______ | || || | | | | || || | | ___|| ___|| | | | | || _ || ___| -| |___ | |___ | | | |_| || |_| || |___ +| |___ | |___ | | | |_| || |_| || |___ | ___|| ___|| |___ | || ___|| ___| -| | | |___ | || || | | |___ +| | | |___ | || || | | |___ |___| |_______||_______||_______||___| |_______| This file is part of felupe. diff --git a/tests/test_constitution_jax.py b/tests/test_constitution_jax.py index 483156ad..6cdf7e3e 100644 --- a/tests/test_constitution_jax.py +++ b/tests/test_constitution_jax.py @@ -1,11 +1,11 @@ # -*- coding: utf-8 -*- """ - _______ _______ ___ __ __ _______ _______ + _______ _______ ___ __ __ _______ _______ | || || | | | | || || | | ___|| ___|| | | | | || _ || ___| -| |___ | |___ | | | |_| || |_| || |___ +| |___ | |___ | | | |_| || |_| || |___ | ___|| ___|| |___ | || ___|| ___| -| | | |___ | || || | | |___ +| | | |___ | || || | | |___ |___| |_______||_______||_______||___| |_______| This file is part of felupe. diff --git a/tests/test_constitution_newton.py b/tests/test_constitution_newton.py index 181c7e40..74a39884 100644 --- a/tests/test_constitution_newton.py +++ b/tests/test_constitution_newton.py @@ -1,11 +1,11 @@ # -*- coding: utf-8 -*- """ - _______ _______ ___ __ __ _______ _______ + _______ _______ ___ __ __ _______ _______ | || || | | | | || || | | ___|| ___|| | | | | || _ || ___| -| |___ | |___ | | | |_| || |_| || |___ +| |___ | |___ | | | |_| || |_| || |___ | ___|| ___|| |___ | || ___|| ___| -| | | |___ | || || | | |___ +| | | |___ | || || | | |___ |___| |_______||_______||_______||___| |_______| This file is part of felupe. diff --git a/tests/test_dof.py b/tests/test_dof.py index ba61ee7d..51745439 100644 --- a/tests/test_dof.py +++ b/tests/test_dof.py @@ -1,11 +1,11 @@ # -*- coding: utf-8 -*- """ - _______ _______ ___ __ __ _______ _______ + _______ _______ ___ __ __ _______ _______ | || || | | | | || || | | ___|| ___|| | | | | || _ || ___| -| |___ | |___ | | | |_| || |_| || |___ +| |___ | |___ | | | |_| || |_| || |___ | ___|| ___|| |___ | || ___|| ___| -| | | |___ | || || | | |___ +| | | |___ | || || | | |___ |___| |_______||_______||_______||___| |_______| This file is part of felupe. diff --git a/tests/test_dtype.py b/tests/test_dtype.py index 08ff2862..22950393 100644 --- a/tests/test_dtype.py +++ b/tests/test_dtype.py @@ -1,11 +1,11 @@ # -*- coding: utf-8 -*- """ - _______ _______ ___ __ __ _______ _______ + _______ _______ ___ __ __ _______ _______ | || || | | | | || || | | ___|| ___|| | | | | || _ || ___| -| |___ | |___ | | | |_| || |_| || |___ +| |___ | |___ | | | |_| || |_| || |___ | ___|| ___|| |___ | || ___|| ___| -| | | |___ | || || | | |___ +| | | |___ | || || | | |___ |___| |_______||_______||_______||___| |_______| This file is part of felupe. diff --git a/tests/test_element.py b/tests/test_element.py index afc92a55..27cd46ed 100644 --- a/tests/test_element.py +++ b/tests/test_element.py @@ -1,11 +1,11 @@ # -*- coding: utf-8 -*- """ - _______ _______ ___ __ __ _______ _______ + _______ _______ ___ __ __ _______ _______ | || || | | | | || || | | ___|| ___|| | | | | || _ || ___| -| |___ | |___ | | | |_| || |_| || |___ +| |___ | |___ | | | |_| || |_| || |___ | ___|| ___|| |___ | || ___|| ___| -| | | |___ | || || | | |___ +| | | |___ | || || | | |___ |___| |_______||_______||_______||___| |_______| This file is part of felupe. diff --git a/tests/test_field.py b/tests/test_field.py index a3d8f2be..5b0dca23 100644 --- a/tests/test_field.py +++ b/tests/test_field.py @@ -1,11 +1,11 @@ # -*- coding: utf-8 -*- """ - _______ _______ ___ __ __ _______ _______ + _______ _______ ___ __ __ _______ _______ | || || | | | | || || | | ___|| ___|| | | | | || _ || ___| -| |___ | |___ | | | |_| || |_| || |___ +| |___ | |___ | | | |_| || |_| || |___ | ___|| ___|| |___ | || ___|| ___| -| | | |___ | || || | | |___ +| | | |___ | || || | | |___ |___| |_______||_______||_______||___| |_______| This file is part of felupe. diff --git a/tests/test_form.py b/tests/test_form.py index 108fc31f..71fccc3c 100644 --- a/tests/test_form.py +++ b/tests/test_form.py @@ -1,11 +1,11 @@ # -*- coding: utf-8 -*- """ - _______ _______ ___ __ __ _______ _______ + _______ _______ ___ __ __ _______ _______ | || || | | | | || || | | ___|| ___|| | | | | || _ || ___| -| |___ | |___ | | | |_| || |_| || |___ +| |___ | |___ | | | |_| || |_| || |___ | ___|| ___|| |___ | || ___|| ___| -| | | |___ | || || | | |___ +| | | |___ | || || | | |___ |___| |_______||_______||_______||___| |_______| This file is part of felupe. diff --git a/tests/test_free_vibration.py b/tests/test_free_vibration.py index 54f69329..05cdd346 100644 --- a/tests/test_free_vibration.py +++ b/tests/test_free_vibration.py @@ -1,11 +1,11 @@ # -*- coding: utf-8 -*- """ - _______ _______ ___ __ __ _______ _______ + _______ _______ ___ __ __ _______ _______ | || || | | | | || || | | ___|| ___|| | | | | || _ || ___| -| |___ | |___ | | | |_| || |_| || |___ +| |___ | |___ | | | |_| || |_| || |___ | ___|| ___|| |___ | || ___|| ___| -| | | |___ | || || | | |___ +| | | |___ | || || | | |___ |___| |_______||_______||_______||___| |_______| This file is part of felupe. diff --git a/tests/test_job.py b/tests/test_job.py index 0ce0a992..8129d0dd 100644 --- a/tests/test_job.py +++ b/tests/test_job.py @@ -57,7 +57,7 @@ def test_job(): def test_job_xdmf(): with pytest.warns(): # solidbodygravity is deprecated field, step = pre(SolidBodyForce=fem.SolidBodyGravity) - + job = fem.Job(steps=[step]) job.evaluate() diff --git a/tests/test_math.py b/tests/test_math.py index 281201d7..ec5b4e25 100644 --- a/tests/test_math.py +++ b/tests/test_math.py @@ -1,11 +1,11 @@ # -*- coding: utf-8 -*- """ - _______ _______ ___ __ __ _______ _______ + _______ _______ ___ __ __ _______ _______ | || || | | | | || || | | ___|| ___|| | | | | || _ || ___| -| |___ | |___ | | | |_| || |_| || |___ +| |___ | |___ | | | |_| || |_| || |___ | ___|| ___|| |___ | || ___|| ___| -| | | |___ | || || | | |___ +| | | |___ | || || | | |___ |___| |_______||_______||_______||___| |_______| This file is part of felupe. diff --git a/tests/test_mechanics.py b/tests/test_mechanics.py index faeac28f..91954b91 100644 --- a/tests/test_mechanics.py +++ b/tests/test_mechanics.py @@ -1,11 +1,11 @@ # -*- coding: utf-8 -*- """ - _______ _______ ___ __ __ _______ _______ + _______ _______ ___ __ __ _______ _______ | || || | | | | || || | | ___|| ___|| | | | | || _ || ___| -| |___ | |___ | | | |_| || |_| || |___ +| |___ | |___ | | | |_| || |_| || |___ | ___|| ___|| |___ | || ___|| ___| -| | | |___ | || || | | |___ +| | | |___ | || || | | |___ |___| |_______||_______||_______||___| |_______| This file is part of felupe. @@ -604,6 +604,26 @@ def test_solidbody_cauchy_stress(): assert np.isclose(field[0].values.max(), 0.971866) +def test_truss(): + + mesh = fem.Mesh( + points=[[0, 0], [1, 1], [2.0, 0]], cells=[[0, 1], [1, 2]], cell_type="line" + ) + region = fem.Region(mesh, fem.Line(), fem.GaussLobatto(order=0, dim=1), grad=False) + field = fem.Field(region, dim=2).as_container() + boundaries = fem.BoundaryDict(fixed=fem.Boundary(field[0], fy=0)) + + umat = fem.LinearElastic1D(E=np.ones(2)) + truss = fem.TrussBody(umat, field, area=np.ones(2)) + load = fem.PointLoad(field, [1]) + + move = fem.math.linsteps([0, -0.1], num=5, axis=1, axes=2) + step = fem.Step(items=[truss, load], ramp={load: move}, boundaries=boundaries) + fem.Job(steps=[step]).evaluate() + + assert np.isclose(field[0].values[1, 1], -0.16302376) + + if __name__ == "__main__": test_simple() test_solidbody() @@ -617,3 +637,4 @@ def test_solidbody_cauchy_stress(): test_load() test_view() test_threefield() + test_truss diff --git a/tests/test_mpc.py b/tests/test_mpc.py index 5e172105..ce02d3ad 100644 --- a/tests/test_mpc.py +++ b/tests/test_mpc.py @@ -1,11 +1,11 @@ # -*- coding: utf-8 -*- """ - _______ _______ ___ __ __ _______ _______ + _______ _______ ___ __ __ _______ _______ | || || | | | | || || | | ___|| ___|| | | | | || _ || ___| -| |___ | |___ | | | |_| || |_| || |___ +| |___ | |___ | | | |_| || |_| || |___ | ___|| ___|| |___ | || ___|| ___| -| | | |___ | || || | | |___ +| | | |___ | || || | | |___ |___| |_______||_______||_______||___| |_______| This file is part of felupe. diff --git a/tests/test_planestrain.py b/tests/test_planestrain.py index fca188a0..bda928bb 100644 --- a/tests/test_planestrain.py +++ b/tests/test_planestrain.py @@ -1,11 +1,11 @@ # -*- coding: utf-8 -*- """ - _______ _______ ___ __ __ _______ _______ + _______ _______ ___ __ __ _______ _______ | || || | | | | || || | | ___|| ___|| | | | | || _ || ___| -| |___ | |___ | | | |_| || |_| || |___ +| |___ | |___ | | | |_| || |_| || |___ | ___|| ___|| |___ | || ___|| ___| -| | | |___ | || || | | |___ +| | | |___ | || || | | |___ |___| |_______||_______||_______||___| |_______| This file is part of felupe. diff --git a/tests/test_plot.py b/tests/test_plot.py index 8a75bca5..3989c285 100644 --- a/tests/test_plot.py +++ b/tests/test_plot.py @@ -1,11 +1,11 @@ # -*- coding: utf-8 -*- """ - _______ _______ ___ __ __ _______ _______ + _______ _______ ___ __ __ _______ _______ | || || | | | | || || | | ___|| ___|| | | | | || _ || ___| -| |___ | |___ | | | |_| || |_| || |___ +| |___ | |___ | | | |_| || |_| || |___ | ___|| ___|| |___ | || ___|| ___| -| | | |___ | || || | | |___ +| | | |___ | || || | | |___ |___| |_______||_______||_______||___| |_______| This file is part of felupe. diff --git a/tests/test_region.py b/tests/test_region.py index 47bdb0a9..cd1455d5 100644 --- a/tests/test_region.py +++ b/tests/test_region.py @@ -1,11 +1,11 @@ # -*- coding: utf-8 -*- """ - _______ _______ ___ __ __ _______ _______ + _______ _______ ___ __ __ _______ _______ | || || | | | | || || | | ___|| ___|| | | | | || _ || ___| -| |___ | |___ | | | |_| || |_| || |___ +| |___ | |___ | | | |_| || |_| || |___ | ___|| ___|| |___ | || ___|| ___| -| | | |___ | || || | | |___ +| | | |___ | || || | | |___ |___| |_______||_______||_______||___| |_______| This file is part of felupe. diff --git a/tests/test_solve.py b/tests/test_solve.py index e5bb7a2b..269117c9 100644 --- a/tests/test_solve.py +++ b/tests/test_solve.py @@ -1,11 +1,11 @@ # -*- coding: utf-8 -*- """ - _______ _______ ___ __ __ _______ _______ + _______ _______ ___ __ __ _______ _______ | || || | | | | || || | | ___|| ___|| | | | | || _ || ___| -| |___ | |___ | | | |_| || |_| || |___ +| |___ | |___ | | | |_| || |_| || |___ | ___|| ___|| |___ | || ___|| ___| -| | | |___ | || || | | |___ +| | | |___ | || || | | |___ |___| |_______||_______||_______||___| |_______| This file is part of felupe. diff --git a/tests/test_tools.py b/tests/test_tools.py index 81aafedc..25475bde 100644 --- a/tests/test_tools.py +++ b/tests/test_tools.py @@ -1,11 +1,11 @@ # -*- coding: utf-8 -*- """ - _______ _______ ___ __ __ _______ _______ + _______ _______ ___ __ __ _______ _______ | || || | | | | || || | | ___|| ___|| | | | | || _ || ___| -| |___ | |___ | | | |_| || |_| || |___ +| |___ | |___ | | | |_| || |_| || |___ | ___|| ___|| |___ | || ___|| ___| -| | | |___ | || || | | |___ +| | | |___ | || || | | |___ |___| |_______||_______||_______||___| |_______| This file is part of felupe. From 7d5f59005ec3f32fbe0ad0f5b2197c82668d6cfb Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Wed, 28 May 2025 16:39:00 +0200 Subject: [PATCH 33/40] Update test_mechanics.py --- tests/test_mechanics.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/tests/test_mechanics.py b/tests/test_mechanics.py index 91954b91..1fd9fbb5 100644 --- a/tests/test_mechanics.py +++ b/tests/test_mechanics.py @@ -623,6 +623,9 @@ def test_truss(): assert np.isclose(field[0].values[1, 1], -0.16302376) + hess = truss.evaluate.hessian(field) + fem.TrussBody(umat, field, area=[1, 1], statevars=[0, 0]) + if __name__ == "__main__": test_simple() From 8f9bdcd36e8b572d5c44631b06d3ed1d54dc5b88 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Wed, 28 May 2025 22:43:40 +0200 Subject: [PATCH 34/40] Revert changes of example 09 --- examples/ex09_numeric-continuation.py | 46 ++++++--------------------- 1 file changed, 9 insertions(+), 37 deletions(-) diff --git a/examples/ex09_numeric-continuation.py b/examples/ex09_numeric-continuation.py index 9d6b36b2..45de3723 100644 --- a/examples/ex09_numeric-continuation.py +++ b/examples/ex09_numeric-continuation.py @@ -25,6 +25,7 @@ import numpy as np import felupe as fem +import felupe.constitution.tensortrax as mat # %% # First, setup a problem as usual (mesh, region, field, boundaries and umat). The @@ -41,24 +42,10 @@ # partition degrees of freedom dof0, dof1 = fem.dof.partition(field, bounds) -import tensortrax.math as tm - # constitutive isotropic hyperelastic material formulation -import felupe.constitution.tensortrax as mat - - -def fun(C): - I3 = tm.linalg.det(C) - - return 0.5 * (I3 ** (-1 / 3) * tm.trace(C) - 3) + 10 * (tm.sqrt(I3) - 1) ** 2 / 2 - return 0.5 * (I3 ** (-1 / 3) * tm.trace(C) - 3) + 10 * (tm.sqrt(I3) - 1) ** 2 / 2 - # return 0.5 * (I3**(-1/3) * tm.trace(C) - 3) + 10 * tm.log(tm.sqrt(I3))**2 / 2 - - -yeoh = umat = fem.Hyperelastic(fun) - +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) -body = fem.SolidBody(yeoh, field) +body = fem.SolidBodyNearlyIncompressible(yeoh, field, bulk=5000) # %% # An external normal force is applied at :math:`x=1` on a quarter model of a cube with @@ -69,7 +56,7 @@ def fun(C): # 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.ones_like(v) * 0, np.ones_like(v) * 0]).T +values_load = np.vstack([v, np.zeros_like(v), np.zeros_like(v)]).T load = fem.PointLoad(field, right, values_load) @@ -129,21 +116,20 @@ def step_to_xdmf(step, res): fun=fun, jac=[dfundx, dfundl], x0=field[0][dof1], - control0=(-1, 1), lpf0=0, - dxmax=0.1, - dlpfmax=0.1, - maxsteps=100, + dxmax=0.06, + dlpfmax=0.02, + maxsteps=35, rebalance=True, overshoot=1.05, callback=step_to_xdmf, - tol=1e-8, + 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) +assert np.isclose(X[-1, 1], -0.3982995) # %% # Finally, the force-displacement curve is plotted. It can be seen that the resulting @@ -156,17 +142,3 @@ def step_to_xdmf(step, res): plt.ylabel(r"load-proportionality-factor $\lambda$ $\longrightarrow$") field.plot("Displacement", component=0).show() - -umat = yeoh -stability = umat.is_stable(field.extract()) -# stability = umat.is_stable(field.extract(), hessian=solid.evaluate.hessian(field)[0]) -field.view(cell_data={"Stability": stability.sum(axis=0) / 8}).plot( - "Stability", cmap=plt.get_cmap("coolwarm_r", 8), clim=[0, 1] -).show() - -import numpy as np - -stretches = fem.math.eigvalsh(field.extract()[0])[:, 1, 0] - -F = np.diag(stretches) -s = umat.is_stable([F.reshape(3, 3, 1, 1)]) From 5ea25a2c1a73f778420fe4d02c9ab6a29b274c2b Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Wed, 28 May 2025 23:13:09 +0200 Subject: [PATCH 35/40] Add `RegionTruss` (#967) --- CHANGELOG.md | 1 + docs/felupe/region.rst | 8 +++++++- src/felupe/__init__.py | 2 ++ src/felupe/mechanics/_truss.py | 25 ++++++++++++++----------- src/felupe/region/__init__.py | 2 ++ src/felupe/region/_templates.py | 13 ++++++++++++- tests/test_mechanics.py | 2 +- 7 files changed, 39 insertions(+), 14 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 91f7972d..45f91ff1 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,7 @@ All notable changes to this project will be documented in this file. The format - Release FElupe on conda-forge, starting with v9.2.0. - Add `ConstitutiveMaterial.is_stable()` which returns a boolean mask of stability for isotropic material model formulations. Note that this will require an additional volumetric part of the strain energy density function for hyperelastic material model formulations without a volumetric part. - Add the linear-elastic material formulation `constitution.LinearElastic1D()` and a truss-body `mechanics.TrussBody()`. +- Add `RegionTruss` with a truss element. This is a line element with a `GaussLobatto(order=0)` quadrature, i.e. with two quadrature points located at the boundaries. ### Changed - Change the required setuptools-version in the build-system table of `pyproject.toml` to match PEP639 (setuptools>=77.0.3). diff --git a/docs/felupe/region.rst b/docs/felupe/region.rst index 5b6e7c51..0ca4a653 100644 --- a/docs/felupe/region.rst +++ b/docs/felupe/region.rst @@ -34,6 +34,7 @@ This module contains the definition of a region as well as a boundary region alo RegionQuadBoundary RegionHexahedronBoundary RegionVertex + RegionTruss **Detailed API Reference** @@ -135,4 +136,9 @@ This module contains the definition of a region as well as a boundary region alo .. autoclass:: felupe.RegionVertex :members: :undoc-members: - :show-inheritance: \ No newline at end of file + :show-inheritance: + +.. autoclass:: felupe.RegionTruss + :members: + :undoc-members: + :show-inheritance: diff --git a/src/felupe/__init__.py b/src/felupe/__init__.py index 9a609023..f0d65c11 100644 --- a/src/felupe/__init__.py +++ b/src/felupe/__init__.py @@ -122,6 +122,7 @@ RegionTriangleMINI, RegionTriQuadraticHexahedron, RegionTriQuadraticHexahedronBoundary, + RegionTruss, RegionVertex, ) from .tools import hello_world, newtonrhapson, project, runs_on, save, topoints @@ -250,6 +251,7 @@ "RegionTriangleMINI", "RegionTriQuadraticHexahedron", "RegionTriQuadraticHexahedronBoundary", + "RegionTruss", "RegionVertex", "newtonrhapson", "project", diff --git a/src/felupe/mechanics/_truss.py b/src/felupe/mechanics/_truss.py index 9c0ae507..1f2722ab 100644 --- a/src/felupe/mechanics/_truss.py +++ b/src/felupe/mechanics/_truss.py @@ -50,31 +50,34 @@ class TrussBody(Solid): Examples -------- .. pyvista-plot:: + :context: :force_static: >>> import felupe as fem >>> >>> mesh = fem.Mesh( - ... points=[[0, 0], [1, 1], [2.0, 0]], cells=[[0, 1], [1, 2]], cell_type="line" + ... points=[[0, 0], [1, 1], [2.0, 0]], + ... cells=[[0, 1], [1, 2]], + ... cell_type="line", ... ) - >>> region = fem.Region(mesh, fem.Line(), fem.GaussLobatto(order=0, dim=1), grad=False) + >>> region = fem.RegionTruss(mesh) >>> field = fem.Field(region, dim=2).as_container() >>> boundaries = fem.BoundaryDict(fixed=fem.Boundary(field[0], fy=0)) >>> >>> umat = fem.LinearElastic1D(E=[1, 1]) >>> truss = fem.TrussBody(umat, field, area=[1, 1]) - >>> load = fem.PointLoad(field, [1]) + >>> load = fem.PointLoad(field, points=[1]) >>> - >>> mesh.plot(plotter=load.plot(plotter=boundaries.plot()), line_width=5).show() - >>> - >>> move = fem.math.linsteps([0, -0.1], num=5, axis=1, axes=2) - >>> step = fem.Step(items=[truss, load], ramp={load: move}, boundaries=boundaries) + >>> table = fem.math.linsteps([0, 1], num=5, axis=1, axes=2) + >>> ramp = {load: table * -0.1} + >>> step = fem.Step(items=[truss, load], ramp=ramp, boundaries=boundaries) >>> job = fem.Job(steps=[step]).evaluate() + >>> + >>> mesh.plot( + ... plotter=load.plot(plotter=boundaries.plot(), scale=0.5), + ... line_width=5 + ... ).show() - See Also - -------- - felupe.SolidBodyNearlyIncompressible : A (nearly) incompressible solid body with - methods for the assembly of sparse vectors/matrices. """ def __init__(self, umat, field, area, statevars=None): diff --git a/src/felupe/region/__init__.py b/src/felupe/region/__init__.py index 4086467b..bf6f2a80 100644 --- a/src/felupe/region/__init__.py +++ b/src/felupe/region/__init__.py @@ -22,6 +22,7 @@ RegionTriangleMINI, RegionTriQuadraticHexahedron, RegionTriQuadraticHexahedronBoundary, + RegionTruss, RegionVertex, ) @@ -49,5 +50,6 @@ "RegionTriangleMINI", "RegionTriQuadraticHexahedron", "RegionTriQuadraticHexahedronBoundary", + "RegionTruss", "RegionVertex", ] diff --git a/src/felupe/region/_templates.py b/src/felupe/region/_templates.py index ec25f8cb..e47c699e 100644 --- a/src/felupe/region/_templates.py +++ b/src/felupe/region/_templates.py @@ -22,6 +22,7 @@ ConstantHexahedron, ConstantQuad, Hexahedron, + Line, Quad, QuadraticHexahedron, QuadraticQuad, @@ -35,7 +36,7 @@ Vertex, ) from ..mesh import Mesh -from ..quadrature import GaussLegendre, GaussLegendreBoundary +from ..quadrature import GaussLegendre, GaussLegendreBoundary, GaussLobatto from ..quadrature import Tetrahedron as TetraQuadrature from ..quadrature import Triangle as TriangleQuadrature from ._boundary import RegionBoundary @@ -863,3 +864,13 @@ def __init__( ): element = Vertex() super().__init__(mesh, element, quadrature, grad=grad, **kwargs) + + +class RegionTruss(Region): + "A region with a truss element." + + def __init__( + self, mesh, quadrature=GaussLobatto(order=0, dim=1), grad=False, **kwargs + ): + element = Line() + super().__init__(mesh, element, quadrature, grad=grad, **kwargs) diff --git a/tests/test_mechanics.py b/tests/test_mechanics.py index 1fd9fbb5..df8eee6f 100644 --- a/tests/test_mechanics.py +++ b/tests/test_mechanics.py @@ -609,7 +609,7 @@ def test_truss(): mesh = fem.Mesh( points=[[0, 0], [1, 1], [2.0, 0]], cells=[[0, 1], [1, 2]], cell_type="line" ) - region = fem.Region(mesh, fem.Line(), fem.GaussLobatto(order=0, dim=1), grad=False) + region = fem.RegionTruss(mesh) field = fem.Field(region, dim=2).as_container() boundaries = fem.BoundaryDict(fixed=fem.Boundary(field[0], fy=0)) From 0f41d8eafabddbbe3a6f1f03df465b7eea0889d3 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Sun, 1 Jun 2025 23:16:01 +0200 Subject: [PATCH 36/40] Enhance the docstring of `TrussBody` --- src/felupe/mechanics/_truss.py | 102 +++++++++++++++++++++++++++++++++ 1 file changed, 102 insertions(+) diff --git a/src/felupe/mechanics/_truss.py b/src/felupe/mechanics/_truss.py index 1f2722ab..e4f95627 100644 --- a/src/felupe/mechanics/_truss.py +++ b/src/felupe/mechanics/_truss.py @@ -47,6 +47,108 @@ class TrussBody(Solid): statevars : ndarray or None, optional Array of initial internal state variables (default is None). + Notes + ----- + For a truss element the stretch may be calculated as given in Eq. + :eq:`truss-stretch`. + + .. math:: + :label: truss-stretch + + \Lambda = \frac{l}{L} = \sqrt{1 + 2 \left( + \frac{\boldsymbol{\Delta X}}{L} + \right)^T \left(\frac{\boldsymbol{\Delta u}}{L}\right) + \left( + \frac{\boldsymbol{\Delta u}}{L} + \right)^T \left(\frac{\boldsymbol{\Delta u}}{L}\right)} + + This follows from Eqs. :eq:`truss-lengths` + + .. math:: + :label: truss-lengths + + l^2 &= \boldsymbol{\Delta x}^T \boldsymbol{\Delta x} \\ + L^2 &= \boldsymbol{\Delta X}^T \boldsymbol{\Delta X} + + and enables the Biot strain measure, see Eq. :eq:`truss-strain`. + + .. math:: + :label: truss-strain + + E_{11} = \Lambda - 1 + + The normal force of a truss depends directly on the geometric exactly defined + strain measure :math:`E_{11}`. For the general case of a nonlinear material + behaviour the normal force is defined as given in Eq. :eq:`truss-force` + + .. math:: + :label: truss-force + + N = S_{11}(E_{11})~A + + and the derivative according to Eq. :eq:`truss-force-derivative`. + + .. math:: + :label: truss-force-derivative + + \frac{\partial N}{\partial E_{11}} = \frac{ + \partial S_{11}(E_{11}) + }{\partial E_{11}}~A + + For the case of a linear elastic material this reduces to + Eqs. :eq:`truss-force-linear`. + + .. math:: + :label: truss-force-linear + + S_{11}(E_{11}) &= E~E_{11} \\ + N &= EA~E_{11} \\ + \frac{\partial N}{\partial E_{11}} &= EA + + The (nonlinear) fixing force column vector may be expressed as a function of the + elemental force :math:`N_k` and the deformed unit vector :math:`\boldsymbol{n}_k`, + see Eq. :eq:`truss-fixing-force`. + + .. math:: + :label: truss-fixing-force + + \boldsymbol{r}_k = \begin{bmatrix} + \boldsymbol{r}_A \\ + \boldsymbol{r}_E + \end{bmatrix} = N_k \begin{pmatrix} + -\boldsymbol{n}_k \\ + \phantom{-}\boldsymbol{n}_k + \end{pmatrix} + + For a truss the stiffness matrix is divided into four block matrices of the + same components but with different signs, see Eq. :eq:`truss-stiffness-matrix`. + + .. math:: + :label: truss-stiffness-matrix + + \boldsymbol{K}_{k~(6,6)} = \begin{bmatrix} + \boldsymbol{K}_{AA} & \boldsymbol{K}_{AE}\\ + \boldsymbol{K}_{EA} & \boldsymbol{K}_{EE} + \end{bmatrix} = \begin{pmatrix} + \phantom{-}\boldsymbol{K}_{EE} & -\boldsymbol{K}_{EE}\\ + -\boldsymbol{K}_{EE} & \phantom{-}\boldsymbol{K}_{EE} + \end{pmatrix} + + A change in the fixing force vector at the end node `E` w.r.t. a small change of + the displacements at node `E` is defined as the tangent stiffnes `EE`, see + Eq. :eq:`truss-stiffness-block`. + + .. math:: + :label: truss-stiffness-block + + \boldsymbol{K}_{EE} &= \frac{ + \partial \boldsymbol{r}_E}{\partial \boldsymbol{U}_E + } \\ + \boldsymbol{K}_{EE} &= \frac{A}{L} ~ \frac{ + \partial S_{11}(E_{11}) + }{\partial E_{11}} \boldsymbol{n} \otimes \boldsymbol{n} + \frac{N}{l} \left( + \boldsymbol{1} - \boldsymbol{n} \otimes \boldsymbol{n} + \right) + Examples -------- .. pyvista-plot:: From 08d9b4c889f90e0efe9a15ece3454abb342be97d Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Mon, 9 Jun 2025 13:55:13 +0200 Subject: [PATCH 37/40] Docs: Add howto/assemble --- docs/howto.rst | 1 + docs/howto/assemble.rst | 66 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 67 insertions(+) create mode 100644 docs/howto/assemble.rst diff --git a/docs/howto.rst b/docs/howto.rst index 0c4953e4..c16966c7 100644 --- a/docs/howto.rst +++ b/docs/howto.rst @@ -17,3 +17,4 @@ How-To Guides howto/umat_hyperelasticity howto/items howto/project + howto/assemble diff --git a/docs/howto/assemble.rst b/docs/howto/assemble.rst new file mode 100644 index 00000000..d4449650 --- /dev/null +++ b/docs/howto/assemble.rst @@ -0,0 +1,66 @@ +Assemble Vectors and Matrices +----------------------------- + +Integral (weak) forms are transformed into sparse vectors and matrices by calling the +assemble-method of an :class:`integral form`. The +:class:`Neo-Hookean ` material model formulation is used to evaluate +both the variation and linerization of its strain energy density function. + +.. plot:: + :context: + + import felupe as fem + + mesh = fem.Cube(n=3) + region = fem.RegionHexahedron(mesh) + umat = fem.NeoHooke(mu=1, bulk=2) + field = fem.Field(region, dim=3).as_container() + + F = field.extract(grad=True, sym=False, add_identity=True) # deformation-gradient + P = umat.gradient(F)[:1] # list with first Piola-Kirchhoff stress tensor + A = umat.hessian(F) # list with fourth-order elasticity tensor + +The code for the variation of the total-potential energy, as given in Eq. +:eq:`total-potential-energy-variation`, is very close to the analytical expression. + +.. math:: + :label: total-potential-energy-variation + + \delta \Pi = \int_V \boldsymbol{P} : \frac{ + \partial \delta \boldsymbol{u} + }{\partial \boldsymbol{X}} \ dV + +.. plot:: + :context: + + δΠ = fem.IntegralForm( # variation of total potential energy + fun=P, + v=field, # container for the test field + dV=region.dV, # differential volume + grad_v=[True], # use gradient of test field + ) + vector = δΠ.assemble() + +For the linearization, see Eq. eq:`total-potential-energy-linearization`. + +.. math:: + :label: total-potential-energy-linearization + + \Delta \delta \Pi = \int_V \frac{ + \partial \delta \boldsymbol{u} + }{\partial \boldsymbol{X}} : \mathbb{A} : \frac{ + \partial \Delta \boldsymbol{u} + }{\partial \boldsymbol{X}} \ dV + +.. plot:: + :context: + + ΔδΠ = fem.IntegralForm( # linearization of total potential energy + fun=A, + v=field, # container for the test field + u=field, # container for the trial field + dV=region.dV, # differential volume + grad_v=[True], # use gradient of test field + grad_u=[True], # use gradient of trial field + ) + matrix = ΔδΠ.assemble() From 488bd34f53a73cf6cc6c8d77d672f94defaff915 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Mon, 9 Jun 2025 22:48:58 +0200 Subject: [PATCH 38/40] Update assemble.rst fix missing : for Eq. 2 --- docs/howto/assemble.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/howto/assemble.rst b/docs/howto/assemble.rst index d4449650..686f25b6 100644 --- a/docs/howto/assemble.rst +++ b/docs/howto/assemble.rst @@ -41,7 +41,7 @@ The code for the variation of the total-potential energy, as given in Eq. ) vector = δΠ.assemble() -For the linearization, see Eq. eq:`total-potential-energy-linearization`. +For the linearization, see Eq. :eq:`total-potential-energy-linearization`. .. math:: :label: total-potential-energy-linearization From 6416c7d6967ef88092ac8c55566c571a195bcda1 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Thu, 12 Jun 2025 15:22:38 +0200 Subject: [PATCH 39/40] Enhance the docstring of `TrussBody` --- src/felupe/mechanics/_truss.py | 14 +++++--------- 1 file changed, 5 insertions(+), 9 deletions(-) diff --git a/src/felupe/mechanics/_truss.py b/src/felupe/mechanics/_truss.py index e4f95627..69ea8073 100644 --- a/src/felupe/mechanics/_truss.py +++ b/src/felupe/mechanics/_truss.py @@ -49,19 +49,15 @@ class TrussBody(Solid): Notes ----- - For a truss element the stretch may be calculated as given in Eq. - :eq:`truss-stretch`. + For a truss element the stretch is evaluated as given in Eq. :eq:`truss-stretch` .. math:: :label: truss-stretch - \Lambda = \frac{l}{L} = \sqrt{1 + 2 \left( - \frac{\boldsymbol{\Delta X}}{L} - \right)^T \left(\frac{\boldsymbol{\Delta u}}{L}\right) + \left( - \frac{\boldsymbol{\Delta u}}{L} - \right)^T \left(\frac{\boldsymbol{\Delta u}}{L}\right)} + \Lambda = \sqrt{\frac{l^2}{L^2}} - This follows from Eqs. :eq:`truss-lengths` + with the squared undeformed and deformed lengths, denoted in + Eqs. :eq:`truss-lengths`. .. math:: :label: truss-lengths @@ -69,7 +65,7 @@ class TrussBody(Solid): l^2 &= \boldsymbol{\Delta x}^T \boldsymbol{\Delta x} \\ L^2 &= \boldsymbol{\Delta X}^T \boldsymbol{\Delta X} - and enables the Biot strain measure, see Eq. :eq:`truss-strain`. + This enables the Biot strain measure, see Eq. :eq:`truss-strain`. .. math:: :label: truss-strain From 4927e7e5e5f5ade62dea210d836a392d363b9685 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Sun, 29 Jun 2025 21:45:43 +0200 Subject: [PATCH 40/40] set version tag to 9.3.0 --- CHANGELOG.md | 2 ++ src/felupe/__about__.py | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 45f91ff1..55c7404f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,6 +3,8 @@ All notable changes to this project will be documented in this file. The format ## [Unreleased] +## [9.3.0] - 2025-06-29 + ### Added - Add `Mesh.add_points(points)` to update the mesh with additional points. - Add `Mesh.clear_points_without_cells()` to clear the list of points without cells (useful for center-points of multi-point constraints). diff --git a/src/felupe/__about__.py b/src/felupe/__about__.py index d155c125..1f28d710 100644 --- a/src/felupe/__about__.py +++ b/src/felupe/__about__.py @@ -1 +1 @@ -__version__ = "9.3.0-dev" +__version__ = "9.3.0"