From 9407533e62f612c68f2786b5f6e0fded1a9db8e4 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Tue, 26 Nov 2024 23:34:56 +0100 Subject: [PATCH 01/47] Add `SolidBody.evaluate.mass(density=1.0)` (#906) * Add `SolidBody.evaluate.mass(density=1.0)` * Add tests and update changelog * format black --- CHANGELOG.md | 3 +++ src/felupe/__about__.py | 2 +- .../jax/models/hyperelastic/_storakers.py | 4 +--- .../models/hyperelastic/_storakers.py | 5 +---- src/felupe/mechanics/_helpers.py | 3 ++- src/felupe/mechanics/_solidbody.py | 17 ++++++++++++++++- .../mechanics/_solidbody_incompressible.py | 17 ++++++++++++++++- tests/test_mechanics.py | 5 +++++ 8 files changed, 45 insertions(+), 11 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 68715dc2..d2eba599 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,6 +3,9 @@ All notable changes to this project will be documented in this file. The format ## [Unreleased] +### Added +- Add `SolidBody.assemble.mass(density=1.0)` and `SolidBodyNearlyIncompressible.assemble.mass(density=1.0)` to assemble the mass matrix. + ## [9.1.0] - 2024-11-23 ### Added diff --git a/src/felupe/__about__.py b/src/felupe/__about__.py index ba280739..982a458b 100644 --- a/src/felupe/__about__.py +++ b/src/felupe/__about__.py @@ -1 +1 @@ -__version__ = "9.1.0" +__version__ = "9.2.0-dev" diff --git a/src/felupe/constitution/jax/models/hyperelastic/_storakers.py b/src/felupe/constitution/jax/models/hyperelastic/_storakers.py index 1d905aff..b5dac2ff 100644 --- a/src/felupe/constitution/jax/models/hyperelastic/_storakers.py +++ b/src/felupe/constitution/jax/models/hyperelastic/_storakers.py @@ -33,6 +33,4 @@ def storakers(C, mu, alpha, beta): α = array(alpha) β = array(beta) - return asum( - 2 * μ / α**2 * (λ1**α + λ2**α + λ3**α - 3 + (J ** (-α * β) - 1) / β) - ) + return asum(2 * μ / α**2 * (λ1**α + λ2**α + λ3**α - 3 + (J ** (-α * β) - 1) / β)) diff --git a/src/felupe/constitution/tensortrax/models/hyperelastic/_storakers.py b/src/felupe/constitution/tensortrax/models/hyperelastic/_storakers.py index 6fac588c..118606ba 100644 --- a/src/felupe/constitution/tensortrax/models/hyperelastic/_storakers.py +++ b/src/felupe/constitution/tensortrax/models/hyperelastic/_storakers.py @@ -116,10 +116,7 @@ def storakers(C, mu, alpha, beta): return tsum( [ - 2 - * μ - / α**2 - * (tsum(λ2 ** (α / 2)) - 3 + (det(C) ** (-α * β / 2) - 1) / β) + 2 * μ / α**2 * (tsum(λ2 ** (α / 2)) - 3 + (det(C) ** (-α * β / 2) - 1) / β) for μ, α, β in zip(mu, alpha, beta) ] ) diff --git a/src/felupe/mechanics/_helpers.py b/src/felupe/mechanics/_helpers.py index 7d1e8965..dc3b5146 100644 --- a/src/felupe/mechanics/_helpers.py +++ b/src/felupe/mechanics/_helpers.py @@ -26,9 +26,10 @@ class Assemble: "A class with methods for assembling vectors and matrices of an Item." - def __init__(self, vector, matrix, multiplier=None): + def __init__(self, vector, matrix, mass=None, multiplier=None): self.vector = vector self.matrix = matrix + self.mass = mass self.multiplier = multiplier diff --git a/src/felupe/mechanics/_solidbody.py b/src/felupe/mechanics/_solidbody.py index e2239f55..e7b9cd7b 100644 --- a/src/felupe/mechanics/_solidbody.py +++ b/src/felupe/mechanics/_solidbody.py @@ -230,7 +230,9 @@ def __init__(self, umat, field, statevars=None): ) ) - self.assemble = Assemble(vector=self._vector, matrix=self._matrix) + self.assemble = Assemble( + vector=self._vector, matrix=self._matrix, mass=self._mass + ) self.evaluate = Evaluate( gradient=self._gradient, @@ -361,3 +363,16 @@ def _cauchy_stress(self, field=None): J = det(F) return dot(P, transpose(F)) / J + + def _mass(self, density=1.0): + field = self.field[0].as_container() + dim = field[0].dim + form = self._form( + fun=[density * np.eye(dim).reshape(dim, dim, 1, 1)], + v=field, + u=field, + dV=field.region.dV, + grad_v=[False], + grad_u=[False], + ) + return form.assemble() diff --git a/src/felupe/mechanics/_solidbody_incompressible.py b/src/felupe/mechanics/_solidbody_incompressible.py index a8a1ed38..dfa5264b 100644 --- a/src/felupe/mechanics/_solidbody_incompressible.py +++ b/src/felupe/mechanics/_solidbody_incompressible.py @@ -355,7 +355,9 @@ def __init__(self, umat, field, bulk, state=None, statevars=None): self.results.state = state self.results.kinematics = self._extract(self.field) - self.assemble = Assemble(vector=self._vector, matrix=self._matrix) + self.assemble = Assemble( + vector=self._vector, matrix=self._matrix, mass=self._mass + ) self.evaluate = Evaluate( gradient=self._gradient, @@ -537,3 +539,16 @@ def _cauchy_stress(self, field=None): J = det(F) return dot(P, transpose(F)) / J + + def _mass(self, density=1.0): + field = self.field[0].as_container() + dim = field[0].dim + form = self._form( + fun=[density * np.eye(dim).reshape(dim, dim, 1, 1)], + v=field, + u=field, + dV=field.region.dV, + grad_v=[False], + grad_u=[False], + ) + return form.assemble() diff --git a/tests/test_mechanics.py b/tests/test_mechanics.py index 96bdbf30..82897e26 100644 --- a/tests/test_mechanics.py +++ b/tests/test_mechanics.py @@ -43,6 +43,7 @@ def test_simple(): r = b.assemble.vector() K = b.assemble.matrix() + M = b.assemble.mass() r = b.assemble.vector(u) F = b.results.kinematics P = b.results.stress @@ -51,6 +52,7 @@ def test_simple(): C = b.results.elasticity assert K.shape == (81, 81) + assert M.shape == (81, 81) assert r.shape == (81, 1) assert F[0].shape == (3, 3, 8, 8) assert P[0].shape == (3, 3, 8, 8) @@ -236,6 +238,9 @@ def test_solidbody_incompressible(): umat=umat, field=u, bulk=5000, state=fem.StateNearlyIncompressible(u) ) + M = b.assemble.mass() + assert M.shape == (81, 81) + for parallel in [False, True]: kwargs = {"parallel": parallel} From e2305db98c12a5f2dfe1c81c31ca9f75c8f5df3d Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Wed, 27 Nov 2024 17:39:21 +0100 Subject: [PATCH 02/47] Add `SolidBody.evaluate.stress()` and autodetect the stress-type in `SolidBody.plot(name)` from `name` (#908) * Autodetect stress-type in `SolidBody.plot()` from string * Update CHANGELOG.md * Update _solidbody.py * Update test_readme.py * Update test_mechanics.py * Update _solidbody.py --- CHANGELOG.md | 5 +++++ examples/ex02_plate-with-hole.py | 6 +++--- examples/ex11_notch-stress.py | 6 +++--- src/felupe/mechanics/_helpers.py | 9 +++++++- src/felupe/mechanics/_solidbody.py | 34 ++++++++++++++++++++++++++++-- src/felupe/view/_solid.py | 12 +++++++---- tests/test_mechanics.py | 12 ++++++----- tests/test_readme.py | 7 ------ 8 files changed, 66 insertions(+), 25 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index d2eba599..1afcc43f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,11 @@ All notable changes to this project will be documented in this file. The format ### Added - Add `SolidBody.assemble.mass(density=1.0)` and `SolidBodyNearlyIncompressible.assemble.mass(density=1.0)` to assemble the mass matrix. +- Add `SolidBody.evaluate.stress(field)` to evaluate the (first Piola-Kirchhoff) stress tensor (engineering stress in linear elasticity). + +### Changed +- The first Piola-Kirchhoff stress tensor is evaluated if `ViewSolid(stress_type=None)`. +- Autodetect the stress-type in `SolidBody.plot(name)` from `name`. ## [9.1.0] - 2024-11-23 diff --git a/examples/ex02_plate-with-hole.py b/examples/ex02_plate-with-hole.py index b7b769ba..e4b600c8 100644 --- a/examples/ex02_plate-with-hole.py +++ b/examples/ex02_plate-with-hole.py @@ -76,7 +76,7 @@ load = fem.SolidBodyPressure(field_boundary, pressure=-100) # %% -# The simulation model is now ready to be solved. The equivalent von Mises Cauchy stress +# The simulation model is now ready to be solved. The equivalent von Mises stress # will be plotted. For the two-dimensional case it is calculated by Eq. :eq:`svM`. # Stresses, located at quadrature-points of cells, are shifted to and averaged at mesh- # points. @@ -89,7 +89,7 @@ step = fem.Step(items=[solid, load], boundaries=boundaries) job = fem.Job(steps=[step]).evaluate() -solid.plot("Equivalent of Cauchy Stress", show_edges=False, project=fem.topoints).show() +solid.plot("Equivalent of Stress", show_edges=False, project=fem.topoints).show() # %% # The normal stress distribution :math:`\sigma_{11}` over the hole at :math:`x=0` is @@ -97,7 +97,7 @@ import matplotlib.pyplot as plt plt.plot( - fem.tools.project(solid.evaluate.cauchy_stress(), region)[:, 0, 0][mesh.x == 0], + fem.tools.project(solid.evaluate.stress(), region)[:, 0, 0][mesh.x == 0], mesh.points[:, 1][mesh.x == 0] / h, "o-", ) diff --git a/examples/ex11_notch-stress.py b/examples/ex11_notch-stress.py index 3dc4d7ea..31b1cfa9 100644 --- a/examples/ex11_notch-stress.py +++ b/examples/ex11_notch-stress.py @@ -62,7 +62,7 @@ job = fem.Job(steps=[step]).evaluate(parallel=True, solver=pypardiso.spsolve) solid.plot( - "Principal Values of Cauchy Stress", + "Principal Values of Stress", show_edges=False, view="xy", project=fem.topoints, @@ -72,7 +72,7 @@ # %% # The number of maximum endurable cycles between zero and the applied displacement is # evaluated with a SN-curve as denoted in Eq. :eq:`sn-curve`. The range of the maximum -# principal value of the Cauchy stress tensor is used to evaluate the fatigue life. +# principal value of the stress tensor is used to evaluate the fatigue life. # For simplicity, the stress is evaluated for the total solid body. To consider only # stresses on points which lie on the surface of the solid body, the cells on faces # :meth:`~felupe.RegionHexahedronBoundary.mesh.cells_faces` must be determined @@ -87,7 +87,7 @@ N_D = 2e6 # cycles k = 5 # slope -S = fem.topoints(fem.math.eigvalsh(solid.evaluate.cauchy_stress())[-1], region) +S = fem.topoints(fem.math.eigvalsh(solid.evaluate.stress())[-1], region) N = N_D * (abs(S) / S_D) ** -k view = solid.view(point_data={"Endurable Cycles": N}) diff --git a/src/felupe/mechanics/_helpers.py b/src/felupe/mechanics/_helpers.py index dc3b5146..f511ea17 100644 --- a/src/felupe/mechanics/_helpers.py +++ b/src/felupe/mechanics/_helpers.py @@ -36,7 +36,9 @@ def __init__(self, vector, matrix, mass=None, multiplier=None): class Evaluate: "A class with evaluate methods of an Item." - def __init__(self, gradient, hessian, cauchy_stress=None, kirchhoff_stress=None): + def __init__( + self, gradient, hessian, stress=None, cauchy_stress=None, kirchhoff_stress=None + ): self.gradient = gradient self.hessian = hessian @@ -44,6 +46,11 @@ def __init__(self, gradient, hessian, cauchy_stress=None, kirchhoff_stress=None) self.cauchy_stress = cauchy_stress self.kirchhoff_stress = kirchhoff_stress + if stress is None: + stress = lambda field=None, **kwargs: gradient(field, **kwargs)[0] + + self.stress = stress + class Results: "A class with intermediate results of an Item." diff --git a/src/felupe/mechanics/_solidbody.py b/src/felupe/mechanics/_solidbody.py index e7b9cd7b..b8c64175 100644 --- a/src/felupe/mechanics/_solidbody.py +++ b/src/felupe/mechanics/_solidbody.py @@ -31,7 +31,14 @@ class Solid: "Base class for solid bodies which provides methods for visualisations." - def view(self, point_data=None, cell_data=None, cell_type=None, project=None): + def view( + self, + point_data=None, + cell_data=None, + cell_type=None, + project=None, + stress_type="Cauchy", + ): """View the solid with optional given dicts of point- and cell-data items. Parameters @@ -44,6 +51,10 @@ def view(self, point_data=None, cell_data=None, cell_type=None, project=None): Cell-type of PyVista (default is None). project : callable or None, optional Project stress at quadrature-points to mesh-points (default is None). + stress_type : str or None, optional + The type of stress which is exported, either "Cauchy", "Kirchhoff" or None. + If None, the first Piola-Kirchhoff stress (engineering stress in linear + elasticity) is used. Default is "Cauchy". Returns ------- @@ -65,6 +76,7 @@ def view(self, point_data=None, cell_data=None, cell_type=None, project=None): cell_data=cell_data, cell_type=cell_type, project=project, + stress_type=stress_type, ) def plot(self, *args, project=None, **kwargs): @@ -76,7 +88,25 @@ def plot(self, *args, project=None, **kwargs): felupe.project: Project given values at quadrature-points to mesh-points. felupe.topoints: Shift given values at quadrature-points to mesh-points. """ - return self.view(project=project).plot(*args, **kwargs) + + stress_type = "Cauchy" + + if len(args) > 0: + name = kwargs.pop("name", args[0]) + + if "Stress" in name: + stress_type = ( + name.lower() + .split("principal values of ")[-1] + .split("equivalent of ")[-1] + .split("stress")[0] + .rstrip() + ) + + if len(stress_type) == 0: + stress_type = None + + return self.view(project=project, stress_type=stress_type).plot(*args, **kwargs) def screenshot( self, diff --git a/src/felupe/view/_solid.py b/src/felupe/view/_solid.py index 7ecfbbd0..2106e052 100644 --- a/src/felupe/view/_solid.py +++ b/src/felupe/view/_solid.py @@ -34,9 +34,10 @@ class ViewSolid(ViewField): The field-container. solid : felupe.SolidBody or felupe.SolidBodyIncompressible or None, optional A solid body to evaluate the (Cauchy) stress (default is None). - stress_type : str, optional - The type of stress which is exported, either "Cauchy" or "Kirchhoff" (default is - "Cauchy"). + stress_type : str or None, optional + The type of stress which is exported, either "Cauchy", "Kirchhoff" or None. If + None, the first Piola-Kirchhoff stress (engineering stress in linear elasticity) + is used. Default is "Cauchy". point_data : dict or None, optional Additional point-data dict (default is None). cell_data : dict or None, optional @@ -101,12 +102,15 @@ def __init__( cell_data_from_solid = {} if solid is not None: + if stress_type is None: + stress_type = "" stress_from_field = { + "": solid.evaluate.stress, "cauchy": solid.evaluate.cauchy_stress, "kirchhoff": solid.evaluate.kirchhoff_stress, } stress = stress_from_field[stress_type.lower()](field) - stress_label = f"{stress_type.title()} Stress" + stress_label = f"{stress_type.title()} Stress".lstrip() if project is None: cell_data_from_solid[stress_label] = tovoigt(stress.mean(-2)).T diff --git a/tests/test_mechanics.py b/tests/test_mechanics.py index 82897e26..4039afc5 100644 --- a/tests/test_mechanics.py +++ b/tests/test_mechanics.py @@ -276,6 +276,10 @@ def test_solidbody_incompressible(): F2 = b._extract(u) assert np.allclose(F1, F2) + p1 = b.evaluate.stress() + p2 = b.evaluate.stress(u) + assert np.allclose(p1, p2) + s1 = b.evaluate.cauchy_stress() s2 = b.evaluate.cauchy_stress(u) assert np.allclose(s1, s2) @@ -512,7 +516,7 @@ def test_view(): field = fem.FieldContainer([fem.FieldPlaneStrain(region, dim=2)]) umat = fem.NeoHooke(mu=1, bulk=2) solid = fem.SolidBody(umat, field) - plotter = solid.plot(off_screen=True) + plotter = solid.plot("Principal Values of Cauchy Stress", off_screen=True) # img = solid.screenshot(transparent_background=True) # ax = solid.imshow() @@ -520,12 +524,10 @@ def test_view(): umat = fem.LinearElasticPlaneStress(E=1, nu=0.3) solid = fem.SolidBody(umat, field) - with pytest.warns(): - plotter = solid.plot(off_screen=True) + plotter = solid.plot("Equivalent of Stress", off_screen=True) solid = fem.SolidBodyNearlyIncompressible(umat, field, bulk=0) - with pytest.warns(): - plotter = solid.plot(off_screen=True) + plotter = solid.plot("Kirchhoff Stress", off_screen=True) def test_threefield(): diff --git a/tests/test_readme.py b/tests/test_readme.py index 77726bac..ebf9f68d 100644 --- a/tests/test_readme.py +++ b/tests/test_readme.py @@ -1,10 +1,3 @@ -# -*- coding: utf-8 -*- -""" -Created on Wed Oct 27 15:33:53 2021 - -@author: adutz -""" - import numpy as np import felupe as fem From c77b69f6bf20e3d1529f74ff978da0978c102c1e Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Thu, 28 Nov 2024 21:17:13 +0100 Subject: [PATCH 03/47] Add `FreeVibration(items, boundaries)` (#909) * Add `FreeVibration(items, boundaries)` with methods to `evaluate()` and `extract()` the n-th eigenvector to a field * Update _free_vibration.py * Create test_free_vibration.py * add test * Update test_mechanics.py * Update test_mechanics.py * add test, format black * Update test_free_vibration.py * update the free-vibration test * Update mechanics.rst --- docs/felupe/mechanics.rst | 6 + src/felupe/__init__.py | 2 + src/felupe/mechanics/__init__.py | 2 + src/felupe/mechanics/_free_vibration.py | 132 ++++++++++++++++++ src/felupe/mechanics/_solidbody.py | 13 +- .../mechanics/_solidbody_incompressible.py | 12 +- tests/test_free_vibration.py | 81 +++++++++++ tests/test_mechanics.py | 8 +- 8 files changed, 250 insertions(+), 6 deletions(-) create mode 100644 src/felupe/mechanics/_free_vibration.py create mode 100644 tests/test_free_vibration.py diff --git a/docs/felupe/mechanics.rst b/docs/felupe/mechanics.rst index 7afe0792..2f9d4d8e 100644 --- a/docs/felupe/mechanics.rst +++ b/docs/felupe/mechanics.rst @@ -22,6 +22,7 @@ Mechanics Step Job CharacteristicCurve + FreeVibration **Point Load and Multi-Point Constraints** @@ -87,6 +88,11 @@ Mechanics :undoc-members: :show-inheritance: +.. autoclass:: felupe.FreeVibration + :members: + :undoc-members: + :show-inheritance: + .. autoclass:: felupe.MultiPointConstraint :members: :undoc-members: diff --git a/src/felupe/__init__.py b/src/felupe/__init__.py index 91427396..dd1c2fcd 100644 --- a/src/felupe/__init__.py +++ b/src/felupe/__init__.py @@ -70,6 +70,7 @@ from .mechanics import ( CharacteristicCurve, FormItem, + FreeVibration, Job, MultiPointConstraint, MultiPointContact, @@ -133,6 +134,7 @@ "tools", "Form", "FormItem", + "FreeVibration", "IntegralForm", "Basis", "Field", diff --git a/src/felupe/mechanics/__init__.py b/src/felupe/mechanics/__init__.py index 9b14d26b..69266ada 100644 --- a/src/felupe/mechanics/__init__.py +++ b/src/felupe/mechanics/__init__.py @@ -1,4 +1,5 @@ from ._curve import CharacteristicCurve +from ._free_vibration import FreeVibration from ._helpers import Assemble, Evaluate, Results, StateNearlyIncompressible from ._item import FormItem from ._job import Job @@ -16,6 +17,7 @@ "Assemble", "CharacteristicCurve", "Evaluate", + "FreeVibration", "FormItem", "StateNearlyIncompressible", "Job", diff --git a/src/felupe/mechanics/_free_vibration.py b/src/felupe/mechanics/_free_vibration.py new file mode 100644 index 00000000..fa149b88 --- /dev/null +++ b/src/felupe/mechanics/_free_vibration.py @@ -0,0 +1,132 @@ +# -*- 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 scipy.sparse import csr_matrix +from scipy.sparse.linalg import eigsh + +from ..dof import partition + + +class FreeVibration: + """A Free-Vibration Step/Job. + + Parameters + ---------- + items : list of SolidBody or SolidBodyNearlyIncompressible + A list of items with methods for the assembly of sparse stiffness and mass + matrices. + boundaries : dict of Boundary, optional + A dict with :class:`~felupe.Boundary` conditions (default is None). + + Notes + ----- + .. note:: + + Boundary conditions with non-zero values are not supported. + + Examples + -------- + .. pyvista-plot:: + + >>> import felupe as fem + >>> import numpy as np + >>> + >>> mesh = fem.Rectangle(b=(5, 1), n=(50, 10)) + >>> region = fem.RegionQuad(mesh) + >>> field = fem.FieldContainer([fem.FieldPlaneStrain(region, dim=2)]) + >>> + >>> boundaries = dict(left=fem.Boundary(field[0], fx=0)) + >>> solid = fem.SolidBody(fem.LinearElastic(E=2.5, nu=0.25), field, density=1.0) + >>> modal = fem.FreeVibration(items=[solid], boundaries=boundaries).evaluate() + >>> + >>> eigenvector, frequency = modal.extract(n=4, inplace=True) + >>> solid.plot("Stress", component=0).show() + """ + + def __init__(self, items, boundaries=None): + self.items = items + + if boundaries is None: + boundaries = {} + + self.boundaries = boundaries + self.eigenvalues = None + self.eigenvectors = None + + def evaluate(self, x0=None, solver=eigsh, parallel=False, **kwargs): + if x0 is not None: + x = x0 + else: + x = self.items[0].field + + # link field of items with global field + [item.field.link(x) for item in self.items] + + # init matrix with shape from global field + shape = (np.sum(x.fieldsizes), np.sum(x.fieldsizes)) + stiffness = csr_matrix(shape) + mass = csr_matrix(shape) + + # assemble matrices + for item in self.items: + K = item.assemble.matrix(parallel=parallel) + M = item.assemble.mass() + + if item.assemble.multiplier is not None: + K *= item.assemble.multiplier + + # check and reshape matrices + if K.shape != stiffness.shape: + K.resize(*shape) + + if M.shape != mass.shape: + M.resize(*shape) + + # add matrices + stiffness += K + mass += M + + dof0, self.dof1 = partition(x, self.boundaries) + + K = stiffness[self.dof1][:, self.dof1] + M = mass[self.dof1][:, self.dof1] + + sigma = kwargs.get("sigma", 0) + self.eigenvalues, self.eigenvectors = solver(A=K, M=M, sigma=sigma, **kwargs) + + return self + + def extract(self, n=0, x0=None, inplace=True): + if x0 is not None: + field = x0 + else: + field = self.items[0].field + + if not inplace: + field = field.copy() + + values = np.zeros(sum(field.fieldsizes)) + values[self.dof1] = self.eigenvectors[:, n] + + frequency = np.sqrt(self.eigenvalues[n]) / (2 * np.pi) + + [f.fill(0) for f in field] + field += values + + return field, frequency diff --git a/src/felupe/mechanics/_solidbody.py b/src/felupe/mechanics/_solidbody.py index b8c64175..fa806639 100644 --- a/src/felupe/mechanics/_solidbody.py +++ b/src/felupe/mechanics/_solidbody.py @@ -162,6 +162,8 @@ class SolidBody(Solid): A field container with one or more fields. statevars : ndarray or None, optional Array of initial internal state variables (default is None). + density : float or None, optional + The density of the solid body. Notes ----- @@ -239,9 +241,10 @@ class SolidBody(Solid): methods for the assembly of sparse vectors/matrices. """ - def __init__(self, umat, field, statevars=None): + def __init__(self, umat, field, statevars=None, density=None): self.umat = umat self.field = field + self.density = density self.results = Results(stress=True, elasticity=True) self.results.kinematics = self._extract(self.field) @@ -394,9 +397,14 @@ def _cauchy_stress(self, field=None): return dot(P, transpose(F)) / J - def _mass(self, density=1.0): + def _mass(self, density=None): + + if density is None: + density = self.density + field = self.field[0].as_container() dim = field[0].dim + form = self._form( fun=[density * np.eye(dim).reshape(dim, dim, 1, 1)], v=field, @@ -405,4 +413,5 @@ def _mass(self, density=1.0): grad_v=[False], grad_u=[False], ) + return form.assemble() diff --git a/src/felupe/mechanics/_solidbody_incompressible.py b/src/felupe/mechanics/_solidbody_incompressible.py index dfa5264b..e52aa69c 100644 --- a/src/felupe/mechanics/_solidbody_incompressible.py +++ b/src/felupe/mechanics/_solidbody_incompressible.py @@ -56,6 +56,8 @@ class SolidBodyNearlyIncompressible(Solid): A valid initial state for a (nearly) incompressible solid (default is None). statevars : ndarray or None, optional Array of initial internal state variables (default is None). + density : float or None, optional + The density of the solid body. Notes ----- @@ -315,10 +317,11 @@ class SolidBodyNearlyIncompressible(Solid): for (nearly) incompressible solid bodies. """ - def __init__(self, umat, field, bulk, state=None, statevars=None): + def __init__(self, umat, field, bulk, state=None, statevars=None, density=None): self.umat = umat self.field = field self.bulk = bulk + self.density = density self._area_change = AreaChange() self._form = IntegralForm @@ -540,9 +543,13 @@ def _cauchy_stress(self, field=None): return dot(P, transpose(F)) / J - def _mass(self, density=1.0): + def _mass(self, density=None): + if density is None: + density = self.density + field = self.field[0].as_container() dim = field[0].dim + form = self._form( fun=[density * np.eye(dim).reshape(dim, dim, 1, 1)], v=field, @@ -551,4 +558,5 @@ def _mass(self, density=1.0): grad_v=[False], grad_u=[False], ) + return form.assemble() diff --git a/tests/test_free_vibration.py b/tests/test_free_vibration.py new file mode 100644 index 00000000..83708194 --- /dev/null +++ b/tests/test_free_vibration.py @@ -0,0 +1,81 @@ +# -*- 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 + +import felupe as fem + + +def test_free_vibration(): + + mesh = fem.Cube(a=(0, 0, 5), b=(50, 100, 30), n=(3, 6, 4)) + region = fem.RegionHexahedron(mesh) + field = fem.FieldContainer([fem.Field(region, dim=3)]) + + umat = fem.NeoHooke(mu=1, bulk=2) + solid = fem.SolidBody(umat=umat, field=field, density=1.5e-9) + + job = fem.FreeVibration([solid]).evaluate() + new_field, frequency = job.extract(n=-1, inplace=False) + + +def test_free_vibration_mixed(): + + meshes = [ + fem.Cube(a=(0, 0, 30), b=(50, 100, 35), n=(3, 6, 2)), + fem.Cube(a=(0, 0, 5), b=(50, 100, 30), n=(3, 6, 4)), + fem.Cube(a=(0, 0, 0), b=(50, 100, 5), n=(3, 6, 2)), + ] + container = fem.MeshContainer(meshes, merge=True) + mesh = container.stack() + + regions = [fem.RegionHexahedron(m) for m in container.meshes] + fields = [ + fem.FieldsMixed(regions[0], n=1), + fem.FieldsMixed(regions[1], n=3), + fem.FieldsMixed(regions[2], n=1), + ] + + region = fem.RegionHexahedron(mesh) + field = fem.FieldContainer([fem.Field(region, dim=3), *fields[1][1:]]) + + boundaries = dict(left=fem.Boundary(field[0], fx=0)) + rubber = fem.ThreeFieldVariation(fem.NeoHooke(mu=1, bulk=5000)) + steel = fem.LinearElasticLargeStrain(2.1e5, 0.3) + solids = [ + fem.SolidBody(umat=steel, field=fields[0], density=7.85e-9), + fem.SolidBody(umat=rubber, field=fields[1], density=1.5e-9), + fem.SolidBody(umat=steel, field=fields[2], density=7.85e-9), + ] + + job = fem.FreeVibration(solids, boundaries).evaluate(x0=field) + new_field, frequency = job.extract(x0=field, n=-1, inplace=False) + + +if __name__ == "__main__": + test_free_vibration() + test_free_vibration_mixed() diff --git a/tests/test_mechanics.py b/tests/test_mechanics.py index 4039afc5..fdf3daa8 100644 --- a/tests/test_mechanics.py +++ b/tests/test_mechanics.py @@ -39,7 +39,7 @@ def test_simple(): r = fem.RegionHexahedron(m, uniform=True) u = fem.FieldContainer([fem.Field(r, dim=3)]) - b = fem.SolidBody(umat, u) + b = fem.SolidBody(umat, u, density=1.0) r = b.assemble.vector() K = b.assemble.matrix() @@ -235,7 +235,11 @@ def test_solidbody_incompressible(): umat = fem.OgdenRoxburgh(fem.NeoHooke(mu=1), r=3, m=1, beta=0) b = fem.SolidBodyNearlyIncompressible( - umat=umat, field=u, bulk=5000, state=fem.StateNearlyIncompressible(u) + umat=umat, + field=u, + bulk=5000, + state=fem.StateNearlyIncompressible(u), + density=1.0, ) M = b.assemble.mass() From 444a12f34ca29bb83c079a0b6a20e033f22745ba Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Thu, 28 Nov 2024 23:27:53 +0100 Subject: [PATCH 04/47] Enhance `hello_world()` with new arguments --- CHANGELOG.md | 2 + src/felupe/tools/_hello_world.py | 156 +++++++++++++++++++++++++++---- 2 files changed, 138 insertions(+), 20 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 1afcc43f..c6753cfc 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,10 +6,12 @@ All notable changes to this project will be documented in this file. The format ### Added - Add `SolidBody.assemble.mass(density=1.0)` and `SolidBodyNearlyIncompressible.assemble.mass(density=1.0)` to assemble the mass matrix. - Add `SolidBody.evaluate.stress(field)` to evaluate the (first Piola-Kirchhoff) stress tensor (engineering stress in linear elasticity). +- Add a free-vibration modal analysis Step/Job `FreeVibration(items, boundaries)` with methods to evaluate `FreeVibration.evaluate()` and to extract `field, frequency = FreeVibration.extract(n)` its n-th result. ### Changed - The first Piola-Kirchhoff stress tensor is evaluated if `ViewSolid(stress_type=None)`. - Autodetect the stress-type in `SolidBody.plot(name)` from `name`. +- Enhance the `hello_world(axisymmetric=False, planestrain=False, curve=False, xdmf=False, container=False)` function with new arguments to customize the generated template script. ## [9.1.0] - 2024-11-23 diff --git a/src/felupe/tools/_hello_world.py b/src/felupe/tools/_hello_world.py index c87fcc8c..120a8f0a 100644 --- a/src/felupe/tools/_hello_world.py +++ b/src/felupe/tools/_hello_world.py @@ -17,8 +17,16 @@ """ -def hello_world(pypardiso=False, parallel=False): - """Print felupe's hyperelastic hello world. +def hello_world( + pypardiso=False, + parallel=False, + axisymmetric=False, + planestrain=False, + curve=False, + xdmf=False, + container=False, +): + """Print FElupe's hyperelastic hello world. Parameters ---------- @@ -27,6 +35,16 @@ def hello_world(pypardiso=False, parallel=False): False). parallel : bool, optional Flag to activate a threaded vector- and matrix-assembly (default is False). + axisymmetric : bool, optional + Flag to create a template for an axisymmetric analysis (default is False). + planestrain : bool, optional + Flag to create a template for an plane strain analysis (default is False). + curve : bool, optional + Flag to use a characteristic-curve job (default is False). + xdmf : bool, optional + Flag to write a XDMF time-series result file (default is False). + container : bool, optional + Flag to use a mesh-container with multiple solid bodies (default is False). """ imports = [ @@ -42,23 +60,121 @@ def hello_world(pypardiso=False, parallel=False): if parallel: kwargs.append("parallel=True") - lines = [ - "mesh = fem.Cube(n=6)", - "region = fem.RegionHexahedron(mesh)", - "field = fem.FieldContainer([fem.Field(region, dim=3)])", - "", - "boundaries, loadcase = fem.dof.uniaxial(field, clamped=True)", - "umat = fem.NeoHooke(mu=1, bulk=2)", - "solid = fem.SolidBody(umat=umat, field=field)", - "", - "move = fem.math.linsteps([0, 1], num=5)", - 'ramp = {boundaries["move"]: move}', - "step = fem.Step(items=[solid], ramp=ramp, boundaries=boundaries)", - "", - "job = fem.Job(steps=[step])", - f'job.evaluate({", ".join(kwargs)})', - "", - 'ax = solid.imshow("Principal Values of Cauchy Stress")', - ] + if xdmf: + kwargs.append('filename="result.xdmf"') + + region = "Region" + field = "Field" + + if axisymmetric: + mesh = "Rectangle" + region += "Quad" + field += "Axisymmetric" + dim = 2 + + elif planestrain: + mesh = "Rectangle" + field += "PlaneStrain" + dim = 2 + + else: + mesh = "Cube" + region += "Hexahedron" + dim = 3 + + job = "Job" + kwargs_job = [] + plot = [] + if curve: + job = "CharacteristicCurve" + kwargs_job.append("") + kwargs_job.append('boundary=boundaries["move"]') + plot = [ + "fig, ax = job.plot(", + r' xlabel=r"Displacement $d$ in mm $\longrightarrow$",', + r' ylabel=r"Normal Force $F$ in N $\longrightarrow$",', + ")", + "", + ] + + kwargs_job = ", ".join(kwargs_job) + plot = "\n".join(plot) + + if not container: + + kwargs = ", ".join(kwargs) + + # fmt: off + lines = [ + f"mesh = fem.{mesh}(n=6)", + f"region = fem.{region}(mesh)", + f"field = fem.FieldContainer([fem.{field}(region, dim={dim})])", + "", + "boundaries, loadcase = fem.dof.uniaxial(field, clamped=True)", + "umat = fem.NeoHooke(mu=1, bulk=2)", + "solid = fem.SolidBody(umat=umat, field=field)", + "", + "move = fem.math.linsteps([0, 1], num=5)", + 'ramp = {boundaries["move"]: move}', + "step = fem.Step(items=[solid], ramp=ramp, boundaries=boundaries)", + "", + f"job = fem.{job}(steps=[step]{kwargs_job})", + f"job.evaluate({kwargs})", + f"{plot}", + 'ax = solid.imshow("Principal Values of Cauchy Stress")', + ] + # fmt: on + + else: + + kwargs.insert(0, "x0=field") + kwargs = ", ".join(kwargs) + + # fmt: off + lines = [ + "meshes = [", + f" fem.{mesh}(n=6),", + f" fem.{mesh}(n=6).translate(1, axis=0),", + "]", + "container = fem.MeshContainer(meshes, merge=True)", + "mesh = container.stack()", + "", + f"region = fem.{region}(mesh)", + f"field = fem.FieldContainer([fem.{field}(region, dim={dim})])", + "", + "regions = [", + f" fem.{region}(container.meshes[0]),", + f" fem.{region}(container.meshes[1]),", + "]", + "fields = [", + f" fem.FieldContainer([fem.{field}(regions[0], dim={dim})]),", + f" fem.FieldContainer([fem.{field}(regions[1], dim={dim})]),", + "]", + "", + "boundaries, loadcase = fem.dof.uniaxial(field, clamped=True)", + "umats = [", + " fem.LinearElasticLargeStrain(E=2.1e5, nu=0.3),", + " fem.NeoHooke(mu=1),", + "]", + "solids = [", + " fem.SolidBody(umat=umats[0], field=fields[0]),", + " fem.SolidBodyNearlyIncompressible(umat=umats[1], field=fields[1], bulk=5000),", + "]", + "", + "move = fem.math.linsteps([0, 1], num=5)", + 'ramp = {boundaries["move"]: move}', + "step = fem.Step(items=solids, ramp=ramp, boundaries=boundaries)", + "", + f"job = fem.{job}(steps=[step]{kwargs_job})", + f"job.evaluate({kwargs})", + f"{plot}", + 'plotter = solids[0].plot(show_undeformed=False, style="wireframe")', + 'solids[1].plot(', + ' "Principal Values of Cauchy Stress",', + ' plotter=plotter,', + ' show_undeformed=False', + ').show()', + ] + # fmt: off print("\n\n".join(["\n".join(imports), "\n".join(lines)])) From 4037031fc62c914cd263bc18ff4762d49e446a13 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Fri, 29 Nov 2024 00:25:05 +0100 Subject: [PATCH 05/47] Update test_tools.py --- tests/test_tools.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/tests/test_tools.py b/tests/test_tools.py index ff49a51f..81aafedc 100644 --- a/tests/test_tools.py +++ b/tests/test_tools.py @@ -105,7 +105,9 @@ def test_solve(): def test_hello_world(): fem.hello_world() - fem.hello_world(pypardiso=True, parallel=True) + fem.hello_world(pypardiso=True, parallel=True, xdmf=True) + fem.hello_world(curve=True, axisymmetric=True) + fem.hello_world(container=True, planestrain=True) def test_solve_mixed(): From 43a7705f232eb3e63b5e1da3568d8be91732b779 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Sat, 30 Nov 2024 21:35:34 +0100 Subject: [PATCH 06/47] `paper.md` (draft): simplify the statement on performance --- paper/paper.md | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/paper/paper.md b/paper/paper.md index 4dca46e7..d2859aa6 100644 --- a/paper/paper.md +++ b/paper/paper.md @@ -35,9 +35,7 @@ numerical solution of nonlinear problems in continuum mechanics of solid bodies. # Statement of need There are well-established Python packages available for finite element analysis. These packages are either distributed as binary packages or need to be compiled on installation, like FEniCSx [@fenicsx], GetFEM [@getfem] or SfePy [@sfepy]. JAX-FEM [@jaxfem], which is built on JAX [@jax], is a pure Python package but requires many dependencies in its recommended environment. `scikit-fem` [@scikitfem] is a pure Python package with minimal dependencies but with a more general scope [@scikitfem]. FElupe is both easy-to-install as well as easy-to-use in its target domain of hyperelastic solid bodies. -The performance of FElupe is good for a non-compiled package but mediocre in comparison to compiled codes. However, it is still well-suited for up to mid-sized problems, i.e. up to $10^5$ degrees of freedom, when basic hyperelastic model formulations are used. A performance benchmark for times spent on stiffness matrix assembly is included in the documentation. Internally, efficient NumPy [@numpy] based math is realized by element-wise operating trailing axes [@scikitfem]. An all-at-once approach per operation is used instead of a cell-by-cell evaluation loop. The constitutive material formulation class is backend agnostic: FElupe provides NumPy-arrays as input arguments and requires NumPy-arrays as return values. This enables backends like JAX [@jax] or PyTorch [@pytorch] to be used. Interactive views of meshes, fields and solid bodies are enabled by PyVista [@pyvista]. The capabilities of FElupe may be enhanced with additional Python packages, e.g. `meshio` [@meshio], `matadi` [@matadi], `tensortrax` [@tensortrax], `hyperelastic` [@hyperelastic] or `feplot` [@feplot]. - - +The performance of FElupe is good for a non-compiled package and it is well-suited for up to mid-sized problems, i.e. up to $10^5$ degrees of freedom, when hyperelastic model formulations are used. A performance benchmark for times spent on stiffness matrix assembly is included in the documentation. Internally, efficient NumPy [@numpy] based math is realized by element-wise operating trailing axes [@scikitfem]. An all-at-once approach per operation is used instead of a cell-by-cell evaluation loop. The constitutive material formulation class is backend agnostic: FElupe provides NumPy-arrays as input arguments and requires NumPy-arrays as return values. This enables backends like JAX [@jax] or PyTorch [@pytorch] to be used. Interactive views of meshes, fields and solid bodies are enabled by PyVista [@pyvista]. The capabilities of FElupe may be enhanced with additional Python packages, e.g. `meshio` [@meshio], `matadi` [@matadi], `tensortrax` [@tensortrax], `hyperelastic` [@hyperelastic] or `feplot` [@feplot]. # Features The essential high-level parts of solving problems with FElupe include a field, a solid body, boundary conditions and a job. With the combination of a mesh, a finite element formulation and a quadrature rule, a numeric region is created. A field for a field container is further created on top of this numeric region, see \autoref{fig:field}. From 088e64f5ffe678628940de0d9dfa478e2a9369f1 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Sat, 30 Nov 2024 22:32:01 +0100 Subject: [PATCH 07/47] Enhance `Boundary`: Add support for multiaxial prescribed `value`-argument (#912) * Enhance `Boundary`: Support multiaxial prescribed values and broadcast them (if necessary) * format black * Update CHANGELOG.md * Update test_dof.py --- CHANGELOG.md | 1 + src/felupe/dof/_tools.py | 15 ++++++++++++++- src/felupe/mechanics/_free_vibration.py | 6 +++--- src/felupe/mechanics/_solidbody.py | 1 - src/felupe/tools/_hello_world.py | 2 -- tests/test_dof.py | 18 ++++++++++++++++++ tests/test_free_vibration.py | 2 -- 7 files changed, 36 insertions(+), 9 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index c6753cfc..2bfea234 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,6 +12,7 @@ All notable changes to this project will be documented in this file. The format - The first Piola-Kirchhoff stress tensor is evaluated if `ViewSolid(stress_type=None)`. - Autodetect the stress-type in `SolidBody.plot(name)` from `name`. - Enhance the `hello_world(axisymmetric=False, planestrain=False, curve=False, xdmf=False, container=False)` function with new arguments to customize the generated template script. +- Enhance `Boundary` with added support for multiaxial prescribed values. ## [9.1.0] - 2024-11-23 diff --git a/src/felupe/dof/_tools.py b/src/felupe/dof/_tools.py index e21af72d..780764f6 100644 --- a/src/felupe/dof/_tools.py +++ b/src/felupe/dof/_tools.py @@ -235,8 +235,21 @@ def apply(field, bounds, dof0=None): # get offset for field-dof of current boundary offset = offsets[[b.field == f for f in field.fields]] + if isinstance(b.value, np.ndarray): + # get size and broadcast the values if necessary + value = b.value + if b.dof.size != value.size: + if len(value.shape) == 1: + value = value.reshape(1, -1) + + value = np.broadcast_to(value, (b.points.size, b.value.shape[-1])) + + value = value.ravel() + else: + value = b.value + # set prescribed values - u.ravel()[b.dof + offset] = b.value + u.ravel()[b.dof + offset] = value if dof0 is None: return u diff --git a/src/felupe/mechanics/_free_vibration.py b/src/felupe/mechanics/_free_vibration.py index fa149b88..4c27ae79 100644 --- a/src/felupe/mechanics/_free_vibration.py +++ b/src/felupe/mechanics/_free_vibration.py @@ -33,11 +33,11 @@ class FreeVibration: matrices. boundaries : dict of Boundary, optional A dict with :class:`~felupe.Boundary` conditions (default is None). - + Notes ----- .. note:: - + Boundary conditions with non-zero values are not supported. Examples @@ -54,7 +54,7 @@ class FreeVibration: >>> boundaries = dict(left=fem.Boundary(field[0], fx=0)) >>> solid = fem.SolidBody(fem.LinearElastic(E=2.5, nu=0.25), field, density=1.0) >>> modal = fem.FreeVibration(items=[solid], boundaries=boundaries).evaluate() - >>> + >>> >>> eigenvector, frequency = modal.extract(n=4, inplace=True) >>> solid.plot("Stress", component=0).show() """ diff --git a/src/felupe/mechanics/_solidbody.py b/src/felupe/mechanics/_solidbody.py index fa806639..e6d6d314 100644 --- a/src/felupe/mechanics/_solidbody.py +++ b/src/felupe/mechanics/_solidbody.py @@ -398,7 +398,6 @@ def _cauchy_stress(self, field=None): return dot(P, transpose(F)) / J def _mass(self, density=None): - if density is None: density = self.density diff --git a/src/felupe/tools/_hello_world.py b/src/felupe/tools/_hello_world.py index 120a8f0a..5986a345 100644 --- a/src/felupe/tools/_hello_world.py +++ b/src/felupe/tools/_hello_world.py @@ -101,7 +101,6 @@ def hello_world( plot = "\n".join(plot) if not container: - kwargs = ", ".join(kwargs) # fmt: off @@ -126,7 +125,6 @@ def hello_world( # fmt: on else: - kwargs.insert(0, "x0=field") kwargs = ", ".join(kwargs) diff --git a/tests/test_dof.py b/tests/test_dof.py index eb473f4d..1b18b2e2 100644 --- a/tests/test_dof.py +++ b/tests/test_dof.py @@ -145,6 +145,24 @@ def test_loadcase(): assert "top" in sh[0] +def test_boundary_multiaxial(): + mesh = fem.Rectangle(n=3) + region = fem.RegionQuad(mesh) + field = fem.FieldContainer([fem.FieldPlaneStrain(region, dim=2)]) + + for value in [1.0, np.arange(1, 3), np.arange(1, 3).reshape(1, 2)]: + boundaries = dict( + left=fem.Boundary(field[0], fx=0), + right=fem.Boundary(field[0], fx=1, value=value), + ) + + dof0, dof1 = fem.dof.partition(field, boundaries) + ext0 = fem.dof.apply(field, boundaries, dof0) + + assert ext0.shape == dof0.shape + + if __name__ == "__main__": test_boundary() + test_boundary_multiaxial() test_loadcase() diff --git a/tests/test_free_vibration.py b/tests/test_free_vibration.py index 83708194..54f69329 100644 --- a/tests/test_free_vibration.py +++ b/tests/test_free_vibration.py @@ -31,7 +31,6 @@ def test_free_vibration(): - mesh = fem.Cube(a=(0, 0, 5), b=(50, 100, 30), n=(3, 6, 4)) region = fem.RegionHexahedron(mesh) field = fem.FieldContainer([fem.Field(region, dim=3)]) @@ -44,7 +43,6 @@ def test_free_vibration(): def test_free_vibration_mixed(): - meshes = [ fem.Cube(a=(0, 0, 30), b=(50, 100, 35), n=(3, 6, 2)), fem.Cube(a=(0, 0, 5), b=(50, 100, 30), n=(3, 6, 4)), From 1b30023ae66013208b394f84b28e8198be9e7c73 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Sat, 30 Nov 2024 22:36:06 +0100 Subject: [PATCH 08/47] Examples: Simplify and remove `ravel()` from nd-prescribed `values` this is now handled automatically --- examples/ex05_rubber-metal-bushing.py | 2 +- examples/ex13_morph-rubber-wheel.py | 2 +- examples/ex17_torsion-gif.py | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/ex05_rubber-metal-bushing.py b/examples/ex05_rubber-metal-bushing.py index bd7996e6..8eabdd34 100644 --- a/examples/ex05_rubber-metal-bushing.py +++ b/examples/ex05_rubber-metal-bushing.py @@ -106,7 +106,7 @@ center=[0, 0, 0], )[0] inner_radial = 8 * np.array([0, 0, 1]) * progress - move.append((inner_radial + inner_rotated - inner).ravel()) + move.append(inner_radial + inner_rotated - inner) # %% # After defining the load step, the simulation model is ready to be solved. diff --git a/examples/ex13_morph-rubber-wheel.py b/examples/ex13_morph-rubber-wheel.py index 92424738..af8106b0 100644 --- a/examples/ex13_morph-rubber-wheel.py +++ b/examples/ex13_morph-rubber-wheel.py @@ -233,7 +233,7 @@ def sigmoid(x): axis=0, center=[0, 0], )[0] - move.append((center_rotated - center).ravel()) + move.append(center_rotated - center) # %% # A nearly-incompressible solid body is created for the rubber. At the bottom, a diff --git a/examples/ex17_torsion-gif.py b/examples/ex17_torsion-gif.py index 76b42061..47b4b818 100644 --- a/examples/ex17_torsion-gif.py +++ b/examples/ex17_torsion-gif.py @@ -47,7 +47,7 @@ axis=2, center=[0, 0, 1], )[0] - move.append((top_rotated - top).ravel()) + move.append(top_rotated - top) # %% # The reaction moment on the centerpoint of the right end face is tracked by a From 6b59077ead82c311bfa7e192caa7d433d1783e63 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Sun, 1 Dec 2024 15:30:26 +0100 Subject: [PATCH 09/47] Add `Boundary.plot()` and fix `Boundary(..., mode="and")` (#913) * Add `Boundary.plot()` * Update CHANGELOG.md * Update extut03_building_blocks.py --- CHANGELOG.md | 4 + .../examples/extut03_building_blocks.py | 4 + src/felupe/dof/_boundary.py | 108 +++++++++--------- tests/test_dof.py | 19 +++ 4 files changed, 79 insertions(+), 56 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 2bfea234..8a7528ce 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,7 @@ All notable changes to this project will be documented in this file. The format - Add `SolidBody.assemble.mass(density=1.0)` and `SolidBodyNearlyIncompressible.assemble.mass(density=1.0)` to assemble the mass matrix. - Add `SolidBody.evaluate.stress(field)` to evaluate the (first Piola-Kirchhoff) stress tensor (engineering stress in linear elasticity). - Add a free-vibration modal analysis Step/Job `FreeVibration(items, boundaries)` with methods to evaluate `FreeVibration.evaluate()` and to extract `field, frequency = FreeVibration.extract(n)` its n-th result. +- Add `Boundary.plot()` to plot the points and prescribed directions of a boundary. ### Changed - The first Piola-Kirchhoff stress tensor is evaluated if `ViewSolid(stress_type=None)`. @@ -14,6 +15,9 @@ All notable changes to this project will be documented in this file. The format - Enhance the `hello_world(axisymmetric=False, planestrain=False, curve=False, xdmf=False, container=False)` function with new arguments to customize the generated template script. - Enhance `Boundary` with added support for multiaxial prescribed values. +### Fixed +- Fix `Boundary(..., mode="and")` by ignoring any undefined axis. + ## [9.1.0] - 2024-11-23 ### Added diff --git a/docs/tutorial/examples/extut03_building_blocks.py b/docs/tutorial/examples/extut03_building_blocks.py index fb2d23e5..3a31347e 100644 --- a/docs/tutorial/examples/extut03_building_blocks.py +++ b/docs/tutorial/examples/extut03_building_blocks.py @@ -176,6 +176,10 @@ def W(C, mu, bulk): boundaries["right"] = fem.Boundary(displacement, fx=f1, skip=(1, 0, 0)) boundaries["move"] = fem.Boundary(displacement, fx=f1, skip=(0, 1, 1), value=0.5) +plotter = boundaries["left"].plot(color="green") +plotter = boundaries["right"].plot(plotter=plotter, color="red") +boundaries["move"].plot(plotter=plotter, color="blue", point_size=5).show() + # %% # Partition of deegrees of freedom # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ diff --git a/src/felupe/dof/_boundary.py b/src/felupe/dof/_boundary.py index acb288e6..b7561cb6 100644 --- a/src/felupe/dof/_boundary.py +++ b/src/felupe/dof/_boundary.py @@ -92,18 +92,10 @@ class Boundary: :context: :force_static: - >>> import pyvista as pv - >>> >>> right = fem.Boundary(displacement, fx=x.max()) >>> right = fem.Boundary(displacement, fx=lambda x: np.isclose(x, x.max())) >>> - >>> plotter = pv.Plotter() - >>> actor = plotter.add_points( - ... mesh.points[right.points], - ... point_size=20, - ... color="red", - ... ) - >>> mesh.plot(plotter=plotter, opacity=0.7).show() + >>> right.plot(color="red", plotter=mesh.plot(opacity=0.7)).show() If ``fx`` and ``fy`` are given, the masks are combined by *logical-or*. @@ -112,14 +104,7 @@ class Boundary: :force_static: >>> axes = fem.Boundary(displacement, fx=0, fy=0, mode="or") - >>> - >>> plotter = pv.Plotter() - >>> actor = plotter.add_points( - ... mesh.points[axes.points], - ... point_size=20, - ... color="red", - ... ) - >>> mesh.plot(plotter=plotter, opacity=0.7).show() + >>> axes.plot(color="red", plotter=mesh.plot(opacity=0.7)).show() This may be changed to *logical-and* if desired. @@ -128,14 +113,7 @@ class Boundary: :force_static: >>> center = fem.Boundary(displacement, fx=0, fy=0, mode="and") - >>> - >>> plotter = pv.Plotter() - >>> actor = plotter.add_points( - ... mesh.points[center.points], - ... point_size=20, - ... color="red", - ... ) - >>> mesh.plot(plotter=plotter, opacity=0.7).show() + >>> center.plot(color="red", plotter=mesh.plot(opacity=0.7)).show() For the most-general case, a user-defined boolean mask for the selection of the mesh-points is provided. While the two upper methods are useful to select @@ -149,13 +127,7 @@ class Boundary: >>> mask = np.logical_and(np.isclose(x**2 + y**2, 1), x >= 0) >>> surface = fem.Boundary(displacement, mask=mask) >>> - >>> plotter = pv.Plotter() - >>> actor = plotter.add_points( - ... mesh.points[surface.points], - ... point_size=20, - ... color="red", - ... ) - >>> mesh.plot(plotter=plotter, opacity=0.7).show() + >>> surface.plot(color="red", plotter=mesh.plot(opacity=0.7)).show() The application of a new mask allows to change the selected points of an existing boundary condition. @@ -167,13 +139,7 @@ class Boundary: >>> new_mask = np.logical_and(mask, y <= 0) >>> surface.apply_mask(new_mask) >>> - >>> plotter = pv.Plotter() - >>> actor = plotter.add_points( - ... mesh.points[surface.points], - ... point_size=20, - ... color="red", - ... ) - >>> mesh.plot(plotter=plotter, opacity=0.7).show() + >>> surface.plot(color="red", plotter=mesh.plot(opacity=0.7)).show() A boundary condition may be skipped on given axes, i.e. if only the x-components of a field should be prescribed on the selected points, then the y-axis must @@ -184,14 +150,7 @@ class Boundary: :force_static: >>> axes_x = fem.Boundary(displacement, fx=0, fy=0, skip=(False, True)) - >>> - >>> plotter = pv.Plotter() - >>> actor = plotter.add_points( - ... mesh.points[axes_x.points], - ... point_size=20, - ... color="red", - ... ) - >>> mesh.plot(plotter=plotter, opacity=0.7).show() + >>> axes_x.plot(color="red", plotter=mesh.plot(opacity=0.7)).show() Values for the prescribed degress of freedom are either applied during creation or by the update-method. @@ -203,13 +162,7 @@ class Boundary: >>> left = fem.Boundary(displacement, fx=x.min(), value=-0.2) >>> left.update(-0.3) >>> - >>> plotter = pv.Plotter() - >>> actor = plotter.add_points( - ... mesh.points[left.points], - ... point_size=20, - ... color="red", - ... ) - >>> mesh.plot(plotter=plotter, opacity=0.7).show() + >>> left.plot(color="red", plotter=mesh.plot(opacity=0.7)).show() Sometimes it is useful to create a boundary with all axes skipped. This boundary has no prescribed degrees of freedom and hence, is without effect. @@ -265,7 +218,7 @@ def __init__( # select the logical combination function "or" or "and" combine = {"or": np.logical_or, "and": np.logical_and}[self.mode] - mask = combine.reduce(masks) + mask = combine.reduce(np.array(masks)[[f != np.isnan for f in self.fun]]) self.apply_mask(mask) @@ -307,4 +260,47 @@ def apply_mask(self, mask): def update(self, value): "Update the value of the boundary in-place." - self.value = value + self.value = value # + + def plot( + self, plotter=None, color="black", scale=0.125, point_size=10, width=3, **kwargs + ): + "Plot the points and their prescribed directions of a boundary condition." + + mesh = self.field.region.mesh + + if plotter is None: + plotter = mesh.plot() + + if self.dim > 1 and self.dim != mesh.dim: + raise ValueError( + " ".join( + [ + "Plotting is not supported.", + "Field and Mesh must have equal dimensions.", + ] + ) + ) + + if len(self.points) > 0: + magnitude = min(mesh.points.max(axis=0) - mesh.points.min(axis=0)) * scale + + _ = plotter.add_points( + mesh.points[self.points], + color=color, + point_size=point_size, + label=self.name, + ) + + for skip, direction in zip(self.skip, np.eye(mesh.dim)): + if not skip: + end = mesh.points[self.points] + direction * magnitude + _ = plotter.add_lines( + np.hstack([mesh.points[self.points], end]).reshape( + -1, mesh.dim + ), + color=color, + width=width, + ) + + return plotter diff --git a/tests/test_dof.py b/tests/test_dof.py index 1b18b2e2..faa6372d 100644 --- a/tests/test_dof.py +++ b/tests/test_dof.py @@ -162,7 +162,26 @@ def test_boundary_multiaxial(): assert ext0.shape == dof0.shape +def test_boundary_plot(): + region = fem.RegionHexahedron(fem.Cube(b=(3, 1, 1), n=2)) + field = fem.Field(region, dim=3).as_container() + boundaries = dict( + left=fem.Boundary(field[0], fx=0, skip=(0, 1, 0)), + right=fem.Boundary(field[0], fx=1, skip=(0, 0, 0)), + ) + _ = boundaries["left"].plot(plotter=boundaries["right"].plot()) + + field = fem.Field(region, dim=2).as_container() + boundaries = dict( + left=fem.Boundary(field[0], fx=0, skip=(0, 1, 0)), + ) + + with pytest.raises(ValueError): + _ = boundaries["left"].plot() + + if __name__ == "__main__": test_boundary() test_boundary_multiaxial() + test_boundary_plot() test_loadcase() From 6c0afa24ab522b5f392ddd8c9b46ba72a1be04a5 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Mon, 2 Dec 2024 08:24:10 +0100 Subject: [PATCH 10/47] Add `BoundaryDict` (#915) * Add `BoundaryDict` with methods to plot, screenshot and imshow * Update test_dof.py * Update test_dof.py * Ensure 3d-data in `Boundary.plot()` for PyVista * Update test_dof.py * Add `Boundary.plot(..., show_points=True, show_lines=True)` --- CHANGELOG.md | 1 + docs/felupe/dof.rst | 6 ++ .../examples/extut03_building_blocks.py | 7 +- src/felupe/__init__.py | 3 +- src/felupe/dof/__init__.py | 2 + src/felupe/dof/_boundary.py | 53 ++++++++----- src/felupe/dof/_dict.py | 79 +++++++++++++++++++ tests/test_dof.py | 14 ++++ 8 files changed, 142 insertions(+), 23 deletions(-) create mode 100644 src/felupe/dof/_dict.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 8a7528ce..10b3767d 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 `SolidBody.evaluate.stress(field)` to evaluate the (first Piola-Kirchhoff) stress tensor (engineering stress in linear elasticity). - Add a free-vibration modal analysis Step/Job `FreeVibration(items, boundaries)` with methods to evaluate `FreeVibration.evaluate()` and to extract `field, frequency = FreeVibration.extract(n)` its n-th result. - Add `Boundary.plot()` to plot the points and prescribed directions of a boundary. +- Add `BoundaryDict` as a subclassed dict with methods to `plot()`, `screenshot()` and `imshow()`. ### Changed - The first Piola-Kirchhoff stress tensor is evaluated if `ViewSolid(stress_type=None)`. diff --git a/docs/felupe/dof.rst b/docs/felupe/dof.rst index 359691c6..a41ebedf 100644 --- a/docs/felupe/dof.rst +++ b/docs/felupe/dof.rst @@ -10,6 +10,7 @@ This module contains the definition of a boundary condition, tools related to th .. autosummary:: Boundary + BoundaryDict **Tools** @@ -37,6 +38,11 @@ This module contains the definition of a boundary condition, tools related to th :undoc-members: :inherited-members: +.. autoclass:: felupe.BoundaryDict + :members: + :undoc-members: + :inherited-members: + .. autofunction:: felupe.dof.partition .. autofunction:: felupe.dof.apply diff --git a/docs/tutorial/examples/extut03_building_blocks.py b/docs/tutorial/examples/extut03_building_blocks.py index 3a31347e..fb5c06da 100644 --- a/docs/tutorial/examples/extut03_building_blocks.py +++ b/docs/tutorial/examples/extut03_building_blocks.py @@ -171,14 +171,13 @@ def W(C, mu, bulk): f0 = lambda x: np.isclose(x, 0) f1 = lambda x: np.isclose(x, 1) -boundaries = {} +boundaries = fem.BoundaryDict() + boundaries["left"] = fem.Boundary(displacement, fx=f0) boundaries["right"] = fem.Boundary(displacement, fx=f1, skip=(1, 0, 0)) boundaries["move"] = fem.Boundary(displacement, fx=f1, skip=(0, 1, 1), value=0.5) -plotter = boundaries["left"].plot(color="green") -plotter = boundaries["right"].plot(plotter=plotter, color="red") -boundaries["move"].plot(plotter=plotter, color="blue", point_size=5).show() +boundaries.plot().show() # %% # Partition of deegrees of freedom diff --git a/src/felupe/__init__.py b/src/felupe/__init__.py index dd1c2fcd..4d4fb7af 100644 --- a/src/felupe/__init__.py +++ b/src/felupe/__init__.py @@ -40,7 +40,7 @@ linear_elastic, linear_elastic_plastic_isotropic_hardening, ) -from .dof import Boundary +from .dof import Boundary, BoundaryDict from .element import ArbitraryOrderLagrange as ArbitraryOrderLagrangeElement from .element import ( BiQuadraticQuad, @@ -168,6 +168,7 @@ "linear_elastic", "linear_elastic_plastic_isotropic_hardening", "Boundary", + "BoundaryDict", "ArbitraryOrderLagrangeElement", "BiQuadraticQuad", "ConstantHexahedron", diff --git a/src/felupe/dof/__init__.py b/src/felupe/dof/__init__.py index 9444a0a6..e5ef8e0a 100644 --- a/src/felupe/dof/__init__.py +++ b/src/felupe/dof/__init__.py @@ -1,9 +1,11 @@ from ._boundary import Boundary +from ._dict import BoundaryDict from ._loadcase import biaxial, shear, symmetry, uniaxial from ._tools import apply, get_dof0, get_dof1, partition __all__ = [ "Boundary", + "BoundaryDict", "biaxial", "shear", "symmetry", diff --git a/src/felupe/dof/_boundary.py b/src/felupe/dof/_boundary.py index b7561cb6..5e2ade69 100644 --- a/src/felupe/dof/_boundary.py +++ b/src/felupe/dof/_boundary.py @@ -263,7 +263,16 @@ def update(self, value): self.value = value # def plot( - self, plotter=None, color="black", scale=0.125, point_size=10, width=3, **kwargs + self, + plotter=None, + color="black", + scale=0.125, + point_size=10, + width=3, + label=None, + show_points=True, + show_lines=True, + **kwargs, ): "Plot the points and their prescribed directions of a boundary condition." @@ -283,24 +292,32 @@ def plot( ) if len(self.points) > 0: - magnitude = min(mesh.points.max(axis=0) - mesh.points.min(axis=0)) * scale + if label is None: + label = self.name + + if show_points or show_lines: + points = np.pad(mesh.points, ((0, 0), (0, 3 - mesh.dim))) + + if show_points: + _ = plotter.add_points( + points[self.points], + color=color, + point_size=point_size, + label=label, + ) - _ = plotter.add_points( - mesh.points[self.points], - color=color, - point_size=point_size, - label=self.name, - ) + if show_lines: + magnitude = ( + min(mesh.points.max(axis=0) - mesh.points.min(axis=0)) * scale + ) - for skip, direction in zip(self.skip, np.eye(mesh.dim)): - if not skip: - end = mesh.points[self.points] + direction * magnitude - _ = plotter.add_lines( - np.hstack([mesh.points[self.points], end]).reshape( - -1, mesh.dim - ), - color=color, - width=width, - ) + for skip, direction in zip(self.skip, np.eye(3)): + if not skip: + end = points[self.points] + direction * magnitude + _ = plotter.add_lines( + np.hstack([points[self.points], end]).reshape(-1, 3), + color=color, + width=width, + ) return plotter diff --git a/src/felupe/dof/_dict.py b/src/felupe/dof/_dict.py new file mode 100644 index 00000000..f99381cd --- /dev/null +++ b/src/felupe/dof/_dict.py @@ -0,0 +1,79 @@ +class BoundaryDict(dict): + "A dict of boundary conditions." + + def plot( + self, + plotter=None, + colors=None, + size=(0.1, 0.1), + show_points=True, + show_lines=True, + **kwargs, + ): + "Plot the boundary conditions." + + if colors is None: + import matplotlib.colors as mcolors + + colors = list(mcolors.TABLEAU_COLORS.values()) + + colors_list = [] + while len(colors_list) < len(self.keys()): + colors_list = [*colors_list, *colors] + + for (key, boundary), color in zip(self.items(), colors_list): + label = key + if boundary.name != "default": + label = boundary.name + + plotter = boundary.plot( + label=label, + color=color, + plotter=plotter, + show_points=show_points, + show_lines=show_lines, + **kwargs, + ) + + plotter.add_legend(size=size, bcolor="white") + + return plotter + + def screenshot( + self, + *args, + filename="boundaries.png", + transparent_background=None, + scale=None, + colors=None, + plotter=None, + **kwargs, + ): + "Take a screenshot of the boundary conditions." + + if plotter is None: + mesh = self[list(self.keys())[0]].field.region.mesh + plotter = mesh.plot(off_screen=True) + + plotter = self.plot(plotter=plotter, colors=colors, **kwargs) + + return plotter.screenshot( + filename=filename, + transparent_background=transparent_background, + scale=scale, + ) + + def imshow(self, *args, ax=None, dpi=None, **kwargs): + """Take a screenshot of the boundary conditions, show the image data in + a figure and return the ax. + """ + + if ax is None: + import matplotlib.pyplot as plt + + fig, ax = plt.subplots(dpi=dpi) + + ax.imshow(self.screenshot(*args, filename=None, **kwargs)) + ax.set_axis_off() + + return ax diff --git a/tests/test_dof.py b/tests/test_dof.py index faa6372d..ba61ee7d 100644 --- a/tests/test_dof.py +++ b/tests/test_dof.py @@ -180,8 +180,22 @@ def test_boundary_plot(): _ = boundaries["left"].plot() +def test_boundary_dict(): + field = fem.FieldPlaneStrain( + fem.RegionQuad(fem.Rectangle(b=(3, 1), n=5)), dim=2 + ).as_container() + boundaries = fem.BoundaryDict( + left=fem.Boundary(field[0], fx=0, skip=(0, 0)), + right=fem.Boundary(field[0], name="my_label", fx=3, skip=(0, 1)), + ) + plotter = boundaries.plot() + # img = boundaries.screenshot() + # ax = boundaries.imshow() + + if __name__ == "__main__": test_boundary() + test_boundary_dict() test_boundary_multiaxial() test_boundary_plot() test_loadcase() From 1985d46c205fad7e98d76e80a61b47a6b3b012e0 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Mon, 2 Dec 2024 08:27:10 +0100 Subject: [PATCH 11/47] Update ex10_poisson-equation.py --- examples/ex10_poisson-equation.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/examples/ex10_poisson-equation.py b/examples/ex10_poisson-equation.py index 1925c3f9..c95db314 100644 --- a/examples/ex10_poisson-equation.py +++ b/examples/ex10_poisson-equation.py @@ -28,7 +28,7 @@ \ d\Omega = \int_\Omega \delta u \cdot f \ d\Omega """ - +# sphinx_gallery_thumbnail_number = -1 import felupe as fem mesh = fem.Rectangle(n=2**5).triangulate() @@ -36,12 +36,13 @@ u = fem.Field(region, dim=1) field = fem.FieldContainer([u]) -boundaries = dict( +boundaries = fem.BoundaryDict( bottom=fem.Boundary(u, fy=0), top=fem.Boundary(u, fy=1), left=fem.Boundary(u, fx=0), right=fem.Boundary(u, fx=1), ) +boundaries.plot(show_lines=False).show() solid = fem.SolidBody(umat=fem.Laplace(), field=field) load = fem.SolidBodyForce(field=field, values=1.0) From 139383dbc844e41735d2a7598617f31e53131f54 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Mon, 2 Dec 2024 08:42:27 +0100 Subject: [PATCH 12/47] Update _loadcase.py --- src/felupe/dof/_loadcase.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/felupe/dof/_loadcase.py b/src/felupe/dof/_loadcase.py index 33b19e56..eebacad7 100644 --- a/src/felupe/dof/_loadcase.py +++ b/src/felupe/dof/_loadcase.py @@ -19,6 +19,7 @@ import numpy as np from ._boundary import Boundary +from ._dict import BoundaryDict from ._tools import apply, partition @@ -108,7 +109,7 @@ def symmetry(field, axes=(True, True, True), x=0.0, y=0.0, z=0.0, bounds=None): ] if bounds is None: - bounds = {} + bounds = BoundaryDict() labels = ["symx", "symy", "symz"] # loop over symmetry conditions and add them to a new dict @@ -559,7 +560,7 @@ def shear( bounds = symmetry(f, axes=sym) else: - bounds = {} + bounds = BoundaryDict() fy = ["fx", "fy", "fz"][axes[1]] From ada2a2f8528dee36b13edcee666a8077709e6a81 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Wed, 4 Dec 2024 17:12:38 +0100 Subject: [PATCH 13/47] Add `PointLoad.plot()` (#916) * Add `PointLoad.plot()` * Update _pointload.py * Update test_mechanics.py --- src/felupe/mechanics/_pointload.py | 39 ++++++++++++++++++++++++++++++ tests/test_mechanics.py | 7 +++++- 2 files changed, 45 insertions(+), 1 deletion(-) diff --git a/src/felupe/mechanics/_pointload.py b/src/felupe/mechanics/_pointload.py index e5f825b5..389a60df 100644 --- a/src/felupe/mechanics/_pointload.py +++ b/src/felupe/mechanics/_pointload.py @@ -109,3 +109,42 @@ def _matrix(self, field=None, parallel=False): self.results.stiffness = csr_matrix(([0.0], ([0], [0])), shape=(n, n)) return self.results.stiffness + + def plot( + self, + plotter=None, + color="red", + scale=0.125, + **kwargs, + ): + "Plot the point load." + + mesh = self.field.region.mesh + + if plotter is None: + plotter = mesh.plot() + + if len(self.points) > 0: + points = np.pad(mesh.points, ((0, 0), (0, 3 - mesh.dim))) + magnitude = min(mesh.points.max(axis=0) - mesh.points.min(axis=0)) * scale + + values = np.atleast_2d(self.values) + + skip = np.zeros(3, dtype=bool) + skip[values.shape[1] :] = True + + if values.shape[1] > 1: + 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) + if not skip_axis: + _ = plotter.add_arrows( + points[self.points], + direction=d * np.sign(values[0, a]), + mag=magnitude, + color=color, + **kwargs, + ) + + return plotter diff --git a/tests/test_mechanics.py b/tests/test_mechanics.py index fdf3daa8..5d0d2c0a 100644 --- a/tests/test_mechanics.py +++ b/tests/test_mechanics.py @@ -503,8 +503,8 @@ def test_load(): values *= 0.025 body = fem.SolidBody(umat, field) - load = fem.PointLoad(field, mask, values=values, axisymmetric=axi) + load = fem.PointLoad(field, mask, values=values, axisymmetric=axi) bounds = {"fix": fem.Boundary(field[0], fx=lambda x: x == 0)} dof0, dof1 = fem.dof.partition(field, bounds) @@ -512,6 +512,11 @@ def test_load(): assert res.success + plotter = load.plot(plotter=field.region.mesh.plot(off_screen=True)) + + load.update([0, 0.1]) + plotter = load.plot(plotter=field.region.mesh.plot(off_screen=True)) + def test_view(): mesh = fem.Rectangle(n=3) From be2f25cacc7689b17665fcbddaab4b0856bfaf98 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Sun, 8 Dec 2024 02:10:36 +0100 Subject: [PATCH 14/47] Enhance `math.linsteps(..., values=0.0)` (#917) * Enhance `math.linsteps(..., values=0.0)` * `BoundaryDict().plot()` pass missing `size` arg * Use `BoundaryDict` in examples --- CHANGELOG.md | 1 + examples/ex01_beam.py | 4 +++- examples/ex02_plate-with-hole.py | 3 +-- examples/ex03_plasticity.py | 4 +++- examples/ex05_rubber-metal-bushing.py | 11 +++++------ examples/ex07_engine-mount.py | 5 +++-- examples/ex08_shear.py | 10 +++++----- src/felupe/dof/_dict.py | 1 + src/felupe/math/_math.py | 25 ++++++++++++++----------- 9 files changed, 36 insertions(+), 28 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 10b3767d..e9b4c5e9 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -15,6 +15,7 @@ All notable changes to this project will be documented in this file. The format - Autodetect the stress-type in `SolidBody.plot(name)` from `name`. - Enhance the `hello_world(axisymmetric=False, planestrain=False, curve=False, xdmf=False, container=False)` function with new arguments to customize the generated template script. - Enhance `Boundary` with added support for multiaxial prescribed values. +- Enhance `math.linsteps(..., values=0)` with default values except for the column `axis` if `axis` is not None. ### Fixed - Fix `Boundary(..., mode="and")` by ignoring any undefined axis. diff --git a/examples/ex01_beam.py b/examples/ex01_beam.py index 67fd84a7..2821f487 100644 --- a/examples/ex01_beam.py +++ b/examples/ex01_beam.py @@ -20,6 +20,7 @@ dimensional vector-valued displacement field is initiated on the region. """ +# sphinx_gallery_thumbnail_number = -1 import felupe as fem cube = fem.Cube(a=(0, 0, 0), b=(2000, 100, 100), n=(101, 6, 6)) @@ -29,7 +30,8 @@ # %% # A fixed boundary condition is applied on the left end of the beam. -boundaries = {"fixed": fem.dof.Boundary(displacement, fx=0)} +boundaries = fem.BoundaryDict(fixed=fem.dof.Boundary(displacement, fx=0)) +boundaries.plot().show() # %% # The material behaviour is defined through a built-in isotropic linear-elastic material diff --git a/examples/ex02_plate-with-hole.py b/examples/ex02_plate-with-hole.py index e4b600c8..5dd16265 100644 --- a/examples/ex02_plate-with-hole.py +++ b/examples/ex02_plate-with-hole.py @@ -47,8 +47,6 @@ mesh = fem.mesh.concatenate([face, face.mirror(normal=[-1, 1, 0]), rect]) mesh = mesh.sweep(decimals=5) -mesh.plot().show() - # %% # A numeric quad-region created on the mesh in combination with a vector-valued # displacement field represents the plate. The Boundary conditions for the symmetry @@ -58,6 +56,7 @@ field = fem.FieldContainer([displacement]) boundaries = fem.dof.symmetry(displacement) +boundaries.plot().show() # %% # The material behaviour is defined through a built-in isotropic linear-elastic material diff --git a/examples/ex03_plasticity.py b/examples/ex03_plasticity.py index 13e94871..c40e0392 100644 --- a/examples/ex03_plasticity.py +++ b/examples/ex03_plasticity.py @@ -23,6 +23,7 @@ numeric region. """ +# sphinx_gallery_thumbnail_number = -1 import numpy as np import felupe as fem @@ -34,7 +35,8 @@ # %% # A fixed boundary condition is applied at :math:`x=0`. -boundaries = {"fixed": fem.dof.Boundary(displacement, fx=0)} +boundaries = fem.BoundaryDict(fixed=fem.dof.Boundary(displacement, fx=0)) +boundaries.plot().show() # %% # The material behaviour is defined through a built-in isotropic linear-elastic plastic diff --git a/examples/ex05_rubber-metal-bushing.py b/examples/ex05_rubber-metal-bushing.py index 8eabdd34..100640b3 100644 --- a/examples/ex05_rubber-metal-bushing.py +++ b/examples/ex05_rubber-metal-bushing.py @@ -57,8 +57,6 @@ mesh = fem.mesh.stack(meshes.meshes) x, y, z = mesh.points.T -mesh.plot().show() - # %% # A global region as well as sub-regions for all materials are generated. The same # applies to the fields, the material formulations as well as the solid bodies. @@ -78,10 +76,11 @@ # %% # The boundary conditions are created on the global displacement field. Masks are # created for both the innermost and the outermost metal sheet faces. -boundaries = { - "inner": fem.dof.Boundary(field[0], mask=np.isclose(np.sqrt(y**2 + z**2), 25)), - "outer": fem.dof.Boundary(field[0], mask=np.isclose(np.sqrt(y**2 + z**2), 65)), -} +boundaries = fem.BoundaryDict( + inner=fem.dof.Boundary(field[0], mask=np.isclose(np.sqrt(y**2 + z**2), 25)), + outer=fem.dof.Boundary(field[0], mask=np.isclose(np.sqrt(y**2 + z**2), 65)), +) +boundaries.plot(show_lines=False).show() # prescribed values for the innermost radial mesh points table = fem.math.linsteps([0, 1], num=3) diff --git a/examples/ex07_engine-mount.py b/examples/ex07_engine-mount.py index 6ea405a6..93e97aa6 100644 --- a/examples/ex07_engine-mount.py +++ b/examples/ex07_engine-mount.py @@ -42,7 +42,6 @@ # sub-meshes with shared points-array and a global mesh meshes = fem.MeshContainer([metal, rubber, air], merge=True) mesh = fem.mesh.concatenate(meshes).sweep() -meshes.plot(colors=["grey", "black", "white"]).show() # %% # A global region as well as sub-regions for all materials are generated. The same @@ -65,11 +64,13 @@ inner = np.logical_and(only_cells_metal, radius <= 45) outer = np.logical_and(only_cells_metal, radius > 45) -boundaries = dict( +boundaries = fem.BoundaryDict( fixed=fem.Boundary(field[0], mask=outer), u_x=fem.Boundary(field[0], mask=inner, skip=(0, 1)), u_y=fem.Boundary(field[0], mask=inner, skip=(1, 0)), ) +plotter = meshes.plot(colors=["grey", "black", "white"]) +boundaries.plot(plotter=plotter, scale=0.02).show() # %% # The material behaviour of the rubberlike solid is defined through a built-in hyperelastic diff --git a/examples/ex08_shear.py b/examples/ex08_shear.py index f30afa2b..31fdfbb9 100644 --- a/examples/ex08_shear.py +++ b/examples/ex08_shear.py @@ -58,10 +58,10 @@ region = fem.RegionQuad(mesh) field = fem.FieldsMixed(region, n=3, planestrain=True) -boundaries = { - "fixed": fem.Boundary(field[0], fy=mesh.y.min()), - "control": fem.Boundary(field[0], fy=mesh.y.max(), skip=(0, 1)), -} +boundaries = fem.BoundaryDict( + fixed=fem.Boundary(field[0], fy=mesh.y.min()), + control=fem.Boundary(field[0], fy=mesh.y.max(), skip=(0, 1)), +) dof0, dof1 = fem.dof.partition(field, boundaries) @@ -97,7 +97,7 @@ centerpoint=mesh.npoints - 1, ) -plotter = mesh.plot() +plotter = boundaries.plot() plotter = mpc.plot(plotter=plotter) plotter.show() diff --git a/src/felupe/dof/_dict.py b/src/felupe/dof/_dict.py index f99381cd..d1fba786 100644 --- a/src/felupe/dof/_dict.py +++ b/src/felupe/dof/_dict.py @@ -30,6 +30,7 @@ def plot( label=label, color=color, plotter=plotter, + size=size, show_points=show_points, show_lines=show_lines, **kwargs, diff --git a/src/felupe/math/_math.py b/src/felupe/math/_math.py index 3bc4f869..d27ff646 100644 --- a/src/felupe/math/_math.py +++ b/src/felupe/math/_math.py @@ -19,7 +19,7 @@ import numpy as np -def linsteps(points, num=10, endpoint=True, axis=None, axes=None): +def linsteps(points, num=10, endpoint=True, axis=None, axes=None, values=0.0): """Return a sequence from batches of evenly spaced samples between pairs of milestone points. @@ -40,6 +40,9 @@ def linsteps(points, num=10, endpoint=True, axis=None, axes=None): axes : int or None, optional Number of output columns if ``axis`` is not None. Requires ``axes > axis`` and will be set to ``axes = axis + 1`` by default (None). + values : array_like, optional + Default values if ``axis`` is not None (not used for column ``axis``). Default + is 0.0. Returns ------- @@ -49,7 +52,7 @@ def linsteps(points, num=10, endpoint=True, axis=None, axes=None): Examples -------- >>> import felupe as fem - + >>> >>> fem.math.linsteps([0, 0.5, 1.5, 3.5], num=2) array([0. , 0.25, 0.5 , 1. , 1.5 , 2.5 , 3.5 ]) @@ -58,14 +61,14 @@ def linsteps(points, num=10, endpoint=True, axis=None, axes=None): >>> fem.math.linsteps([0, 0.5, 1.5, 3.5], num=2, endpoint=False) array([0. , 0.25, 0.5 , 1. , 1.5 , 2.5 ]) - >>> fem.math.linsteps([0, 0.5, 1.5, 3.5], num=2, axis=1) - array([[0. , 0. ], - [0. , 0.25], - [0. , 0.5 ], - [0. , 1. ], - [0. , 1.5 ], - [0. , 2.5 ], - [0. , 3.5 ]]) + >>> fem.math.linsteps([0, 0.5, 1.5, 3.5], num=2, axis=1, values=[-1, 0]) + array([[-1. , 0. ], + [-1. , 0.25], + [-1. , 0.5 ], + [-1. , 1. ], + [-1. , 1.5 ], + [-1. , 2.5 ], + [-1. , 3.5 ]]) Output with four columns: @@ -120,7 +123,7 @@ def linsteps(points, num=10, endpoint=True, axis=None, axes=None): axes = axis + 1 steps_1d = steps - steps = np.zeros((len(steps_1d), axes)) + steps = np.ones((len(steps_1d), axes)) * np.atleast_2d(values) steps[:, axis] = steps_1d return steps From e28bf5b13afe622ad0a7155698e3777df6490f4b Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Tue, 10 Dec 2024 16:19:28 +0100 Subject: [PATCH 15/47] Enhance `FieldContainer.link()` (#918) --- CHANGELOG.md | 1 + src/felupe/field/_container.py | 10 +++++++--- tests/test_field.py | 17 +++++++++++++++++ 3 files changed, 25 insertions(+), 3 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index e9b4c5e9..138602fd 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 - Enhance the `hello_world(axisymmetric=False, planestrain=False, curve=False, xdmf=False, container=False)` function with new arguments to customize the generated template script. - Enhance `Boundary` with added support for multiaxial prescribed values. - Enhance `math.linsteps(..., values=0)` with default values except for the column `axis` if `axis` is not None. +- Link all field-values to the values of the first field if no other field is given in `FieldContainer.link()`. ### Fixed - Fix `Boundary(..., mode="and")` by ignoring any undefined axis. diff --git a/src/felupe/field/_container.py b/src/felupe/field/_container.py index ffc15d1d..c85c02b0 100644 --- a/src/felupe/field/_container.py +++ b/src/felupe/field/_container.py @@ -210,10 +210,14 @@ def copy(self): "Return a copy of the field." return deepcopy(self) - def link(self, other_field): + def link(self, other_field=None): "Link value array of other field." - for field, newfield in zip(self.fields, other_field.fields): - field.values = newfield.values + if other_field is None: + for u in self.fields[1:]: # link n-to-1 (all-to-first) + u.values = self.fields[0].values + else: + for field, newfield in zip(self.fields, other_field.fields): # link 1-to-1 + field.values = newfield.values def view(self, point_data=None, cell_data=None, cell_type=None, project=None): """View the field with optional given dicts of point- and cell-data items. diff --git a/tests/test_field.py b/tests/test_field.py index d6a27d86..7b6d0135 100644 --- a/tests/test_field.py +++ b/tests/test_field.py @@ -247,9 +247,26 @@ def test_view(): # ax = mesh.imshow() +def test_link(): + mesh = fem.Cube(n=2) + region = fem.RegionHexahedron(mesh) + field = fem.Field(region, dim=3) & fem.Field(region, dim=1) + field.link() + + assert field[0].values is field[1].values + + field1 = fem.Field(region, dim=3) & fem.Field(region, dim=1) + field2 = fem.Field(region, dim=3) & fem.Field(region, dim=1) + field1.link(field2) + + assert field1[0].values is field2[0].values + assert field1[1].values is field2[1].values + + if __name__ == "__main__": test_axi() test_3d() test_3d_mixed() test_mixed_lagrange() test_view() + test_link() From 8607b63bd0cf29a51d6d4ca41430bc405cb5c993 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Wed, 11 Dec 2024 09:32:03 +0100 Subject: [PATCH 16/47] Update _tensor.py --- src/felupe/math/_tensor.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/felupe/math/_tensor.py b/src/felupe/math/_tensor.py index 2d4cd151..e6791408 100644 --- a/src/felupe/math/_tensor.py +++ b/src/felupe/math/_tensor.py @@ -401,7 +401,7 @@ def inv(A, determinant=None, full_output=False, sym=False, out=None): " ".join( [ "Wrong shape of first two axes.", - "Must be (1, 1), (2, 2) or (3, 3) but {A.shape[:2]} is given.", + f"Must be (1, 1), (2, 2) or (3, 3) but {A.shape[:2]} is given.", ] ) ) From a666ecca6e9df0dcc1191252906bf0f3921c6871 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Tue, 17 Dec 2024 16:05:40 +0100 Subject: [PATCH 17/47] Add `SolidBody(..., block=True, apply=None)` (#920) * Add `SolidBody(..., block=True, apply=None)` to apply a final callable in `SolidBody.assemble.vector()` or `SolidBody.assemble.matrix()`. * Add `SolidBody(..., multiplier=None)` * Update test_mechanics.py --- CHANGELOG.md | 2 + src/felupe/assembly/_integral.py | 7 ++-- src/felupe/mechanics/_solidbody.py | 61 ++++++++++++++++++++++++++++-- tests/test_mechanics.py | 10 ++++- 4 files changed, 70 insertions(+), 10 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 138602fd..56240aae 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 - Add a free-vibration modal analysis Step/Job `FreeVibration(items, boundaries)` with methods to evaluate `FreeVibration.evaluate()` and to extract `field, frequency = FreeVibration.extract(n)` its n-th result. - Add `Boundary.plot()` to plot the points and prescribed directions of a boundary. - Add `BoundaryDict` as a subclassed dict with methods to `plot()`, `screenshot()` and `imshow()`. +- Add a new argument to apply a callable on the assembled vector/matrix of a solid body, `SolidBody(..., apply=None)`. This may be used to sum the list of sub-blocks instead of stacking them together, `SolidBody(..., block=False, apply=None)`. This is useful for mixed formulations where both the deformation gradient and the displacement values are required. ### Changed - The first Piola-Kirchhoff stress tensor is evaluated if `ViewSolid(stress_type=None)`. @@ -17,6 +18,7 @@ All notable changes to this project will be documented in this file. The format - Enhance `Boundary` with added support for multiaxial prescribed values. - Enhance `math.linsteps(..., values=0)` with default values except for the column `axis` if `axis` is not None. - Link all field-values to the values of the first field if no other field is given in `FieldContainer.link()`. +- Change the default arguments from `block=True` to `block=None` in `SolidBody.assemble.vector(block=None)` and `SolidBody.assemble.matrix(block=None)` and define `block` on creation, `SolidBody(..., block=True)` instead. ### Fixed - Fix `Boundary(..., mode="and")` by ignoring any undefined axis. diff --git a/src/felupe/assembly/_integral.py b/src/felupe/assembly/_integral.py index ab97717b..a46aef93 100644 --- a/src/felupe/assembly/_integral.py +++ b/src/felupe/assembly/_integral.py @@ -273,13 +273,12 @@ def assemble(self, values=None, parallel=False, block=True, out=None): if i != j: K[j, i] = res[a].T - return bmat(K).tocsr() + res = bmat(K).tocsr() if block and self.mode == 1: - return vstack(res).tocsr() + res = vstack(res).tocsr() - else: - return res + return res def integrate(self, parallel=False, out=None): if out is None: diff --git a/src/felupe/mechanics/_solidbody.py b/src/felupe/mechanics/_solidbody.py index e6d6d314..6824ed75 100644 --- a/src/felupe/mechanics/_solidbody.py +++ b/src/felupe/mechanics/_solidbody.py @@ -164,6 +164,13 @@ class SolidBody(Solid): Array of initial internal state variables (default is None). density : float or None, optional The density of the solid body. + block : bool, optional + Assemble a sparse matrix from sparse sub-blocks or assemble a sparse vector by + stacking sparse matrices vertically (row wise). Default is True. + apply : callable or None, optional + Apply a callable on the assembled vectors and sparse matrices. Default is None. + multiplier : float or None, optional + A scale factor for the assembled vector and matrix. Default is None. Notes ----- @@ -241,10 +248,21 @@ class SolidBody(Solid): methods for the assembly of sparse vectors/matrices. """ - def __init__(self, umat, field, statevars=None, density=None): + def __init__( + self, + umat, + field, + statevars=None, + density=None, + block=True, + apply=None, + multiplier=None, + ): self.umat = umat self.field = field self.density = density + self.block = block + self.apply = apply self.results = Results(stress=True, elasticity=True) self.results.kinematics = self._extract(self.field) @@ -264,7 +282,10 @@ def __init__(self, umat, field, statevars=None, density=None): ) self.assemble = Assemble( - vector=self._vector, matrix=self._matrix, mass=self._mass + vector=self._vector, + matrix=self._matrix, + mass=self._mass, + multiplier=multiplier, ) self.evaluate = Evaluate( @@ -279,7 +300,14 @@ def __init__(self, umat, field, statevars=None, density=None): self._form = IntegralForm def _vector( - self, field=None, parallel=False, items=None, args=(), kwargs=None, block=True + self, + field=None, + parallel=False, + items=None, + args=(), + kwargs=None, + block=None, + apply=None, ): if kwargs is None: kwargs = {} @@ -287,6 +315,12 @@ def _vector( if field is not None: self.field = field + if block is None: + block = self.block + + if apply is None: + apply = self.apply + self.results.stress = self._gradient(field, args=args, kwargs=kwargs) self.results.force = self._form( fun=self.results.stress[slice(items)], @@ -294,10 +328,20 @@ def _vector( dV=self.field.region.dV, ).assemble(parallel=parallel, block=block) + if apply is not None: + self.results.force = apply(self.results.force) + return self.results.force def _matrix( - self, field=None, parallel=False, items=None, args=(), kwargs=None, block=True + self, + field=None, + parallel=False, + items=None, + args=(), + kwargs=None, + block=None, + apply=None, ): if kwargs is None: kwargs = {} @@ -305,6 +349,12 @@ def _matrix( if field is not None: self.field = field + if block is None: + block = self.block + + if apply is None: + apply = self.apply + self.results.elasticity = self._hessian(field, args=args, kwargs=kwargs) form = self._form( @@ -322,6 +372,9 @@ def _matrix( values=self.results.stiffness_values, block=block ) + if apply is not None: + self.results.stiffness = apply(self.results.stiffness) + return self.results.stiffness def _extract(self, field): diff --git a/tests/test_mechanics.py b/tests/test_mechanics.py index 5d0d2c0a..6b5e0072 100644 --- a/tests/test_mechanics.py +++ b/tests/test_mechanics.py @@ -293,7 +293,7 @@ def test_solidbody_incompressible(): assert np.allclose(t1, t2) -def test_solidbody_axi(): +def test_solidbody_axi_incompressible(): umat, u = pre_axi(bulk=None) b = fem.SolidBodyNearlyIncompressible(umat=umat, field=u, bulk=5000) b = fem.SolidBodyNearlyIncompressible( @@ -338,7 +338,7 @@ def test_solidbody_axi(): assert np.allclose(t1, t2) -def test_solidbody_axi_incompressible(): +def test_solidbody_axi(): umat, u = pre_axi() b = fem.SolidBody(umat=umat, field=u) @@ -351,12 +351,18 @@ def test_solidbody_axi_incompressible(): r2 = b.assemble.vector(**kwargs) assert np.allclose(r1.toarray(), r2.toarray()) + r3 = b.assemble.vector(**kwargs, block=False, apply=sum) + assert np.allclose(r1.toarray(), r3.toarray()) + K1 = b.assemble.matrix(u, **kwargs) assert K1.shape == (18, 18) K2 = b.assemble.matrix(**kwargs) assert np.allclose(K1.toarray(), K2.toarray()) + K3 = b.assemble.matrix(**kwargs, block=False, apply=sum) + assert np.allclose(K1.toarray(), K3.toarray()) + P1 = b.results.stress P2 = b.evaluate.gradient() P2 = b.evaluate.gradient(u) From 67e10df1918508d9ae2f1e4ca0baf4742326f053 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Wed, 18 Dec 2024 14:48:08 +0100 Subject: [PATCH 18/47] Enhance `IntegralForm`: Non-symmetric mixed forms and `None` (#922) * Enhance `IntegralForm`: Non-symmetric mixed forms and `None` * Update test_mechanics.py * Add more tests --- CHANGELOG.md | 2 ++ src/felupe/assembly/_cartesian.py | 32 +++++++++++++++++++------------ src/felupe/assembly/_integral.py | 27 +++++++++++++++++++++----- src/felupe/constitution/_mixed.py | 6 +++--- tests/test_form.py | 14 ++++++++++++++ tests/test_mechanics.py | 5 ++++- 6 files changed, 65 insertions(+), 21 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 56240aae..5e47b42d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,6 +10,7 @@ All notable changes to this project will be documented in this file. The format - Add `Boundary.plot()` to plot the points and prescribed directions of a boundary. - Add `BoundaryDict` as a subclassed dict with methods to `plot()`, `screenshot()` and `imshow()`. - Add a new argument to apply a callable on the assembled vector/matrix of a solid body, `SolidBody(..., apply=None)`. This may be used to sum the list of sub-blocks instead of stacking them together, `SolidBody(..., block=False, apply=None)`. This is useful for mixed formulations where both the deformation gradient and the displacement values are required. +- Add support for non-symmetric bilinear mixed forms in `IntegralForm`. ### Changed - The first Piola-Kirchhoff stress tensor is evaluated if `ViewSolid(stress_type=None)`. @@ -19,6 +20,7 @@ All notable changes to this project will be documented in this file. The format - Enhance `math.linsteps(..., values=0)` with default values except for the column `axis` if `axis` is not None. - Link all field-values to the values of the first field if no other field is given in `FieldContainer.link()`. - Change the default arguments from `block=True` to `block=None` in `SolidBody.assemble.vector(block=None)` and `SolidBody.assemble.matrix(block=None)` and define `block` on creation, `SolidBody(..., block=True)` instead. +- Integrate and assemble `None`-items in the `fun`-argument of `IntegralForm`. On integration, `None` is returned. `None` will be assembled to an emtpy sparse matrix. ### Fixed - Fix `Boundary(..., mode="and")` by ignoring any undefined axis. diff --git a/src/felupe/assembly/_cartesian.py b/src/felupe/assembly/_cartesian.py index 02e59bd8..19ea1ae2 100644 --- a/src/felupe/assembly/_cartesian.py +++ b/src/felupe/assembly/_cartesian.py @@ -104,7 +104,10 @@ class IntegralFormCartesian: """ def __init__(self, fun, v, dV, u=None, grad_v=False, grad_u=False): - self.fun = np.ascontiguousarray(fun) + self.fun = fun + if self.fun is not None: + self.fun = np.ascontiguousarray(self.fun) + self.dV = dV self.v = v @@ -137,20 +140,22 @@ def assemble(self, values=None, parallel=False, out=None): if values is None: values = self.integrate(parallel=parallel, out=out) - permute = np.append(len(values.shape) - 1, range(len(values.shape) - 1)).astype( - int - ) + if values is not None: + permute = np.append( + len(values.shape) - 1, range(len(values.shape) - 1) + ).astype(int) - # broadcast values of a uniform grid mesh - if values.size < self.indices[0].size: - new_shape = (*values.shape[:-1], self.v.region.mesh.ncells) - values = np.broadcast_to(values, new_shape) + # broadcast values of a uniform grid mesh + if values.size < self.indices[0].size: + new_shape = (*values.shape[:-1], self.v.region.mesh.ncells) + values = np.broadcast_to(values, new_shape) - res = sparsematrix( - (values.transpose(permute).ravel(), self.indices), shape=self.shape - ) + return sparsematrix( + (values.transpose(permute).ravel(), self.indices), shape=self.shape + ) - return res + else: + return sparsematrix(self.shape) def integrate(self, parallel=False, out=None): "Return evaluated (but not assembled) integrals." @@ -160,6 +165,9 @@ def integrate(self, parallel=False, out=None): dV = self.dV fun = self.fun + if fun is None: + return None + # plane strain # trim 3d vector-valued functions to the dimension of the field function_dimension = len(fun.shape) - 2 diff --git a/src/felupe/assembly/_integral.py b/src/felupe/assembly/_integral.py index a46aef93..9d9eb791 100644 --- a/src/felupe/assembly/_integral.py +++ b/src/felupe/assembly/_integral.py @@ -239,8 +239,8 @@ def __init__(self, fun, v, dV, u=None, grad_v=None, grad_u=None): f = IntForm(fun=fun, v=v, dV=self.dV, grad_v=grad_v) self.forms.append(f) - elif len(fun) == np.sum(1 + np.arange(self.nv)) and u is not None: - # BilinearForm + elif u is not None and len(fun) == np.sum(1 + np.arange(self.nv)): + # BilinearForm (symmetric, upper-triangle entries) self.mode = 2 self.i, self.j = np.triu_indices(self.nv) @@ -254,6 +254,23 @@ def __init__(self, fun, v, dV, u=None, grad_v=None, grad_u=None): grad_u=self.grad_u[j], ) self.forms.append(f) + + elif u is not None and len(fun) == self.nv * self.nu: + # BilinearForm (non-symmetric) + self.mode = 3 + self.i, self.j = [i.ravel() for i in np.indices((self.nv, self.nu))] + + for a, (i, j) in enumerate(zip(self.i, self.j)): + f = IntForm( + self.fun[a], + v=self.v[i], + dV=self.dV, + u=self.u[j], + grad_v=self.grad_v[i], + grad_u=self.grad_u[j], + ) + self.forms.append(f) + else: raise ValueError("Unknown input format.") @@ -266,11 +283,11 @@ def assemble(self, values=None, parallel=False, block=True, out=None): for val, form in zip(values, self.forms): res.append(form.assemble(val, parallel=parallel, out=out)) - if block and self.mode == 2: - K = np.zeros((self.nv, self.nv), dtype=object) + if block and (self.mode == 2 or self.mode == 3): + K = np.zeros((self.nv, self.nu), dtype=object) for a, (i, j) in enumerate(zip(self.i, self.j)): K[i, j] = res[a] - if i != j: + if self.mode == 2 and i != j: K[j, i] = res[a].T res = bmat(K).tocsr() diff --git a/src/felupe/constitution/_mixed.py b/src/felupe/constitution/_mixed.py index c4e24c1c..06894059 100644 --- a/src/felupe/constitution/_mixed.py +++ b/src/felupe/constitution/_mixed.py @@ -214,8 +214,8 @@ def hessian(self, x, out=None): d2WdFdF = self.material.hessian([F, statevars], **kwargs)[0] d2WdFdF += p * detF * (dya(iFT, iFT) - cdya_il(iFT, iFT)) d2WdFdp = detF * iFT - d2WdFdJ = np.zeros_like(F) - d2Wdpdp = np.zeros_like(p) + d2WdFdJ = None # np.zeros_like(F) + d2Wdpdp = None # np.zeros_like(p) d2WdpdJ = -np.ones_like(J) d2WdJdJ = self.d2UdJdJ(J, self.bulk) * np.ones_like(J) return [d2WdFdF, d2WdFdp, d2WdFdJ, d2Wdpdp, d2WdpdJ, d2WdJdJ] @@ -623,7 +623,7 @@ def _hessian_pp(self, F, p, J): ndarray p,p - part of hessian """ - return np.zeros_like(p) + return None # np.zeros_like(p) def _hessian_JJ(self, F, p, J): """Linearization w.r.t. volume ratio of variation of diff --git a/tests/test_form.py b/tests/test_form.py index a4551a7b..108fc31f 100644 --- a/tests/test_form.py +++ b/tests/test_form.py @@ -252,6 +252,20 @@ def test_mixed(): assert b.shape == (z, 1) + L = fem.IntegralForm([f[0], None, f[2]], v, r.dV) + x = L.integrate() + b = L.assemble(x) + + assert b.shape == (z, 1) + + a = fem.IntegralForm( + [A[0], A[1], A[2], A[1], A[3], A[4], A[2], A[4], A[5]], v, r.dV, v + ) + y = a.integrate() + K = a.assemble(y) + + assert K.shape == (z, z) + r, v, f, A = pre_axi_mixed() for parallel in [False, True]: diff --git a/tests/test_mechanics.py b/tests/test_mechanics.py index 6b5e0072..faeac28f 100644 --- a/tests/test_mechanics.py +++ b/tests/test_mechanics.py @@ -465,7 +465,10 @@ def test_solidbody_mixed(): A2 = b.evaluate.hessian() A2 = b.evaluate.hessian(u) for a1, a2 in zip(A1, A2): - assert np.allclose(a1, a2) + if a1 is None: + assert a2 is None + else: + assert np.allclose(a1, a2) F1 = b.results.kinematics F2 = b._extract(u) From 4aaf105caaa1b80c674319e9e6aa569f2de57d73 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Thu, 2 Jan 2025 23:16:22 +0100 Subject: [PATCH 19/47] Fix Ex20: Merge points of top-level mesh container also rename `gamma` to `gamma0` and add the parameter `kr` --- examples/ex20_third-medium-contact.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/examples/ex20_third-medium-contact.py b/examples/ex20_third-medium-contact.py index 352e9f72..1d4658ce 100644 --- a/examples/ex20_third-medium-contact.py +++ b/examples/ex20_third-medium-contact.py @@ -48,12 +48,13 @@ # down material parameters. G = 5 / 14 K = 5 / 3 -gamma = 2e-6 +kr = 1e-6 +gamma0 = 1e-6 -container = fem.MeshContainer([body, medium], merge=True) +container = fem.MeshContainer([body, medium], merge=True, decimals=6) regions = [fem.RegionQuad(m, hess=True) for m in container.meshes] fields = [fem.FieldPlaneStrain(r, dim=2).as_container() for r in regions] -umats = [fem.NeoHooke(mu=G, bulk=K), fem.NeoHooke(mu=G * gamma, bulk=K * gamma)] +umats = [fem.NeoHooke(mu=G, bulk=K), fem.NeoHooke(mu=G * gamma0, bulk=K * gamma0)] solids = [fem.SolidBody(umat, f) for umat, f in zip(umats, fields)] # %% @@ -97,15 +98,15 @@ def linearform(): regularization = fem.FormItem( bilinearform=bilinearform, linearform=linearform, - kwargs={"kr": K * L**2 * gamma}, + kwargs={"kr": kr * K * L**2}, ) # %% # The prescribed displacement is ramped up to the maximum value and released until zero. -move = fem.math.linsteps([0, 0.6, 1, 0.6, 0], num=15) +move = fem.math.linsteps([0, 1, 0], num=30) step = fem.Step( items=[*solids, regularization], - ramp={bounds["move"]: move * -0.5 * L}, + ramp={bounds["move"]: -0.4 * move * L}, boundaries=bounds, ) From 2ab6697114ee9cbe8319577e5d37939887a511b4 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Fri, 3 Jan 2025 22:52:08 +0100 Subject: [PATCH 20/47] Simplify init of base-class `element.Element` (#924) * Simplify `Element` and add `element.Element` to the top-level namespace * Update _lagrange.py * Don't use `Element.shape` in `"tetra10"` and `"triangle6"` --- CHANGELOG.md | 4 ++++ src/felupe/__init__.py | 2 ++ src/felupe/element/_base.py | 13 +++++------- src/felupe/element/_hexahedron.py | 5 ----- src/felupe/element/_lagrange.py | 6 +++--- src/felupe/element/_line.py | 1 - src/felupe/element/_quad.py | 5 ----- src/felupe/element/_tetra.py | 5 +---- src/felupe/element/_triangle.py | 5 +---- tests/test_element.py | 35 ------------------------------- 10 files changed, 16 insertions(+), 65 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 5e47b42d..878b32db 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,7 @@ All notable changes to this project will be documented in this file. The format - Add `BoundaryDict` as a subclassed dict with methods to `plot()`, `screenshot()` and `imshow()`. - Add a new argument to apply a callable on the assembled vector/matrix of a solid body, `SolidBody(..., apply=None)`. This may be used to sum the list of sub-blocks instead of stacking them together, `SolidBody(..., block=False, apply=None)`. This is useful for mixed formulations where both the deformation gradient and the displacement values are required. - Add support for non-symmetric bilinear mixed forms in `IntegralForm`. +- Add `element.Element` to the top-level package namespace. ### Changed - The first Piola-Kirchhoff stress tensor is evaluated if `ViewSolid(stress_type=None)`. @@ -25,6 +26,9 @@ All notable changes to this project will be documented in this file. The format ### Fixed - Fix `Boundary(..., mode="and")` by ignoring any undefined axis. +### Removed +- Remove the unused `shape`-argument in `element.Element(shape)`. Adopt the arbitrary-lagrange element to use its own `dim`-argument. This simplifies the creation of custom finite element formulations. + ## [9.1.0] - 2024-11-23 ### Added diff --git a/src/felupe/__init__.py b/src/felupe/__init__.py index 4d4fb7af..a5e16468 100644 --- a/src/felupe/__init__.py +++ b/src/felupe/__init__.py @@ -46,6 +46,7 @@ BiQuadraticQuad, ConstantHexahedron, ConstantQuad, + Element, Hexahedron, Line, Quad, @@ -173,6 +174,7 @@ "BiQuadraticQuad", "ConstantHexahedron", "ConstantQuad", + "Element", "Hexahedron", "Line", "Quad", diff --git a/src/felupe/element/_base.py b/src/felupe/element/_base.py index 676696e6..0fd3c1af 100644 --- a/src/felupe/element/_base.py +++ b/src/felupe/element/_base.py @@ -20,9 +20,6 @@ class Element: - def __init__(self, shape): - self.shape = shape - self.dim = self.shape[1] def view(self, point_data=None, cell_data=None, cell_type=None): """View the element with optional given dicts of point- and cell-data items. @@ -91,7 +88,7 @@ def plot( ) if add_point_labels: plotter.add_point_labels( - points=np.pad(self.points, ((0, 0), (0, 3 - self.shape[1]))), + points=np.pad(self.points, ((0, 0), (0, 3 - self.points.shape[1]))), labels=[f"{a}" for a in np.arange(len(self.points))], font_size=font_size, show_points=show_points, @@ -108,19 +105,19 @@ def plot( actor.SetNormalizedShaftLength((0.9, 0.9, 0.9)) actor.SetNormalizedTipLength((0.1, 0.1, 0.1)) - if self.shape[1] == 3: + if self.points.shape[1] == 3: actor.SetTotalLength([1.3, 1.3, 1.3]) - elif self.shape[1] == 2: + elif self.points.shape[1] == 2: actor.SetZAxisLabelText("") actor.SetTotalLength([1.3, 1.3, 0]) - elif self.shape[1] == 1: + elif self.points.shape[1] == 1: actor.SetYAxisLabelText("") actor.SetZAxisLabelText("") actor.SetTotalLength([1.3, 0, 0]) plotter.camera.zoom(0.7) - if self.shape[1] == 3: + if self.points.shape[1] == 3: plotter.camera.azimuth = -17 return plotter diff --git a/src/felupe/element/_hexahedron.py b/src/felupe/element/_hexahedron.py index d8264ca4..2721a044 100644 --- a/src/felupe/element/_hexahedron.py +++ b/src/felupe/element/_hexahedron.py @@ -69,7 +69,6 @@ def __init__(self): ) self.cells = np.arange(len(self.points)).reshape(1, -1) self.cell_type = "hexahedron" - super().__init__(shape=(1, 3)) def function(self, rst): "Return the shape functions at given coordinates (r, s, t)." @@ -140,7 +139,6 @@ def __init__(self): ) self.cells = np.arange(len(self.points)).reshape(1, -1) self.cell_type = "hexahedron" - super().__init__(shape=(8, 3)) def function(self, rst): "Return the shape functions at given coordinates (r, s, t)." @@ -263,7 +261,6 @@ def __init__(self): ) self.cells = np.arange(len(self.points)).reshape(1, -1) self.cell_type = "hexahedron20" - super().__init__(shape=(20, 3)) def function(self, rst): "Return the shape functions at given coordinates (r, s, t)." @@ -471,8 +468,6 @@ class TriQuadraticHexahedron(Element): """ def __init__(self): - super().__init__(shape=(27, 3)) - self.points = np.array( [ [-1, -1, -1], diff --git a/src/felupe/element/_lagrange.py b/src/felupe/element/_lagrange.py index fef8fe89..5d30036c 100644 --- a/src/felupe/element/_lagrange.py +++ b/src/felupe/element/_lagrange.py @@ -261,14 +261,14 @@ def __init__(self, order, dim, interval=(-1, 1), permute=True): self._nbasis = self._npoints self._interval = interval + self.dim = dim + self.permute = None if permute: self.permute = [None, lagrange_line, lagrange_quad, lagrange_hexahedron][ dim ](order) - super().__init__(shape=(self._npoints, dim)) - # init curve-parameter matrix n = self._nshape self._AT = np.linalg.inv( @@ -278,7 +278,7 @@ def __init__(self, order, dim, interval=(-1, 1), permute=True): # indices for outer product in einstein notation # idx = ["a", "b", "c", ...][:dim] # subscripts = "a,b,c -> abc" - self._idx = [letter for letter in alphabet][: self.dim] + self._idx = [letter for letter in alphabet][:dim] self._subscripts = ",".join(self._idx) + "->" + "".join(self._idx) # init points diff --git a/src/felupe/element/_line.py b/src/felupe/element/_line.py index 12c44aaf..2b9b3d75 100644 --- a/src/felupe/element/_line.py +++ b/src/felupe/element/_line.py @@ -45,7 +45,6 @@ class Line(Element): """ def __init__(self): - super().__init__(shape=(2, 1)) self.points = np.array([-1, 1], dtype=float).reshape(-1, 1) self.cells = np.arange(len(self.points)).reshape(1, -1) self.cell_type = "line" diff --git a/src/felupe/element/_quad.py b/src/felupe/element/_quad.py index 3491181f..65ecc52e 100644 --- a/src/felupe/element/_quad.py +++ b/src/felupe/element/_quad.py @@ -52,7 +52,6 @@ class ConstantQuad(Element): """ def __init__(self): - super().__init__(shape=(1, 2)) self.points = np.array([[-1, -1], [1, -1], [1, 1], [-1, 1]], dtype=float) self.cells = np.arange(len(self.points)).reshape(1, -1) self.cell_type = "quad" @@ -106,7 +105,6 @@ class Quad(Element): """ def __init__(self): - super().__init__(shape=(4, 2)) self.points = np.array([[-1, -1], [1, -1], [1, 1], [-1, 1]], dtype=float) self.cells = np.arange(len(self.points)).reshape(1, -1) self.cell_type = "quad" @@ -188,7 +186,6 @@ class QuadraticQuad(Element): """ def __init__(self): - super().__init__(shape=(8, 2)) self.points = np.array( [ [-1, -1], @@ -289,8 +286,6 @@ class BiQuadraticQuad(Element): """ def __init__(self): - super().__init__(shape=(9, 2)) - self._lagrange = ArbitraryOrderLagrange(order=2, dim=2, permute=False) self._vertices = np.array([0, 2, 8, 6]) diff --git a/src/felupe/element/_tetra.py b/src/felupe/element/_tetra.py index 58e315c1..2e30974c 100644 --- a/src/felupe/element/_tetra.py +++ b/src/felupe/element/_tetra.py @@ -57,7 +57,6 @@ class Tetra(Element): """ def __init__(self): - super().__init__(shape=(4, 3)) self.points = np.array( [[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1]], dtype=float ) @@ -117,7 +116,6 @@ class TetraMINI(Element): """ def __init__(self, bubble_multiplier=1.0): - super().__init__(shape=(5, 3)) self.points = np.array( [[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1], [1 / 3, 1 / 3, 1 / 3]], dtype=float, @@ -237,8 +235,7 @@ class QuadraticTetra(Element): """ def __init__(self): - super().__init__(shape=(10, 3)) - self.points = np.zeros(self.shape) + self.points = np.zeros((10, 3)) self.points[:4] = np.array( [[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1]], dtype=float ) diff --git a/src/felupe/element/_triangle.py b/src/felupe/element/_triangle.py index deeaf797..8285c8ea 100644 --- a/src/felupe/element/_triangle.py +++ b/src/felupe/element/_triangle.py @@ -58,7 +58,6 @@ class Triangle(Element): """ def __init__(self): - super().__init__(shape=(3, 2)) self.points = np.array([[0, 0], [1, 0], [0, 1]], dtype=float) self.cells = np.arange(len(self.points)).reshape(1, -1) self.cell_type = "triangle" @@ -115,7 +114,6 @@ class TriangleMINI(Element): """ def __init__(self, bubble_multiplier=1.0): - super().__init__(shape=(4, 2)) self.points = np.array([[0, 0], [1, 0], [0, 1], [1 / 3, 1 / 3]], dtype=float) self.cells = np.arange(len(self.points) - 1).reshape(1, -1) self.cell_type = "triangle" @@ -198,8 +196,7 @@ class QuadraticTriangle(Element): """ def __init__(self): - super().__init__(shape=(6, 2)) - self.points = np.zeros(self.shape) + self.points = np.zeros((6, 2)) self.points[:3] = np.array([[0, 0], [1, 0], [0, 1]], dtype=float) self.points[3] = np.mean(self.points[[0, 1]], axis=0) self.points[4] = np.mean(self.points[[1, 2]], axis=0) diff --git a/tests/test_element.py b/tests/test_element.py index f90aaac8..a8c78571 100644 --- a/tests/test_element.py +++ b/tests/test_element.py @@ -43,8 +43,6 @@ def test_line2(): assert np.all(dhdr[0] == -0.5) assert np.all(d2hdrdr == 0.0) - assert line2.shape == dhdr.shape - line2.plot(off_screen=True) @@ -59,8 +57,6 @@ def test_line_lagrange(): assert np.isclose(h[0], 1) assert np.isclose(dhdr[0, 0], -5.70833333) - assert line6.shape == dhdr.shape - line6.plot(off_screen=True) @@ -77,8 +73,6 @@ def test_quad0(): assert np.all(dhdr[0] == 0) assert np.all(d2hdrdr == 0.0) - assert quad0.shape == dhdr.shape - def test_quad4(): quad4 = fem.element.Quad() @@ -93,8 +87,6 @@ def test_quad4(): assert np.all(dhdr[0] == -0.5) assert d2hdrdr.shape == (4, 2, 2) - assert quad4.shape == dhdr.shape - quad4.plot(off_screen=True) @@ -111,8 +103,6 @@ def test_quad8(): assert np.all(dhdr[0] == -1.5) assert d2hdrdr.shape == (8, 2, 2) - assert quad8.shape == dhdr.shape - def test_quad9(): quad9 = fem.element.BiQuadraticQuad() @@ -125,8 +115,6 @@ def test_quad9(): assert h[0] == 1 assert np.all(dhdr[0] == -1.5) - assert quad9.shape == dhdr.shape - def test_hex0(): hex0 = fem.element.ConstantHexahedron() @@ -141,8 +129,6 @@ def test_hex0(): assert np.all(dhdr[0] == 0) assert np.all(d2hdrdr == 0) - assert hex0.shape == dhdr.shape - def test_hex8(): hex8 = fem.element.Hexahedron() @@ -157,8 +143,6 @@ def test_hex8(): assert np.all(dhdr[0] == -0.5) assert d2hdrdr.shape == (8, 3, 3) - assert hex8.shape == dhdr.shape - def test_hex20(): hex20 = fem.element.QuadraticHexahedron() @@ -171,8 +155,6 @@ def test_hex20(): assert h[0] == 1 assert np.all(dhdr[0] == -1.5) - assert hex20.shape == dhdr.shape - hex20.plot(off_screen=True) @@ -187,8 +169,6 @@ def test_hex27(): assert h[0] == 1 assert np.all(dhdr[0] == -1.5) - assert hex27.shape == dhdr.shape - hex27.plot(off_screen=True) @@ -205,9 +185,6 @@ def test_tri3(): assert np.all(dhdr[0] == -1) assert np.all(d2hdrdr == 0) - assert tri3.shape == dhdr.shape - - def test_tri6(): tri6 = fem.element.QuadraticTriangle() @@ -219,8 +196,6 @@ def test_tri6(): assert h[0] == 1 assert np.all(dhdr[0] == -3) - assert tri6.shape == dhdr.shape - def test_tri_mini(): trim = fem.element.TriangleMINI() @@ -236,8 +211,6 @@ def test_tri_mini(): assert np.all(dhdr[0] == -1) assert d2hdrdr.shape == (4, 2, 2) - assert trim.shape == dhdr.shape - def test_tet4(): tet4 = fem.element.Tetra() @@ -252,8 +225,6 @@ def test_tet4(): assert np.all(dhdr[0] == -1) assert np.all(d2hdrdr == 0) - assert tet4.shape == dhdr.shape - def test_tet10(): tet10 = fem.element.QuadraticTetra() @@ -266,8 +237,6 @@ def test_tet10(): assert h[0] == 1 assert np.all(dhdr[0] == -3) - assert tet10.shape == dhdr.shape - def test_tet_mini(): tetm = fem.element.TetraMINI() @@ -283,8 +252,6 @@ def test_tet_mini(): assert np.all(dhdr[0] == -1) assert d2hdrdr.shape == (5, 3, 3) - assert tetm.shape == dhdr.shape - def test_aol(): aol32 = fem.element.ArbitraryOrderLagrange(order=3, dim=2) @@ -296,7 +263,6 @@ def test_aol(): dhdr = aol32.gradient(r) assert h[0] == 1 - assert aol32.shape == dhdr.shape r = [-1, -1, -1] @@ -304,7 +270,6 @@ def test_aol(): dhdr = aol23.gradient(r) assert h[0] == 1 - assert aol23.shape == dhdr.shape if __name__ == "__main__": From 1132065b463ce5ff9b6b98d5843e08225bb56ddb Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Fri, 3 Jan 2025 23:03:25 +0100 Subject: [PATCH 21/47] Add a docstring for `Element` --- docs/felupe/element.rst | 13 ++++++++++-- src/felupe/element/_base.py | 41 +++++++++++++++++++++++++++++++++++++ 2 files changed, 52 insertions(+), 2 deletions(-) diff --git a/docs/felupe/element.rst b/docs/felupe/element.rst index d65c2528..ec91e7c2 100644 --- a/docs/felupe/element.rst +++ b/docs/felupe/element.rst @@ -5,10 +5,14 @@ Element This module provides classes for the finite element formulations. -**Linear Elements** - .. currentmodule:: felupe +.. autosummary:: + + Element + +**Linear Elements** + .. autosummary:: Line @@ -50,6 +54,11 @@ This module provides classes for the finite element formulations. **Detailed API Reference** +.. autoclass:: felupe.Element + :members: + :undoc-members: + :show-inheritance: + .. autoclass:: felupe.Line :members: :undoc-members: diff --git a/src/felupe/element/_base.py b/src/felupe/element/_base.py index 0fd3c1af..877e5d8c 100644 --- a/src/felupe/element/_base.py +++ b/src/felupe/element/_base.py @@ -20,6 +20,47 @@ class Element: + """ + Base-class for a finite element which provides methods for plotting. + + Examples + -------- + This example shows how to implement a hexahedron element. + + .. pyvista-plot:: + :force_static: + + >>> import felupe as fem + >>> 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 function(self, rst): + ... r, s, t = rst + ... a, b, c = self.points.T + ... ar, bs, ct = 1 + a * r, 1 + b * s, 1 + c * t + ... return (ar * bs * ct) / 8 + ... + ... def gradient(self, rst): + ... r, s, t = rst + ... a, b, c = self.points.T + ... ar, bs, ct = 1 + a * r, 1 + b * s, 1 + c * t + ... return np.stack([a * bs * ct, ar * b * ct, ar * bs * c], axis=1) + >>> + >>> mesh = fem.Cube(n=6) + >>> element = Hexahedron() + >>> quadrature = fem.GaussLegendre(order=1, dim=3) + >>> region = fem.Region(mesh, element, quadrature) + """ def view(self, point_data=None, cell_data=None, cell_type=None): """View the element with optional given dicts of point- and cell-data items. From ea697b874832cc433ab63b94619cfd1a4e21cad9 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Wed, 22 Jan 2025 14:31:09 +0100 Subject: [PATCH 22/47] Fix `hello_word(planestrain=True)` (#928) --- CHANGELOG.md | 1 + src/felupe/tools/_hello_world.py | 1 + 2 files changed, 2 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 878b32db..ef10f40d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -25,6 +25,7 @@ All notable changes to this project will be documented in this file. The format ### Fixed - Fix `Boundary(..., mode="and")` by ignoring any undefined axis. +- Fix `tools.hello_world(planestrain=True)` with the correct region `RegionQuad` for the plane-strain template. ### Removed - Remove the unused `shape`-argument in `element.Element(shape)`. Adopt the arbitrary-lagrange element to use its own `dim`-argument. This simplifies the creation of custom finite element formulations. diff --git a/src/felupe/tools/_hello_world.py b/src/felupe/tools/_hello_world.py index 5986a345..1f070882 100644 --- a/src/felupe/tools/_hello_world.py +++ b/src/felupe/tools/_hello_world.py @@ -74,6 +74,7 @@ def hello_world( elif planestrain: mesh = "Rectangle" + region += "Quad" field += "PlaneStrain" dim = 2 From c5d0e2a083d82ae19207bf64af5b22ed788538cb Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Thu, 23 Jan 2025 09:24:32 +0100 Subject: [PATCH 23/47] Update .readthedocs.yaml --- .readthedocs.yaml | 22 ++++++++++++---------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/.readthedocs.yaml b/.readthedocs.yaml index 3fd675c2..8731da19 100644 --- a/.readthedocs.yaml +++ b/.readthedocs.yaml @@ -3,21 +3,23 @@ version: 2 build: os: ubuntu-lts-latest tools: - python: "3.12" + python: "3.13" apt_packages: - libgl1-mesa-dev - xvfb - pandoc + jobs: + create_environment: + - asdf plugin add uv + - asdf install uv latest + - asdf global uv latest + - uv venv + install: + - uv pip install .[all,docs,examples] + build: + html: + - uv run sphinx-build -T -b html docs $READTHEDOCS_OUTPUT/html sphinx: configuration: docs/conf.py fail_on_warning: false - -python: - install: - - method: pip - path: . - extra_requirements: - - all - - docs - - examples From 586cab1a49668c66a2a69025d4747328036f50da Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Thu, 23 Jan 2025 09:37:59 +0100 Subject: [PATCH 24/47] Revert "Update .readthedocs.yaml" This reverts commit c5d0e2a083d82ae19207bf64af5b22ed788538cb. --- .readthedocs.yaml | 22 ++++++++++------------ 1 file changed, 10 insertions(+), 12 deletions(-) diff --git a/.readthedocs.yaml b/.readthedocs.yaml index 8731da19..3fd675c2 100644 --- a/.readthedocs.yaml +++ b/.readthedocs.yaml @@ -3,23 +3,21 @@ version: 2 build: os: ubuntu-lts-latest tools: - python: "3.13" + python: "3.12" apt_packages: - libgl1-mesa-dev - xvfb - pandoc - jobs: - create_environment: - - asdf plugin add uv - - asdf install uv latest - - asdf global uv latest - - uv venv - install: - - uv pip install .[all,docs,examples] - build: - html: - - uv run sphinx-build -T -b html docs $READTHEDOCS_OUTPUT/html sphinx: configuration: docs/conf.py fail_on_warning: false + +python: + install: + - method: pip + path: . + extra_requirements: + - all + - docs + - examples From a6c23cfc33b89a6f24fb016b257b2cba5b3203fa Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Thu, 23 Jan 2025 10:02:37 +0100 Subject: [PATCH 25/47] Update README.md --- README.md | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/README.md b/README.md index 4031d26f..2fe1a4bc 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) [![PyPI version shields.io](https://img.shields.io/pypi/v/felupe.svg)](https://pypi.python.org/pypi/felupe/) [![Documentation Status](https://readthedocs.org/projects/felupe/badge/?version=latest)](https://felupe.readthedocs.io/en/latest/?badge=latest) [![License: GPL v3](https://img.shields.io/badge/License-GPLv3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0) [![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](https://img.shields.io/badge/%F0%9F%94%8D-FElupe-white)](https://felupe.readthedocs.io) [![PyPI version shields.io](https://img.shields.io/pypi/v/felupe.svg)](https://pypi.python.org/pypi/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) [![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. @@ -50,7 +50,7 @@ Install Python, fire up 🔥 a terminal and run 🏃 pip install felupe[all] ``` -The [documentation](https://felupe.readthedocs.io/) covers more details, like required and optional dependencies and how to install the latest development version. +The [documentation](https://felupe.readthedocs.io/en/stable/) covers more details, like required and optional dependencies and how to install the latest development version. # Getting Started This tutorial covers the essential high-level parts of creating and solving problems with FElupe. As an introductory example 👨‍🏫, a quarter model of a solid cube with hyperelastic material behaviour is subjected to a uniaxial elongation applied at a clamped end-face. @@ -63,7 +63,7 @@ An isotropic pseudo-elastic Ogden-Roxburgh Mullins-softening model formulation i A step generates the consecutive substep-movements of a given boundary condition. The step is further added to a list of steps of a job 👩‍💻 (here, a characteristic-curve 📈 job is used). During evaluation ⏳, each substep of each step is solved by an iterative Newton-Rhapson procedure ⚖️. The solution is exported after each completed substep as a time-series ⌚ XDMF file. Finally, the result of the last completed substep is plotted. -For more details beside this high-level code snippet, please have a look at the 📝 [documentation](https://felupe.readthedocs.io). +For more details beside this high-level code snippet, please have a look at the 📝 [documentation](https://felupe.readthedocs.io/en/stable/). ```python import felupe as fem @@ -96,7 +96,7 @@ solid.plot("Principal Values of Cauchy Stress").show()

# Documentation -The documentation is located [here](https://felupe.readthedocs.io/en/latest/?badge=latest). +The documentation is located [here](https://felupe.readthedocs.io/en/stable/). # Extension Packages The capabilities of FElupe may be enhanced with extension packages created by the community. From ee0f15386f868025593651ae122cb002f30f49ad Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Thu, 23 Jan 2025 15:40:54 +0100 Subject: [PATCH 26/47] Add more badges repostatus, supported Python versions --- README.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index 2fe1a4bc..1ed16c5e 100644 --- a/README.md +++ b/README.md @@ -3,7 +3,8 @@

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) [![PyPI version shields.io](https://img.shields.io/pypi/v/felupe.svg)](https://pypi.python.org/pypi/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) [![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](https://img.shields.io/badge/%F0%9F%94%8D-FElupe-white)](https://felupe.readthedocs.io) [![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) + [![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 fb1bb162ee5013412ec2bbb9d1e81f4e877fd218 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Thu, 6 Feb 2025 22:29:13 +0100 Subject: [PATCH 27/47] Update README.md --- README.md | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/README.md b/README.md index 1ed16c5e..1ec2d59f 100644 --- a/README.md +++ b/README.md @@ -120,6 +120,8 @@ tox # Scientific Publications A list of articles in which [![FElupe](https://img.shields.io/badge/%F0%9F%94%8D-FElupe-white)](https://felupe.readthedocs.io) is involved. If you use FElupe in your scientific work, please star this repository, cite it [![DOI](https://zenodo.org/badge/360657894.svg)](https://zenodo.org/badge/latestdoi/360657894) and add your publication to this list. +
Expand the list... + - A. Dutzler, C. Buzzi, and M. Leitner, "Nondimensional translational characteristics of elastomer components", [Journal of Applied Engineering Design and Simulation](https://jaeds.uitm.edu.my/index.php/jaeds), vol. 1, no. 1. UiTM Press, Universiti Teknologi MARA, Sep. 21, 2021. doi: [10.24191/jaeds.v1i1.20](https://doi.org/10.24191/jaeds.v1i1.20). [![medium-story](https://img.shields.io/badge/medium-story-white)](https://medium.com/@adtzlr/nonlinear-force-displacement-curves-of-rubber-metal-parts-ab7c48448e96) - C. Buzzi, A. Dutzler, T. Faethe, J. Lassacher, M. Leitner, and F.-J. Weber, "Development of a tool for estimating @@ -128,6 +130,8 @@ Bogies and Running Gears](https://transportation.bme.hu/2022/11/30/bogie22/), 20 - J. Torggler, A. Dutzler, B. Oberdorfer, T. Faethe, H. Müller, C. Buzzi, and M. Leitner, "Investigating Damage Mechanisms in Cord-Rubber Composite Air Spring Bellows of Rail Vehicles and Representative Specimen Design", [Applied Composite Materials](https://www.springer.com/journal/10443). Springer Science and Business Media LLC, Aug. 22, 2023. doi: [10.1007/s10443-023-10157-1](https://link.springer.com/article/10.1007/s10443-023-10157-1). Simulation-related Python scripts are available on GitHub at [adtzlr/fiberreinforcedrubber](https://github.com/adtzlr/fiberreinforcedrubber). Open In Colab +
+ # Changelog All notable changes to this project will be documented in [this file](CHANGELOG.md). The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). From 98cbaf061ad42a7318d071faf8ed33fc953b6ca9 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Thu, 6 Feb 2025 22:32:27 +0100 Subject: [PATCH 28/47] Update README.md --- README.md | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index 1ec2d59f..9a4ed0a9 100644 --- a/README.md +++ b/README.md @@ -54,6 +54,9 @@ pip install felupe[all] The [documentation](https://felupe.readthedocs.io/en/stable/) covers more details, like required and optional dependencies and how to install the latest development version. # Getting Started + +
Expand the description... + This tutorial covers the essential high-level parts of creating and solving problems with FElupe. As an introductory example 👨‍🏫, a quarter model of a solid cube with hyperelastic material behaviour is subjected to a uniaxial elongation applied at a clamped end-face. First, let’s import FElupe and create a meshed cube out of hexahedron cells with a given number of points per axis. A numeric region, pre-defined for hexahedrons, is created on the mesh. A vector-valued displacement field is initiated on the region. Next, a field container is created on top of this field. @@ -66,6 +69,8 @@ A step generates the consecutive substep-movements of a given boundary condition For more details beside this high-level code snippet, please have a look at the 📝 [documentation](https://felupe.readthedocs.io/en/stable/). +
+ ```python import felupe as fem @@ -118,7 +123,7 @@ tox ``` # Scientific Publications -A list of articles in which [![FElupe](https://img.shields.io/badge/%F0%9F%94%8D-FElupe-white)](https://felupe.readthedocs.io) is involved. If you use FElupe in your scientific work, please star this repository, cite it [![DOI](https://zenodo.org/badge/360657894.svg)](https://zenodo.org/badge/latestdoi/360657894) and add your publication to this list. +A list of articles in which [![FElupe](https://img.shields.io/badge/%F0%9F%94%8D-FElupe-white)](https://felupe.readthedocs.io/stable) is involved. If you use FElupe in your scientific work, please star this repository, cite it [![DOI](https://zenodo.org/badge/360657894.svg)](https://zenodo.org/badge/latestdoi/360657894) and add your publication to this list.
Expand the list... From 7f40c66f336bd3133e516743a65d57385eccfa1a Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Fri, 7 Feb 2025 23:46:26 +0100 Subject: [PATCH 29/47] Fix `RegionQuadraticQuadBoundary` (#933) * Fix `RegionQuadraticQuadBoundary` fix the cells-array for midpoints of edges * Capture some warnings in the tests * Optionally close figures in tests --- CHANGELOG.md | 1 + src/felupe/region/_boundary.py | 6 +-- tests/test_constitution.py | 96 +++++++++++++++++++++++++++------- tests/test_job.py | 4 +- tests/test_region.py | 10 ++++ 5 files changed, 93 insertions(+), 24 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index ef10f40d..3f682575 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -26,6 +26,7 @@ All notable changes to this project will be documented in this file. The format ### Fixed - Fix `Boundary(..., mode="and")` by ignoring any undefined axis. - Fix `tools.hello_world(planestrain=True)` with the correct region `RegionQuad` for the plane-strain template. +- Fix the cells-array in `RegionQuadraticQuadBoundary` for midpoints of edges. ### Removed - Remove the unused `shape`-argument in `element.Element(shape)`. Adopt the arbitrary-lagrange element to use its own `dim`-argument. This simplifies the creation of custom finite element formulations. diff --git a/src/felupe/region/_boundary.py b/src/felupe/region/_boundary.py index 15076b4c..59971aef 100644 --- a/src/felupe/region/_boundary.py +++ b/src/felupe/region/_boundary.py @@ -71,9 +71,9 @@ def boundary_cells_quad8(mesh): # complementary midpoints of edges for the creation of "boundary" quadratic quads # (rotated quads with 1st edge as n-th edge of one original quad) - n = [5, 6, 7, 4] - p = [4, 5, 6, 7] - q = [6, 7, 4, 5] + n = [4, 6, 5, 7] + p = [5, 7, 6, 4] + q = [6, 4, 7, 5] cells = np.dstack( ( diff --git a/tests/test_constitution.py b/tests/test_constitution.py index 05891b88..5bdf470e 100644 --- a/tests/test_constitution.py +++ b/tests/test_constitution.py @@ -302,7 +302,7 @@ def elasticity(x, mu, lmbda): dsde = linear_elastic.hessian([F, None]) -def test_umat_hyperelastic_statevars(): +def test_umat_hyperelastic_statevars(close_figs=True): r, x = pre(sym=False, add_identity=True) F = x[0] @@ -329,8 +329,11 @@ def test_umat_hyperelastic_statevars(): ux = fem.math.linsteps([1, 2, 1], num=10) ax = umat.plot(ux=ux, bx=None, ps=None, incompressible=True) + if close_figs: + plt.close(ax.get_figure()) + -def test_umat_hyperelastic(): +def test_umat_hyperelastic(close_figs=True): r, x = pre(sym=False, add_identity=True) F = x[0] @@ -433,18 +436,34 @@ def neo_hooke(C, mu=1): for incompressible in [False, True]: ax = umat.plot(incompressible=incompressible) + + if close_figs: + plt.close(ax.get_figure()) + ax = umat.screenshot(incompressible=incompressible) + if close_figs: + plt.close(ax.get_figure()) + umat = fem.Hyperelastic(fem.neo_hooke, mu=np.nan) with pytest.raises(ValueError): - umat.plot(bx=None, ps=None) + ax = umat.plot(bx=None, ps=None) + + if close_figs: + plt.close(ax.get_figure()) with pytest.raises(ValueError): - umat.plot(ux=None, bx=None) + ax = umat.plot(ux=None, bx=None) + + if close_figs: + plt.close(ax.get_figure()) with pytest.raises(ValueError): - umat.plot(ux=None, ps=None) + ax = umat.plot(ux=None, ps=None) + + if close_figs: + plt.close(ax.get_figure()) def test_umat_hyperelastic2(): @@ -476,7 +495,7 @@ def neo_hooke(F, mu=1): assert np.allclose(dsde, dsde2) -def test_umat_viscoelastic(): +def test_umat_viscoelastic(close_figs=True): r, x = pre(sym=False, add_identity=True, add_random=True) F = x[0] @@ -503,13 +522,16 @@ def viscoelastic(C, Cin, mu, eta, dtime): umat = fem.Hyperelastic( fem.constitution.finite_strain_viscoelastic, nstatevars=6, **kwargs ) - umat.plot( + ax = umat.plot( ux=fem.math.linsteps([1, 1.5, 1, 2, 1, 2.5, 1], num=15), ps=None, bx=None, incompressible=True, ) + if close_figs: + plt.close(ax.get_figure()) + s2, statevars_new = umat.gradient([F, statevars]) dsde2 = umat.hessian([F, statevars]) @@ -598,7 +620,7 @@ def test_elpliso(): dsde = umat.hessian([F, statevars]) -def test_composite(): +def test_composite(close_figs=True): r, x = pre(sym=False, add_identity=True) F = x[0] @@ -608,11 +630,14 @@ def test_composite(): ax = umat.plot() + if close_figs: + plt.close(ax.get_figure()) + dWdF, statevars_new = umat.gradient([F, None]) (d2WdFdF,) = umat.hessian([F, None]) -def test_optimize(): +def test_optimize(close_figs=True): stretches, stresses = ( np.array( [ @@ -646,11 +671,24 @@ def test_optimize(): * np.array([1.0, 0.0980665]) ).T + for model in [ + fem.extended_tube, + ]: + umat = fem.Hyperelastic(model) + with pytest.warns(): + umat_new, res = umat.optimize(ux=[stretches, stresses], incompressible=True) + + ux = np.linspace(stretches.min(), stretches.max(), num=200) + ax = umat_new.plot(incompressible=True, ux=ux, bx=None, ps=None) + ax.plot(stretches, stresses, "C0x") + + if close_figs: + plt.close(ax.get_figure()) + for model in [ fem.neo_hooke, fem.ogden, fem.third_order_deformation, - fem.extended_tube, fem.alexander, ]: umat = fem.Hyperelastic(model) @@ -660,6 +698,9 @@ def test_optimize(): ax = umat_new.plot(incompressible=True, ux=ux, bx=None, ps=None) ax.plot(stretches, stresses, "C0x") + if close_figs: + plt.close(ax.get_figure()) + for model in [ fem.neo_hooke, fem.miehe_goektepe_lulei, @@ -673,6 +714,9 @@ def test_optimize(): ax = umat_new.plot(incompressible=True, ux=ux, bx=None, ps=None) ax.plot(stretches, stresses, "C0x") + if close_figs: + plt.close(ax.get_figure()) + def test_lagrange(): r, x = pre(sym=False, add_identity=True) @@ -766,22 +810,34 @@ def test_laplace(): assert A[0].shape == (3, 3, 3, 3, 1, 1) -def test_plot_negative_stretches(): +def test_plot_negative_stretches(close_figs=True): stretches = np.linspace(-0.5, 1, 16) umat = fem.NeoHooke(mu=1.0, bulk=2.0) for incompressible in [False, True]: with pytest.raises(ValueError): - umat.plot(ux=stretches, incompressible=incompressible) + ax = umat.plot(ux=stretches, incompressible=incompressible) + + if close_figs: + plt.close(ax.get_figure()) with pytest.raises(ValueError): - umat.plot(bx=stretches, incompressible=incompressible) + ax = umat.plot(bx=stretches, incompressible=incompressible) + + if close_figs: + plt.close(ax.get_figure()) with pytest.raises(ValueError): - umat.plot(ps=stretches, incompressible=incompressible) + ax = umat.plot(ps=stretches, incompressible=incompressible) + + if close_figs: + plt.close(ax.get_figure()) if __name__ == "__main__": + + close_figs = True + test_nh() test_linear() test_linear_orthotropic() @@ -789,17 +845,17 @@ def test_plot_negative_stretches(): test_linear_planestrain() test_kinematics() test_umat() - test_umat_hyperelastic() + test_umat_hyperelastic(close_figs=close_figs) test_umat_hyperelastic2() - test_umat_hyperelastic_statevars() - test_umat_viscoelastic() + test_umat_hyperelastic_statevars(close_figs=close_figs) + test_umat_viscoelastic(close_figs=close_figs) test_umat_viscoelastic2() test_umat_strain() test_umat_strain_plasticity() test_elpliso() - test_composite() - test_optimize() + test_composite(close_figs=close_figs) + test_optimize(close_figs=close_figs) test_lagrange() test_lagrange_statevars() test_laplace() - test_plot_negative_stretches() + test_plot_negative_stretches(close_figs=True) diff --git a/tests/test_job.py b/tests/test_job.py index 854be430..0ce0a992 100644 --- a/tests/test_job.py +++ b/tests/test_job.py @@ -55,7 +55,9 @@ def test_job(): def test_job_xdmf(): - field, step = pre(SolidBodyForce=fem.SolidBodyGravity) + with pytest.warns(): # solidbodygravity is deprecated + field, step = pre(SolidBodyForce=fem.SolidBodyGravity) + job = fem.Job(steps=[step]) job.evaluate() diff --git a/tests/test_region.py b/tests/test_region.py index d09f6628..93b1ccb1 100644 --- a/tests/test_region.py +++ b/tests/test_region.py @@ -166,5 +166,15 @@ def test_region(): assert not np.any(r.dV <= 0) +def test_region_negative_volumes_cells(): + + mesh = fem.Rectangle() + mesh.cells = mesh.cells[:, ::-1] + + with pytest.warns(): # negative volumes of cells + region = fem.RegionQuad(mesh) + + if __name__ == "__main__": test_region() + test_region_negative_volumes_cells() From 0fcab9cc52f7890fdca664e703e369f27c0c30e6 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Sat, 8 Feb 2025 08:28:45 +0100 Subject: [PATCH 30/47] Don't create a figure if error in `umat.plot()` (#934) --- src/felupe/constitution/_view.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/felupe/constitution/_view.py b/src/felupe/constitution/_view.py index 5a34932e..2ac94cba 100644 --- a/src/felupe/constitution/_view.py +++ b/src/felupe/constitution/_view.py @@ -72,11 +72,12 @@ 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() - data = self.evaluate() for stretch, force, label in data: ax.plot(stretch, force, label=label, **kwargs) From 0b3a63713d26012a276b2e991cd28a4780b775ae Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Mon, 10 Feb 2025 16:45:01 +0100 Subject: [PATCH 31/47] Add docstrings to `solve.solve()` and `solve.partition()` --- src/felupe/mechanics/_solidbody.py | 1 + src/felupe/solve/_solve.py | 103 +++++++++++++++++++++++++++-- 2 files changed, 100 insertions(+), 4 deletions(-) diff --git a/src/felupe/mechanics/_solidbody.py b/src/felupe/mechanics/_solidbody.py index 6824ed75..4fd0419f 100644 --- a/src/felupe/mechanics/_solidbody.py +++ b/src/felupe/mechanics/_solidbody.py @@ -221,6 +221,7 @@ class SolidBody(Solid): Examples -------- .. pyvista-plot:: + :force_static: >>> import felupe as fem >>> diff --git a/src/felupe/solve/_solve.py b/src/felupe/solve/_solve.py index 868eb788..be2b0b0a 100644 --- a/src/felupe/solve/_solve.py +++ b/src/felupe/solve/_solve.py @@ -25,7 +25,67 @@ def partition(v, K, dof1, dof0, r=None): """Perform partitioning of field values (unknowns), (stiffness) matrix and (residuals) vector with given lists of active (dof1) and - prescribed degrees of freedom (dof0).""" + prescribed degrees of freedom (dof0). + + Parameters + ---------- + v : FieldContainer + A field container with the fields. It is used to extract and concatenate the + 1d-array of unknows. + K : scipy.sparse.spmatrix + A two-dimensional sparse (stiffness) matrix. + dof1 : list of int + List of active degrees of freedom (zero-indexed). + dof0 : list of int + List of prescribed degrees of freedom (zero-indexed). + r : 1d-array or None, optional + The 1d-vector of residuals. If None, the residuals are set to zero. Default is + None. + + Returns + ------- + u : 1d-array + The concatenated full 1d-array of unknowns, extracted from the field container. + u0 : 1d-array + The prescribed unknowns. + K11 : scipy.sparse.spmatrix + The active-active block of the two-dimensional sparse (stiffness) matrix. + K10 : scipy.sparse.spmatrix + The active-prescribed block of the two-dimensional sparse (stiffness) matrix. + dof1 : list of int + List of active degrees of freedom (zero-indexed). + dof0 : list of int + List of prescribed degrees of freedom (zero-indexed). + r1 : 1d-array + The 1d-vector of active residuals. + + Examples + -------- + .. pyvista-plot:: + :force_static: + + >>> import felupe as fem + >>> + >>> mesh = fem.Rectangle(n=3) + >>> region = fem.RegionQuad(mesh) + >>> field = fem.FieldPlaneStrain(region, dim=2).as_container() + >>> boundaries, loadcase = fem.dof.uniaxial(field, clamped=True) + >>> umat = fem.NeoHooke(mu=1.0, bulk=2.0) + >>> solid = fem.SolidBody(umat, field) + >>> + >>> K = solid.assemble.matrix() + >>> r = solid.assemble.vector() + >>> + >>> dof0 = loadcase["dof0"] + >>> dof1 = loadcase["dof1"] + >>> ext0 = loadcase["ext0"] + >>> + >>> system = fem.solve.partition(field, K, dof1, dof0, r) + + See Also + -------- + felupe.solve.solve : Linear solution of equation system. + """ # extract values u = values(v) @@ -48,11 +108,46 @@ def partition(v, K, dof1, dof0, r=None): def solve(u, u0, K11, K10, dof1, dof0, r1=None, ext0=None, solver=spsolve): - """Linear solution of equation system with optional given values of + r"""Linear solution of equation system with optional given values of unknowns at prescribed deegrees of freedom. - K_11 du_1 = -r1 - K10 (u0_ext - u0) - + Notes + ----- + + .. math:: + + \boldsymbol{K}_{11}\ d\boldsymbol{u}_1 = + -\boldsymbol{r}_1 - \boldsymbol{K}_{10} ( + \boldsymbol{u}_{0,\text{ext}} - \boldsymbol{u}_0 + ) + + Examples + -------- + .. pyvista-plot:: + :force_static: + + >>> import felupe as fem + >>> + >>> mesh = fem.Rectangle(n=3) + >>> region = fem.RegionQuad(mesh) + >>> field = fem.FieldPlaneStrain(region, dim=2).as_container() + >>> boundaries, loadcase = fem.dof.uniaxial(field, clamped=True) + >>> umat = fem.NeoHooke(mu=1.0, bulk=2.0) + >>> solid = fem.SolidBody(umat, field) + >>> + >>> K = solid.assemble.matrix() + >>> r = solid.assemble.vector() + >>> + >>> dof0 = loadcase["dof0"] + >>> dof1 = loadcase["dof1"] + >>> ext0 = loadcase["ext0"] + >>> + >>> system = fem.solve.partition(field, K, dof1, dof0, r) + >>> field += fem.solve.solve(*system, ext0=ext0) + + See Also + -------- + felupe.solve.partition : Perform a partitioning into active and prescribed dof. """ # init active residuals From 46b50934b83087158e947f91d13e876ee685c191 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Sat, 15 Feb 2025 23:08:23 +0100 Subject: [PATCH 32/47] Raise a TypeError if a container is used in a region (#939) --- CHANGELOG.md | 1 + src/felupe/region/_region.py | 6 ++++++ src/felupe/region/_templates.py | 30 +++++++++++++++++++++++++++++ tests/test_region.py | 34 +++++++++++++++++++++++++++++++++ 4 files changed, 71 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 3f682575..1bbf3676 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,6 +12,7 @@ All notable changes to this project will be documented in this file. The format - Add a new argument to apply a callable on the assembled vector/matrix of a solid body, `SolidBody(..., apply=None)`. This may be used to sum the list of sub-blocks instead of stacking them together, `SolidBody(..., block=False, apply=None)`. This is useful for mixed formulations where both the deformation gradient and the displacement values are required. - Add support for non-symmetric bilinear mixed forms in `IntegralForm`. - Add `element.Element` to the top-level package namespace. +- Raise a TypeError if a mesh container is used as the `mesh`-argument in a region. The error message explains that the container is not supported. A mesh must be used instead. ### Changed - The first Piola-Kirchhoff stress tensor is evaluated if `ViewSolid(stress_type=None)`. diff --git a/src/felupe/region/_region.py b/src/felupe/region/_region.py index 6f4de780..71ba0842 100644 --- a/src/felupe/region/_region.py +++ b/src/felupe/region/_region.py @@ -310,6 +310,12 @@ def reload( region = self if mesh is not None: + + if "container" in type(mesh).__name__.lower(): + raise TypeError( + "A mesh container is not supported by a region, use a mesh instead." + ) + region.mesh = mesh if element is not None: diff --git a/src/felupe/region/_templates.py b/src/felupe/region/_templates.py index a316fdc4..3cdac385 100644 --- a/src/felupe/region/_templates.py +++ b/src/felupe/region/_templates.py @@ -82,6 +82,11 @@ class RegionQuad(Region): def __init__( self, mesh, quadrature=GaussLegendre(order=1, dim=2), grad=True, **kwargs ): + if "container" in type(mesh).__name__.lower(): + raise TypeError( + "A mesh container is not supported by a region, use a mesh instead." + ) + element = Quad() if len(mesh.cells.T) > 4: @@ -117,6 +122,11 @@ class RegionQuadraticQuad(Region): def __init__( self, mesh, quadrature=GaussLegendre(order=2, dim=2), grad=True, **kwargs ): + if "container" in type(mesh).__name__.lower(): + raise TypeError( + "A mesh container is not supported by a region, use a mesh instead." + ) + element = QuadraticQuad() if len(mesh.cells.T) > 8: @@ -358,6 +368,11 @@ class RegionHexahedron(Region): def __init__( self, mesh, quadrature=GaussLegendre(order=1, dim=3), grad=True, **kwargs ): + if "container" in type(mesh).__name__.lower(): + raise TypeError( + "A mesh container is not supported by a region, use a mesh instead." + ) + element = Hexahedron() if len(mesh.cells.T) > 8: @@ -443,6 +458,11 @@ class RegionQuadraticHexahedron(Region): def __init__( self, mesh, quadrature=GaussLegendre(order=2, dim=3), grad=True, **kwargs ): + if "container" in type(mesh).__name__.lower(): + raise TypeError( + "A mesh container is not supported by a region, use a mesh instead." + ) + element = QuadraticHexahedron() if len(mesh.cells.T) > 20: @@ -647,6 +667,11 @@ class RegionTriangle(Region): def __init__( self, mesh, quadrature=TriangleQuadrature(order=1), grad=True, **kwargs ): + if "container" in type(mesh).__name__.lower(): + raise TypeError( + "A mesh container is not supported by a region, use a mesh instead." + ) + element = Triangle() if len(mesh.cells.T) > 3: @@ -682,6 +707,11 @@ class RegionTetra(Region): """ def __init__(self, mesh, quadrature=TetraQuadrature(order=1), grad=True, **kwargs): + if "container" in type(mesh).__name__.lower(): + raise TypeError( + "A mesh container is not supported by a region, use a mesh instead." + ) + element = Tetra() if len(mesh.cells.T) > 4: diff --git a/tests/test_region.py b/tests/test_region.py index 93b1ccb1..b0cd7180 100644 --- a/tests/test_region.py +++ b/tests/test_region.py @@ -175,6 +175,40 @@ def test_region_negative_volumes_cells(): region = fem.RegionQuad(mesh) +def test_container_warning(): + + mesh = fem.MeshContainer([fem.Rectangle()]) + element = fem.Quad() + quadrature = fem.GaussLegendre(order=1, dim=2) + with pytest.raises(TypeError): + region = fem.Region(mesh, element, quadrature) + + mesh = fem.MeshContainer([fem.Rectangle()]) + with pytest.raises(TypeError): + region = fem.RegionQuad(mesh) + + mesh = fem.MeshContainer([fem.Rectangle().add_midpoints_edges()]) + with pytest.raises(TypeError): + region = fem.RegionQuadraticQuad(mesh) + + mesh = fem.MeshContainer([fem.Rectangle().triangulate()]) + with pytest.raises(TypeError): + region = fem.RegionTriangle(mesh) + + mesh = fem.MeshContainer([fem.Cube()]) + with pytest.raises(TypeError): + region = fem.RegionHexahedron(mesh) + + mesh = fem.MeshContainer([fem.Cube().add_midpoints_edges()]) + with pytest.raises(TypeError): + region = fem.RegionQuadraticHexahedron(mesh) + + mesh = fem.MeshContainer([fem.Cube().triangulate()]) + with pytest.raises(TypeError): + region = fem.RegionTetra(mesh) + + if __name__ == "__main__": test_region() test_region_negative_volumes_cells() + test_container_warning() From f0617295e751829311c8f8ad9bcd9103dd912718 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Sun, 16 Feb 2025 21:53:40 +0100 Subject: [PATCH 33/47] Add `MeshContainer.as_vertex_mesh()` and `Field.from_mesh_container()` (#940) * Add `MeshContainer.as_vertex_mesh()` * Fix unique cells * Update __init__.py * Update test_region.py * Update test_element.py * Update CHANGELOG.md * Add `Field.from_mesh_container(container)` to create a top-level vertex-mesh. * Update CHANGELOG.md * Add `"vertex"` celltype for PyVista conversion * Set the cell-type to `"vertex"` in `MeshContainer.as_vertex_mesh()` * Update test_field.py * Update test_field.py --- CHANGELOG.md | 4 +++ src/felupe/__init__.py | 5 +++- src/felupe/element/__init__.py | 2 ++ src/felupe/element/_vertex.py | 47 +++++++++++++++++++++++++++++++++ src/felupe/field/_base.py | 13 +++++++++ src/felupe/mesh/_container.py | 5 ++++ src/felupe/mesh/_convert.py | 1 + src/felupe/region/__init__.py | 2 ++ src/felupe/region/_templates.py | 11 ++++++++ tests/test_element.py | 17 ++++++++++++ tests/test_field.py | 29 ++++++++++++++++++++ tests/test_mesh.py | 4 +++ tests/test_region.py | 6 +++++ 13 files changed, 145 insertions(+), 1 deletion(-) create mode 100644 src/felupe/element/_vertex.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 1bbf3676..20223e39 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,6 +13,10 @@ All notable changes to this project will be documented in this file. The format - Add support for non-symmetric bilinear mixed forms in `IntegralForm`. - Add `element.Element` to the top-level package namespace. - Raise a TypeError if a mesh container is used as the `mesh`-argument in a region. The error message explains that the container is not supported. A mesh must be used instead. +- Add the `Vertex` element formulation. +- Add a vertex-region `RegionVertex`. +- Add `MeshContainer.as_vertex_mesh()` to create a merged vertex mesh from the meshes of the mesh container. +- Add `Field.from_mesh_container(container)` to create a top-level field on a vertex mesh. ### Changed - The first Piola-Kirchhoff stress tensor is evaluated if `ViewSolid(stress_type=None)`. diff --git a/src/felupe/__init__.py b/src/felupe/__init__.py index a5e16468..e1b623a6 100644 --- a/src/felupe/__init__.py +++ b/src/felupe/__init__.py @@ -58,7 +58,7 @@ TetraMINI, Triangle, TriangleMINI, - TriQuadraticHexahedron, + Vertex, ) from .field import ( Field, @@ -113,6 +113,7 @@ RegionTriangleMINI, RegionTriQuadraticHexahedron, RegionTriQuadraticHexahedronBoundary, + RegionVertex, ) from .tools import hello_world, newtonrhapson, project, runs_on, save, topoints from .view import ViewField, ViewMesh @@ -187,6 +188,7 @@ "Triangle", "TriangleMINI", "TriQuadraticHexahedron", + "Vertex", "Circle", "Cube", "Grid", @@ -235,6 +237,7 @@ "RegionTriangleMINI", "RegionTriQuadraticHexahedron", "RegionTriQuadraticHexahedronBoundary", + "RegionVertex", "newtonrhapson", "project", "save", diff --git a/src/felupe/element/__init__.py b/src/felupe/element/__init__.py index 9ecbc8f7..461d54ff 100644 --- a/src/felupe/element/__init__.py +++ b/src/felupe/element/__init__.py @@ -15,6 +15,7 @@ from ._quad import BiQuadraticQuad, ConstantQuad, Quad, QuadraticQuad from ._tetra import QuadraticTetra, Tetra, TetraMINI from ._triangle import QuadraticTriangle, Triangle, TriangleMINI +from ._vertex import Vertex __all__ = [ "Element", @@ -37,4 +38,5 @@ "lagrange_line", "lagrange_quad", "lagrange_hexahedron", + "Vertex", ] diff --git a/src/felupe/element/_vertex.py b/src/felupe/element/_vertex.py new file mode 100644 index 00000000..b180c353 --- /dev/null +++ b/src/felupe/element/_vertex.py @@ -0,0 +1,47 @@ +# -*- 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 Element + + +class Vertex(Element): + r"""A vertex element formulation with constant shape functions. + + Notes + ----- + The vertex element is defined by one point. + """ + + def __init__(self): + self.points = np.array([[0.0]]) + self.cells = np.arange(len(self.points)).reshape(1, -1) + self.cell_type = "vertex" + + def function(self, r): + "Return the shape functions at given coordinate (r)." + return np.ones(1) + + def gradient(self, r): + "Return the gradient of shape functions at given coordinate (r)." + return np.zeros((1, 1)) + + def hessian(self, rs): + "Return the hessian of shape functions at given coordinate (r)." + return np.zeros((1, 1, 1)) diff --git a/src/felupe/field/_base.py b/src/felupe/field/_base.py index 835f56b1..6b53714d 100644 --- a/src/felupe/field/_base.py +++ b/src/felupe/field/_base.py @@ -22,6 +22,7 @@ from ..math import identity from ..math import sym as symmetric +from ..region import RegionVertex from ._container import FieldContainer from ._indices import Indices @@ -138,6 +139,18 @@ def _indices_per_cell(self, cells, dim): return cai, ai + @classmethod + def from_mesh_container(cls, mesh_container, dim=None, values=0.0): + "Create a Field on a vertex mesh from a mesh container." + + mesh = mesh_container.as_vertex_mesh() + region = RegionVertex(mesh) + + if dim is None: + dim = mesh.dim + + return cls(region, dim=dim, values=values) + def grad(self, sym=False, dtype=None, out=None, order="C"): r"""Gradient as partial derivative of field values w.r.t. undeformed coordinates, evaluated at the integration points of all cells in the region. diff --git a/src/felupe/mesh/_container.py b/src/felupe/mesh/_container.py index 17baa2ad..8cf234ff 100644 --- a/src/felupe/mesh/_container.py +++ b/src/felupe/mesh/_container.py @@ -279,6 +279,11 @@ def as_meshio(self, combined=True, **kwargs): return meshio.Mesh(self.points, cells, **kwargs) + def as_vertex_mesh(self): + "Return a merged vertex-mesh." + cells = np.unique(np.concatenate([mesh.cells.ravel() for mesh in self.meshes])) + return Mesh(self.points, cells.reshape(-1, 1), cell_type="vertex") + def copy(self): "Return a deepcopy of the mesh container." return deepcopy(self) diff --git a/src/felupe/mesh/_convert.py b/src/felupe/mesh/_convert.py index 4b978647..3c8377f7 100644 --- a/src/felupe/mesh/_convert.py +++ b/src/felupe/mesh/_convert.py @@ -59,6 +59,7 @@ def cell_types(): import pyvista as pv cell_types = [ + ("vertex", pv.CellType.VERTEX), ("line", pv.CellType.LINE), ("triangle", pv.CellType.TRIANGLE), ("triangle6", pv.CellType.QUADRATIC_TRIANGLE), diff --git a/src/felupe/region/__init__.py b/src/felupe/region/__init__.py index 8569957b..4086467b 100644 --- a/src/felupe/region/__init__.py +++ b/src/felupe/region/__init__.py @@ -22,6 +22,7 @@ RegionTriangleMINI, RegionTriQuadraticHexahedron, RegionTriQuadraticHexahedronBoundary, + RegionVertex, ) __all__ = [ @@ -48,4 +49,5 @@ "RegionTriangleMINI", "RegionTriQuadraticHexahedron", "RegionTriQuadraticHexahedronBoundary", + "RegionVertex", ] diff --git a/src/felupe/region/_templates.py b/src/felupe/region/_templates.py index 3cdac385..ec25f8cb 100644 --- a/src/felupe/region/_templates.py +++ b/src/felupe/region/_templates.py @@ -32,6 +32,7 @@ Triangle, TriangleMINI, TriQuadraticHexahedron, + Vertex, ) from ..mesh import Mesh from ..quadrature import GaussLegendre, GaussLegendreBoundary @@ -852,3 +853,13 @@ class RegionQuadraticTetra(Region): def __init__(self, mesh, quadrature=TetraQuadrature(order=2), grad=True, **kwargs): element = QuadraticTetra() super().__init__(mesh, element, quadrature, grad=grad, **kwargs) + + +class RegionVertex(Region): + "A region with a vertex element." + + def __init__( + self, mesh, quadrature=GaussLegendre(order=0, dim=1), grad=False, **kwargs + ): + element = Vertex() + super().__init__(mesh, element, quadrature, grad=grad, **kwargs) diff --git a/tests/test_element.py b/tests/test_element.py index a8c78571..afc92a55 100644 --- a/tests/test_element.py +++ b/tests/test_element.py @@ -30,6 +30,20 @@ import felupe as fem +def test_vertex1(): + vertex1 = fem.element.Vertex() + + r = [0] + + h = vertex1.function(r) + dhdr = vertex1.gradient(r) + d2hdrdr = vertex1.hessian(r) + + assert h[0] == 1.0 + assert np.all(dhdr[0] == 0.0) + assert np.all(d2hdrdr == 0.0) + + def test_line2(): line2 = fem.element.Line() @@ -185,6 +199,7 @@ def test_tri3(): assert np.all(dhdr[0] == -1) assert np.all(d2hdrdr == 0) + def test_tri6(): tri6 = fem.element.QuadraticTriangle() @@ -273,6 +288,8 @@ def test_aol(): if __name__ == "__main__": + test_vertex1() + test_line2() test_line_lagrange() diff --git a/tests/test_field.py b/tests/test_field.py index 7b6d0135..a3d8f2be 100644 --- a/tests/test_field.py +++ b/tests/test_field.py @@ -263,6 +263,34 @@ def test_link(): assert field1[1].values is field2[1].values +def test_toplevel(): + meshes = [ + fem.Rectangle(n=3), + fem.Rectangle(n=3).translate(1, axis=0).triangulate(), + ] + container = fem.MeshContainer(meshes, merge=True) + field = fem.Field.from_mesh_container(container).as_container() + regions = [ + fem.RegionQuad(container.meshes[0]), + fem.RegionTriangle(container.meshes[1]), + ] + fields = [ + fem.FieldContainer([fem.FieldPlaneStrain(regions[0], dim=2)]), + fem.FieldContainer([fem.FieldPlaneStrain(regions[1], dim=2)]), + ] + boundaries, loadcase = fem.dof.uniaxial(field, clamped=True) + umats = [ + fem.LinearElastic(E=2.1e5, nu=0.3), + fem.LinearElastic(E=1.0, nu=0.3), + ] + solids = [ + fem.SolidBody(umat=umats[0], field=fields[0]), + fem.SolidBody(umat=umats[1], field=fields[1]), + ] + step = fem.Step(items=solids, boundaries=boundaries) + fem.Job(steps=[step]).evaluate(x0=field) + + if __name__ == "__main__": test_axi() test_3d() @@ -270,3 +298,4 @@ def test_link(): test_mixed_lagrange() test_view() test_link() + test_toplevel() diff --git a/tests/test_mesh.py b/tests/test_mesh.py index 4c66604c..835cb286 100644 --- a/tests/test_mesh.py +++ b/tests/test_mesh.py @@ -374,6 +374,10 @@ def test_container(): assert np.allclose(container_2[0].points, mesh_1.points) assert np.allclose(container_2[0].cells, mesh_1.cells) + mesh = container.as_vertex_mesh() + + assert mesh.cells.shape[1] == 1 + def test_read(filename="tests/mesh.bdf"): mesh = fem.mesh.read(filename=filename, dim=2)[0] diff --git a/tests/test_region.py b/tests/test_region.py index b0cd7180..47bdb0a9 100644 --- a/tests/test_region.py +++ b/tests/test_region.py @@ -32,6 +32,12 @@ def test_region(): + + mesh = fem.Point() + r = fem.RegionVertex(mesh) + + assert r.h.flags["C_CONTIGUOUS"] + mesh = fem.Rectangle() r = fem.RegionQuad(mesh, uniform=True) From e252eaa90106c3dc9254ba3fc5bee13922054b67 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Sun, 16 Feb 2025 22:13:19 +0100 Subject: [PATCH 34/47] Update _hello_world.py with the new `Field.from_mesh_container()` --- src/felupe/tools/_hello_world.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/src/felupe/tools/_hello_world.py b/src/felupe/tools/_hello_world.py index 1f070882..7333e299 100644 --- a/src/felupe/tools/_hello_world.py +++ b/src/felupe/tools/_hello_world.py @@ -136,10 +136,7 @@ def hello_world( f" fem.{mesh}(n=6).translate(1, axis=0),", "]", "container = fem.MeshContainer(meshes, merge=True)", - "mesh = container.stack()", - "", - f"region = fem.{region}(mesh)", - f"field = fem.FieldContainer([fem.{field}(region, dim={dim})])", + "field = fem.Field.from_mesh_container(container).as_container()", "", "regions = [", f" fem.{region}(container.meshes[0]),", From 3fb4c488282e8acfa2b80d2577e3508628b47c1a Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Mon, 17 Feb 2025 09:16:22 +0100 Subject: [PATCH 35/47] Add `Vertex` element and region to the API docs --- docs/felupe/element.rst | 6 ++++++ docs/felupe/region.rst | 6 ++++++ 2 files changed, 12 insertions(+) diff --git a/docs/felupe/element.rst b/docs/felupe/element.rst index ec91e7c2..0c72a6e0 100644 --- a/docs/felupe/element.rst +++ b/docs/felupe/element.rst @@ -15,6 +15,7 @@ This module provides classes for the finite element formulations. .. autosummary:: + Vertex Line Quad Hexahedron @@ -59,6 +60,11 @@ This module provides classes for the finite element formulations. :undoc-members: :show-inheritance: +.. autoclass:: felupe.Vertex + :members: + :undoc-members: + :show-inheritance: + .. autoclass:: felupe.Line :members: :undoc-members: diff --git a/docs/felupe/region.rst b/docs/felupe/region.rst index c51352df..5b6e7c51 100644 --- a/docs/felupe/region.rst +++ b/docs/felupe/region.rst @@ -33,6 +33,7 @@ This module contains the definition of a region as well as a boundary region alo RegionLagrange RegionQuadBoundary RegionHexahedronBoundary + RegionVertex **Detailed API Reference** @@ -127,6 +128,11 @@ This module contains the definition of a region as well as a boundary region alo :show-inheritance: .. autoclass:: felupe.RegionLagrange + :members: + :undoc-members: + :show-inheritance: + +.. autoclass:: felupe.RegionVertex :members: :undoc-members: :show-inheritance: \ No newline at end of file From 0f93bc5fc0e30b569dd40400044ab52169d27360 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Mon, 17 Feb 2025 09:39:08 +0100 Subject: [PATCH 36/47] Update composite.rst --- docs/howto/composite.rst | 81 ++++++++++++++++++++++------------------ 1 file changed, 45 insertions(+), 36 deletions(-) diff --git a/docs/howto/composite.rst b/docs/howto/composite.rst index 4efd00bf..f6774486 100644 --- a/docs/howto/composite.rst +++ b/docs/howto/composite.rst @@ -1,95 +1,104 @@ Composite Regions with Solid Bodies ----------------------------------- -This section demonstrates how to set up a problem with two regions, each associated to a seperated solid body. +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. .. pyvista-plot:: :context: import felupe as fem - import numpy as np - m = fem.Rectangle(n=21) + inner = fem.Rectangle(a=(-1, -1), b=(1, 1), n=(5, 5)).triangulate() + + lower = fem.Rectangle(a=(-3, -3), b=(3, -1), n=(13, 5)) + upper = fem.Rectangle(a=(-3, 1), b=(3, 3), n=(13, 5)) + left = fem.Rectangle(a=(-3, -1), b=(-1, 1), n=(5, 5)) + right = fem.Rectangle(a=(1, -1), b=(3, 1), n=(5, 5)) + + outer = fem.MeshContainer([lower, upper, left, right], merge=True).stack() + + container = fem.MeshContainer([inner, outer], merge=True) + +A top-level (vertex) field, which contains all the unknowns, is required for the +definition of the boundary conditions as well as for the evaluation of the job. -In a second step, sub-meshes are created. +.. note:: + Ensure to init the mesh container with ``merge=True``, otherwise the points-array of + the container will be empty. .. pyvista-plot:: :context: - # take some points from the inside for the fiber-reinforced area - eps = 1e-3 - mask = np.arange(m.npoints)[np.logical_and.reduce([ - m.points[:, 0] >= 0.3, - m.points[:, 0] <= 0.7 + eps, - m.points[:, 1] >= 0.3, - m.points[:, 1] <= 0.7 + eps, - ])] - - # copies of the mesh - mesh = [m.copy(), m.copy()] - - # create sub-meshes (fiber, matrix) - mesh[0].update(cells=m.cells[ np.all(np.isin(m.cells, mask), axis=1)]) - mesh[1].update(cells=m.cells[~np.all(np.isin(m.cells, mask), axis=1)]) + container = fem.MeshContainer([inner, outer], merge=True) + field = fem.Field.from_mesh_container(container).as_container() -This is followed by the creation of a global region/field and two sub-regions/sub-fields. +The sub-meshes are available in the global mesh container, on which the sub-fields are +created. .. pyvista-plot:: :context: - - # a global and two sub-regions - region = fem.RegionQuad(m) - regions = [fem.RegionQuad(me) for me in mesh] - - # a global and two sub-fields - field = fem.FieldContainer([fem.FieldPlaneStrain(region, dim=2)]) + + regions = [ + fem.RegionTriangle(container.meshes[0]), + fem.RegionQuad(container.meshes[1]), + ] fields = [ fem.FieldContainer([fem.FieldPlaneStrain(regions[0], dim=2)]), fem.FieldContainer([fem.FieldPlaneStrain(regions[1], dim=2)]), ] -The displacement boundaries are created on the total field. +The displacement boundaries are created on the top-level field. .. pyvista-plot:: :context: boundaries = dict( - fixed=fem.Boundary(field[0], fx=0), - move=fem.Boundary(field[0], fx=1), + fixed=fem.Boundary(field[0], fx=field.region.mesh.x.min()), + move=fem.Boundary(field[0], fx=field.region.mesh.x.max()), ) -The rubber is associated to a Neo-Hookean material formulation whereas the steel is modeled by a linear elastic material formulation. Due to the large rotation, its large-strain formulation is required. For each material a solid body is created. +The rubber is associated to a Neo-Hookean material formulation whereas the steel is +modeled by a linear elastic material formulation. Due to the large rotation, its +large-strain formulation is required. For each material a solid body is created. .. pyvista-plot:: :context: # two material model formulations - neo_hooke = fem.NeoHooke(mu=1, bulk=1) linear_elastic = fem.LinearElasticLargeStrain(E=100, nu=0.3) + neo_hooke = fem.NeoHooke(mu=1, bulk=1) # the solid bodies fiber = fem.SolidBody(umat=linear_elastic, field=fields[0]) matrix = fem.SolidBody(umat=neo_hooke, field=fields[1]) -A step is created and further added to a job. The global field must be passed to the ``x0`` argument during the evaluation of the job. Internally, all field values are linked automatically, i.e. they share their ``values`` attribute. +A step is created and further added to a job. The global field must be passed as the +``x0`` argument during the evaluation of the job. Internally, all field values are +linked automatically, i.e. they share their ``values`` array. .. pyvista-plot:: :context: :force_static: # prepare a step with substeps - move = fem.math.linsteps([0, 0.5], num=10) + move = fem.math.linsteps([0, 3], num=10) step = fem.Step( items=[matrix, fiber], ramp={boundaries["move"]: move}, - boundaries=boundaries + boundaries=boundaries, ) # take care of the x0-argument job = fem.Job(steps=[step]) job.evaluate(x0=field) - field.plot("Principal Values of Logarithmic Strain").show() + plotter = fields[0].plot( + "Principal Values of Logarithmic Strain", show_undeformed=False + ) + fields[1].plot( + "Principal Values of Logarithmic Strain", show_undeformed=False, plotter=plotter + ).show() From b5349e72bec71b4c55a35a3e46bf57e2e0bd0186 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Mon, 17 Feb 2025 10:15:27 +0100 Subject: [PATCH 37/47] fix badge links --- README.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 9a4ed0a9..fbf3c830 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) [![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) [![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. @@ -123,7 +123,7 @@ tox ``` # Scientific Publications -A list of articles in which [![FElupe](https://img.shields.io/badge/%F0%9F%94%8D-FElupe-white)](https://felupe.readthedocs.io/stable) is involved. If you use FElupe in your scientific work, please star this repository, cite it [![DOI](https://zenodo.org/badge/360657894.svg)](https://zenodo.org/badge/latestdoi/360657894) and add your publication to this list. +A list of articles in which [![FElupe](https://img.shields.io/badge/%F0%9F%94%8D-FElupe-white)](https://felupe.readthedocs.io/en/stable) is involved. If you use FElupe in your scientific work, please star this repository, cite it [![DOI](https://zenodo.org/badge/360657894.svg)](https://zenodo.org/badge/latestdoi/360657894) and add your publication to this list.
Expand the list... From 339eae80e3095d81ef5614bcd99405bd16465654 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Mon, 17 Feb 2025 10:34:32 +0100 Subject: [PATCH 38/47] Fix `FieldContainer.plot()` for `RegionVertex` (no gradient available) (#942) --- CHANGELOG.md | 1 + src/felupe/view/_field.py | 60 +++++++++++++++++++++------------------ 2 files changed, 33 insertions(+), 28 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 20223e39..383c08fc 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -32,6 +32,7 @@ All notable changes to this project will be documented in this file. The format - Fix `Boundary(..., mode="and")` by ignoring any undefined axis. - Fix `tools.hello_world(planestrain=True)` with the correct region `RegionQuad` for the plane-strain template. - Fix the cells-array in `RegionQuadraticQuadBoundary` for midpoints of edges. +- Fix `FieldContainer.plot()` for regions without the shape-function gradients available. ### Removed - Remove the unused `shape`-argument in `element.Element(shape)`. Adopt the arbitrary-lagrange element to use its own `dim`-argument. This simplifies the creation of custom finite element formulations. diff --git a/src/felupe/view/_field.py b/src/felupe/view/_field.py index dbc3575d..3c367958 100644 --- a/src/felupe/view/_field.py +++ b/src/felupe/view/_field.py @@ -81,34 +81,38 @@ def __init__( point_data_from_field = {} cell_data_from_field = {} - if project is None: - cell_data_from_field = { - "Deformation Gradient": field.evaluate.deformation_gradient() - .mean(-2) - .T, - "Logarithmic Strain": field.evaluate.strain(tensor=True, asvoigt=True) - .mean(-2) - .T, - "Principal Values of Logarithmic Strain": field.evaluate.strain( - tensor=False - ) - .mean(-2) - .T, - } - elif callable(project): - point_data_from_field = { - "Deformation Gradient": project( - field.evaluate.deformation_gradient(), field.region - ), - "Logarithmic Strain": project( - field.evaluate.strain(tensor=True, asvoigt=True), field.region - ), - "Principal Values of Logarithmic Strain": project( - field.evaluate.strain(tensor=False), field.region - ), - } - else: - raise TypeError("The project-argument must be callable or None.") + if hasattr(field.region, "dhdX"): + + if project is None: + cell_data_from_field = { + "Deformation Gradient": field.evaluate.deformation_gradient() + .mean(-2) + .T, + "Logarithmic Strain": field.evaluate.strain( + tensor=True, asvoigt=True + ) + .mean(-2) + .T, + "Principal Values of Logarithmic Strain": field.evaluate.strain( + tensor=False + ) + .mean(-2) + .T, + } + elif callable(project): + point_data_from_field = { + "Deformation Gradient": project( + field.evaluate.deformation_gradient(), field.region + ), + "Logarithmic Strain": project( + field.evaluate.strain(tensor=True, asvoigt=True), field.region + ), + "Principal Values of Logarithmic Strain": project( + field.evaluate.strain(tensor=False), field.region + ), + } + else: + raise TypeError("The project-argument must be callable or None.") point_data_from_field["Displacement"] = displacement(field) From bfc38f4672e62a98a6d69077deb845130ef4f612 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Mon, 17 Feb 2025 10:38:48 +0100 Subject: [PATCH 39/47] Update composite.rst --- docs/howto/composite.rst | 3 +++ 1 file changed, 3 insertions(+) diff --git a/docs/howto/composite.rst b/docs/howto/composite.rst index f6774486..6b45f4ad 100644 --- a/docs/howto/composite.rst +++ b/docs/howto/composite.rst @@ -30,9 +30,12 @@ definition of the boundary conditions as well as for the evaluation of the job. .. pyvista-plot:: :context: + :force_static: container = fem.MeshContainer([inner, outer], merge=True) field = fem.Field.from_mesh_container(container).as_container() + + field.plot(plotter=container.plot(), point_size=15, color="red").show() The sub-meshes are available in the global mesh container, on which the sub-fields are created. From 9543b097775aba183cec30d004fbfbf887772b36 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Mon, 17 Feb 2025 10:48:12 +0100 Subject: [PATCH 40/47] Update composite.rst --- docs/howto/composite.rst | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/docs/howto/composite.rst b/docs/howto/composite.rst index 6b45f4ad..da100e5c 100644 --- a/docs/howto/composite.rst +++ b/docs/howto/composite.rst @@ -99,9 +99,10 @@ linked automatically, i.e. they share their ``values`` array. job = fem.Job(steps=[step]) job.evaluate(x0=field) - plotter = fields[0].plot( - "Principal Values of Logarithmic Strain", show_undeformed=False - ) fields[1].plot( - "Principal Values of Logarithmic Strain", show_undeformed=False, plotter=plotter + "Principal Values of Logarithmic Strain", + show_undeformed=False, + plotter=fields[0].plot( + "Principal Values of Logarithmic Strain", show_undeformed=False + ), ).show() From f43d34a5656eca795227f3436e7498d512a9e511 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Mon, 17 Feb 2025 11:19:12 +0100 Subject: [PATCH 41/47] Update composite.rst --- docs/howto/composite.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/howto/composite.rst b/docs/howto/composite.rst index da100e5c..1d749569 100644 --- a/docs/howto/composite.rst +++ b/docs/howto/composite.rst @@ -1,5 +1,5 @@ -Composite Regions with Solid Bodies ------------------------------------ +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. From 6aa973c2cb7cc3c9f2d27176a9429774cf6880bc Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Mon, 17 Feb 2025 14:41:00 +0100 Subject: [PATCH 42/47] Add `GaussLobatto` quadrature scheme (#943) * Add `GaussLobatto` quadrature * Update test_quadrature.py * Update test_quadrature.py * Update _gauss_lobatto.py --- CHANGELOG.md | 1 + docs/felupe/quadrature.rst | 12 ++ src/felupe/__init__.py | 10 +- src/felupe/quadrature/__init__.py | 5 +- .../{_gausslegendre.py => _gauss_legendre.py} | 0 src/felupe/quadrature/_gauss_lobatto.py | 171 ++++++++++++++++++ tests/test_quadrature.py | 32 +++- 7 files changed, 226 insertions(+), 5 deletions(-) rename src/felupe/quadrature/{_gausslegendre.py => _gauss_legendre.py} (100%) create mode 100644 src/felupe/quadrature/_gauss_lobatto.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 383c08fc..002a9fbe 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -17,6 +17,7 @@ All notable changes to this project will be documented in this file. The format - Add a vertex-region `RegionVertex`. - Add `MeshContainer.as_vertex_mesh()` to create a merged vertex mesh from the meshes of the mesh container. - Add `Field.from_mesh_container(container)` to create a top-level field on a vertex mesh. +- Add `GaussLobatto` quadrature. The order-argument is equal to the 1d-sample points minus two. This ensures identical minimum order-arguments for Gauss-Legendre and Gauss-Lobatto schemes. ### Changed - The first Piola-Kirchhoff stress tensor is evaluated if `ViewSolid(stress_type=None)`. diff --git a/docs/felupe/quadrature.rst b/docs/felupe/quadrature.rst index d5132ab1..c781370e 100644 --- a/docs/felupe/quadrature.rst +++ b/docs/felupe/quadrature.rst @@ -11,6 +11,8 @@ This module contains quadrature (numeric integration) schemes for different fini GaussLegendre GaussLegendreBoundary + GaussLobatto + GaussLobattoBoundary **Triangles and Tetrahedrons** @@ -44,6 +46,16 @@ This module contains quadrature (numeric integration) schemes for different fini :undoc-members: :show-inheritance: +.. autoclass:: felupe.GaussLobatto + :members: + :undoc-members: + :show-inheritance: + +.. autoclass:: felupe.GaussLobattoBoundary + :members: + :undoc-members: + :show-inheritance: + .. autoclass:: felupe.quadrature.Triangle :members: :undoc-members: diff --git a/src/felupe/__init__.py b/src/felupe/__init__.py index e1b623a6..58108f9c 100644 --- a/src/felupe/__init__.py +++ b/src/felupe/__init__.py @@ -86,7 +86,13 @@ Step, ) from .mesh import Circle, Cube, Grid, Mesh, MeshContainer, Point, Rectangle -from .quadrature import BazantOh, GaussLegendre, GaussLegendreBoundary +from .quadrature import ( + BazantOh, + GaussLegendre, + GaussLegendreBoundary, + GaussLobatto, + GaussLobattoBoundary, +) from .quadrature import Tetrahedron as TetrahedronQuadrature from .quadrature import Triangle as TriangleQuadrature from .region import ( @@ -211,6 +217,8 @@ "MultiPointContact", "GaussLegendre", "GaussLegendreBoundary", + "GaussLobatto", + "GaussLobattoBoundary", "TetrahedronQuadrature", "TriangleQuadrature", "BazantOh", diff --git a/src/felupe/quadrature/__init__.py b/src/felupe/quadrature/__init__.py index 8ca019ab..158ff9b9 100644 --- a/src/felupe/quadrature/__init__.py +++ b/src/felupe/quadrature/__init__.py @@ -1,4 +1,5 @@ -from ._gausslegendre import GaussLegendre, GaussLegendreBoundary +from ._gauss_legendre import GaussLegendre, GaussLegendreBoundary +from ._gauss_lobatto import GaussLobatto, GaussLobattoBoundary from ._scheme import Scheme from ._sphere import BazantOh from ._tetra import Tetrahedron @@ -8,6 +9,8 @@ "Scheme", "GaussLegendre", "GaussLegendreBoundary", + "GaussLobatto", + "GaussLobattoBoundary", "Tetrahedron", "Triangle", "BazantOh", diff --git a/src/felupe/quadrature/_gausslegendre.py b/src/felupe/quadrature/_gauss_legendre.py similarity index 100% rename from src/felupe/quadrature/_gausslegendre.py rename to src/felupe/quadrature/_gauss_legendre.py diff --git a/src/felupe/quadrature/_gauss_lobatto.py b/src/felupe/quadrature/_gauss_lobatto.py new file mode 100644 index 00000000..c1223922 --- /dev/null +++ b/src/felupe/quadrature/_gauss_lobatto.py @@ -0,0 +1,171 @@ +# -*- 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 . +""" + +from string import ascii_lowercase + +import numpy as np + +from ._scheme import Scheme + + +def gauss_lobatto(deg): + r"""Gauss-Lobatto quadrature. + + Computes the sample points and weights for Gauss-Lobatto quadrature. These sample + points and weights will correctly integrate polynomials of degree + :math:`2 \cdot deg - 3` or less over the interval :math:`[-1, 1]` with the weight + function :math:`f(x) = 1`. + + Parameters + ---------- + deg : int + Number of sample points and weights. It must be >= 2. + + Returns + ------- + x : ndarray + 1-D ndarray containing the sample points. + y : ndarray + 1-D ndarray containing the weights. + """ + if deg == 2: + x = np.array([-1.0, 1.0]) + y = np.array([1.0, 1.0]) + + elif deg == 3: + x = np.array([-1.0, 0.0, 1.0]) + y = np.array([1.0, 4.0, 1.0]) / 3 + + elif deg == 4: + a = np.sqrt(0.2) + x = np.array([-1.0, -a, a, 1.0]) + y = np.array([1.0, 5.0, 5.0, 1.0]) / 6 + + elif deg == 5: + a = np.sqrt(3 / 7) + x = np.array([-1.0, -a, 0.0, a, 1.0]) + y = np.array([0.1, 49 / 90, 32 / 45, 49 / 90, 0.1]) + + elif deg == 6: + a = np.sqrt(1 / 3 - 2 * np.sqrt(7) / 21) + b = np.sqrt(1 / 3 + 2 * np.sqrt(7) / 21) + c = (14 + np.sqrt(7)) / 30 + d = (14 - np.sqrt(7)) / 30 + x = np.array([-1.0, -b, -a, a, b, 1.0]) + y = np.array([1 / 15, d, c, c, d, 1 / 15]) + + elif deg == 7: + a = np.sqrt(5 / 11 - 2 / 11 * np.sqrt(5 / 3)) + b = np.sqrt(5 / 11 + 2 / 11 * np.sqrt(5 / 3)) + c = (124 + 7 * np.sqrt(15)) / 350 + d = (124 - 7 * np.sqrt(15)) / 350 + x = np.array([-1.0, -b, -a, 0, a, b, 1.0]) + y = np.array([1 / 21, d, c, 256 / 525, c, d, 1 / 21]) + + else: + raise ValueError("deg must be a positive integer (2 <= deg <= 7)") + + return x, y + + +class GaussLobatto(Scheme): + r"""An arbitrary-`order` Gauss-Lobatto quadrature rule of dimension 1, 2 or 3 on + the interval :math:`[-1, 1]`. + + Parameters + ---------- + order : int + The number of sample points :math:`n` minus two. The quadrature rule integrates + degree :math:`2n-3` polynomials exactly. + dim : int + The dimension of the quadrature region. + permute : bool, optional + Permute the quadrature points according to the cell point orderings (default is + True). This is supported for two and three dimensions as well as first and + second order schemes. Otherwise this flag is silently ignored. + + Notes + ----- + + The approximation is given by + + .. math:: + + \int_{-1}^1 f(x) dx \approx \sum f(x_q) w_q + + with quadrature points :math:`x_q` and corresponding weights :math:`w_q`. + + """ + + def __init__(self, order: int, dim: int): + if dim not in [1, 2, 3]: + raise ValueError("Wrong dimension.") + + x, w = gauss_lobatto(2 + order) + + points = ( + np.stack(np.meshgrid(*([x] * dim), indexing="ij"))[::-1].reshape(dim, -1).T + ) + + idx = list(ascii_lowercase)[:dim] + weights = np.einsum(", ".join(idx), *([w] * dim)).ravel() + + super().__init__(points, weights) + + +class GaussLobattoBoundary(GaussLobatto): + r"""An arbitrary-`order` Gauss-Lobatto quadrature rule of `dim` 1, 2 or 3 on the + interval ``[-1, 1]``. + + Parameters + ---------- + order : int + The number of sample points :math:`n` minus two. The quadrature rule integrates + degree :math:`2n-3` polynomials exactly. + dim : int + The dimension of the quadrature region. + permute : bool, optional + Permute the quadrature points according to the cell point orderings (default is + True). + + Notes + ----- + + The approximation is given by + + .. math:: + + \int_{-1}^1 f(x) dx \approx \sum f(x_q) w_q + + with quadrature points :math:`x_q` and corresponding weights :math:`w_q`. + + """ + + def __init__(self, order: int, dim: int): + super().__init__(order=order, dim=dim - 1) + + # reset dimension + self.dim = dim + + if self.dim == 2 or self.dim == 3: + # quadrature points projected onto first edge of a quad + # or onto first face of a hexahedron + self.points = np.hstack((self.points, -np.ones((len(self.points), 1)))) + + else: + raise ValueError("Given dimension not implemented (must be 2 or 3).") diff --git a/tests/test_quadrature.py b/tests/test_quadrature.py index 721127f5..1e2a9e4f 100644 --- a/tests/test_quadrature.py +++ b/tests/test_quadrature.py @@ -1,11 +1,11 @@ # -*- coding: utf-8 -*- """ - _______ _______ ___ __ __ _______ _______ + _______ _______ ___ __ __ _______ _______ | || || | | | | || || | | ___|| ___|| | | | | || _ || ___| -| |___ | |___ | | | |_| || |_| || |___ +| |___ | |___ | | | |_| || |_| || |___ | ___|| ___|| |___ | || ___|| ___| -| | | |___ | || || | | |___ +| | | |___ | || || | | |___ |___| |_______||_______||_______||___| |_______| This file is part of felupe. @@ -86,6 +86,30 @@ def test_gausslegendre_boundary(): fem.GaussLegendreBoundary(order=1, dim=4) +def test_gausslobatto(): + for order in [0, 1, 2, 3, 4, 5]: + for dim in [1, 2, 3]: + fem.GaussLobatto(order=order, dim=dim) + + with pytest.raises(ValueError): + fem.GaussLobatto(order=1, dim=4) + + with pytest.raises(ValueError): + fem.GaussLobatto(order=9, dim=4) + + +def test_gausslobatto_boundary(): + for order in [0, 1, 2, 3, 4, 5]: + for dim in [2, 3]: + fem.GaussLobattoBoundary(order=order, dim=dim) + + with pytest.raises(ValueError): + fem.GaussLobattoBoundary(order=1, dim=4) + + with pytest.raises(ValueError): + fem.GaussLobattoBoundary(order=9, dim=4) + + def test_triangle(): q = fem.TriangleQuadrature(order=1) assert q.points.shape == (1, 2) @@ -129,6 +153,8 @@ def test_sphere(): if __name__ == "__main__": test_gausslegendre() test_gausslegendre_boundary() + test_gausslobatto_boundary() + test_gausslobatto() test_triangle() test_tetra() test_sphere() From 77917713bd544a607b0fbb30d77b83d1d0b63a2d Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Tue, 18 Feb 2025 13:23:54 +0100 Subject: [PATCH 43/47] Enhance docs/howto source formatting --- docs/howto/axi.rst | 54 ++++++++++++++++++++++------- docs/howto/items.rst | 24 +++++++++---- docs/howto/meshgen.rst | 6 ---- docs/howto/mixed.rst | 26 +++++++++++--- docs/howto/solid.rst | 49 +++++++++++++++++++------- docs/howto/umat.rst | 12 +++++-- docs/howto/umat_hyperelasticity.rst | 14 ++++---- 7 files changed, 134 insertions(+), 51 deletions(-) diff --git a/docs/howto/axi.rst b/docs/howto/axi.rst index ce50ab2d..0e84c8fd 100644 --- a/docs/howto/axi.rst +++ b/docs/howto/axi.rst @@ -1,24 +1,43 @@ -Axisymmetric Problems ---------------------- +Two-dimensional Problems +------------------------ -For axisymmetric analyses an axisymmetric vector-valued field has to be created for the in-plane displacements. +For plane-strain and axisymmetric problems a vector-valued field has to be created for +the two-dimensional in-plane displacement components. -.. code-block:: python +.. tab:: Axisymmetric + + .. code-block:: python + + import felupe as fem + + mesh = fem.Rectangle(n=3) + region = fem.RegionQuad(mesh) + displacement = fem.FieldAxisymmetric(region, dim=2) + +.. tab:: Plane-Strain + + .. code-block:: python - import felupe as fem + import felupe as fem - mesh = fem.Rectangle(n=3) - region = fem.RegionQuad(mesh) - u = fem.FieldAxisymmetric(region, dim=2) - field = fem.FieldContainer([u]) + mesh = fem.Rectangle(n=3) + region = fem.RegionQuad(mesh) + displacement = fem.FieldPlaneStrain(region, dim=2) + -Now it gets important: The 3x3 deformation gradient for an axisymmetric problem is obtained with :meth:`~felupe.FieldAxisymmetric.grad` or :meth:`~felupe.FieldAxisymmetric.extract` methods. For a two-dimensional :class:`~felupe.FieldAxisymmetric` the gradient is modified to return a three-dimensional gradient. +The 3x3 deformation gradient for axisymmetric and plane-strain two-dimensional problems +is obtained by the :meth:`~felupe.FieldAxisymmetric.grad` or +:meth:`~felupe.FieldAxisymmetric.extract` methods (same for +:class:`~felupe.FieldPlaneStrain`). For these two-dimensional fields the gradient is +modified to return a three-dimensional gradient. .. code-block:: python + field = fem.FieldContainer([displacement]) F = field.extract(grad=True, sym=False, add_identity=True) -For simplicity, let's use the isotropic hyperelastic :class:`~felupe.NeoHooke` material model formulation. +For simplicity, let's use the isotropic hyperelastic :class:`~felupe.NeoHooke` material +model formulation. .. code-block:: python @@ -26,7 +45,14 @@ For simplicity, let's use the isotropic hyperelastic :class:`~felupe.NeoHooke` m .. note:: - Internally, FElupe provides an adopted low-level :class:`~felupe.assembly.IntegralFormAxisymmetric` class for the integration and the sparse matrix assemblage of axisymmetric problems. It uses the additional information (e.g. radial coordinates at integration points) stored in :class:`~felupe.FieldAxisymmetric` to provide a consistent interface in comparison to :class:`~felupe.assembly.IntegralFormCartesian`. The top-level :class:`~felupe.IntegralForm` chooses the appropriate low-level integral form based on the kind of field inside the field container. + Internally, FElupe provides an adopted low-level + :class:`~felupe.assembly.IntegralFormAxisymmetric` class for the integration and the + sparse matrix assemblage of axisymmetric problems. It uses the additional + information (e.g. radial coordinates at integration points) stored in + :class:`~felupe.FieldAxisymmetric` to provide a consistent interface in comparison + to :class:`~felupe.assembly.IntegralFormCartesian`. The top-level + :class:`~felupe.IntegralForm` chooses the appropriate low-level integral form based + on the kind of field inside the field container. .. code-block:: python @@ -35,4 +61,6 @@ For simplicity, let's use the isotropic hyperelastic :class:`~felupe.NeoHooke` m r = fem.IntegralForm(umat.gradient(F), field, dA).assemble() K = fem.IntegralForm(umat.hessian(F), field, dA, field).assemble() -To sum up, for axisymmetric problems use :class:`~felupe.FieldAxisymmetric`. Of course, mixed-field formulations may also be used with axisymmetric (displacement) fields. +To sum up, for axisymmetric problems use :class:`~felupe.FieldAxisymmetric` and for +plane-strain problems use :class:`~felupe.FieldPlaneStrain`. Of course, mixed-field +formulations may also be used with axisymmetric or plane-strain (displacement) fields. diff --git a/docs/howto/items.rst b/docs/howto/items.rst index 0ec68a59..51088bfd 100644 --- a/docs/howto/items.rst +++ b/docs/howto/items.rst @@ -1,10 +1,20 @@ Create Custom Items ~~~~~~~~~~~~~~~~~~~ -Items are supported to be used in a :class:`~felupe.Step` and in :func:`~felupe.newtonrhapson`. They provide methods to assemble a vector and a sparse matrix. An item has to be created as a class which takes a :class:`~felupe.FieldContainer` as input argument. In its ``__init__(self, field)`` method, the helpers :class:`~felupe.mechanics.Assemble` and :class:`~felupe.mechanics.Results` have to be evaluated and added as attributes. Internal methods which assemble the sparse vector and matrix, optionally with an updated :class:`~felupe.FieldContainer` provided by the `field` argument, have to be linked to :class:`~felupe.mechanics.Assemble`. +Items are supported to be used in a :class:`~felupe.Step` and in +:func:`~felupe.newtonrhapson`. They provide methods to assemble a vector and a sparse +matrix. An item has to be created as a class which takes a +:class:`~felupe.FieldContainer` as input argument. In its ``__init__(self, field)`` +method, the helpers :class:`~felupe.mechanics.Assemble` and +:class:`~felupe.mechanics.Results` have to be evaluated and added as attributes. +Internal methods which assemble the sparse vector and matrix, optionally with an updated +:class:`~felupe.FieldContainer` provided by the `field` argument, have to be linked to +:class:`~felupe.mechanics.Assemble`. .. note:: - This is only a minimal working example for an item. For more details, see the sources of :class:`~felupe.SolidBody` or :class:`~felupe.SolidBodyNearlyIncompressible`. + This is only a minimal working example for an item. For more details, see the + sources of :class:`~felupe.SolidBody` or + :class:`~felupe.SolidBodyNearlyIncompressible`. .. pyvista-plot:: :context: @@ -25,11 +35,13 @@ Items are supported to be used in a :class:`~felupe.Step` and in :func:`~felupe. def _matrix(self, field=None, **kwargs): return csr_matrix(([0.0], ([0], [0])), shape=(1, 1)) -This item may be added to a basic script. +This item is now added to a basic script. .. note:: - The vector- and matrix-shapes are automatically increased to match the shape of the other items if necessary. Hence, a sparse matrix with shape ``(1, 1)`` is a valid choice for a vector or a matrix filled with zeros. + The vector- and matrix-shapes are automatically increased to match the shape of the + other items if necessary. Hence, a sparse matrix with shape ``(1, 1)`` is a valid + choice for a vector or a matrix filled with zeros. .. pyvista-plot:: :context: @@ -39,7 +51,7 @@ This item may be added to a basic script. boundaries, loadcase = fem.dof.uniaxial(field, clamped=True, move=1.0) solid = fem.SolidBody(umat=fem.NeoHooke(mu=1, bulk=2), field=field) - item = MyItem(field=field) + my_item = MyItem(field=field) - step = fem.Step(items=[solid, item], boundaries=boundaries) + step = fem.Step(items=[solid, my_item], boundaries=boundaries) job = fem.Job(steps=[step]).evaluate() \ No newline at end of file diff --git a/docs/howto/meshgen.rst b/docs/howto/meshgen.rst index 7b1fabd6..e24025ae 100644 --- a/docs/howto/meshgen.rst +++ b/docs/howto/meshgen.rst @@ -29,12 +29,6 @@ FElupe provides a simple mesh generation module :ref:`mesh `. A # view the mesh in an interactive window mesh.plot().show() - # take a screenshot of an off-screen view - img = mesh.screenshot( - filename="mesh.png", - transparent_background=True, - ) - .. pyvista-plot:: :context: :nofigs: diff --git a/docs/howto/mixed.rst b/docs/howto/mixed.rst index a21d1ad6..ada1e0a1 100644 --- a/docs/howto/mixed.rst +++ b/docs/howto/mixed.rst @@ -1,7 +1,13 @@ Mixed-Field Problems ~~~~~~~~~~~~~~~~~~~~ -FElupe supports mixed-field formulations in a similar way it can handle (default) single-field formulations. The definition of a mixed-field formulation is shown for the hydrostatic-volumetric selective three-field-variation with independend fields for displacements :math:`\boldsymbol{u}`, pressure :math:`p` and volume ratio :math:`J`. The total potential energy for nearly-incompressible hyperelasticity is formulated with a determinant-modified deformation gradient. The built-in Neo-Hookean material model is used as an argument of :class:`~felupe.ThreeFieldVariation` for mixed-field problems. +FElupe supports mixed-field formulations in a similar way it can handle (default) +single-field formulations. The definition of a mixed-field formulation is shown for the +hydrostatic-volumetric selective three-field-variation with independend fields for +displacements :math:`\boldsymbol{u}`, pressure :math:`p` and volume ratio :math:`J`. +The total potential energy for nearly-incompressible hyperelasticity is formulated with +a determinant-modified deformation gradient. The built-in Neo-Hookean material model is +used as an argument of :class:`~felupe.ThreeFieldVariation` for mixed-field problems. .. pyvista-plot:: :context: @@ -11,7 +17,15 @@ FElupe supports mixed-field formulations in a similar way it can handle (default neohooke = fem.constitution.NeoHooke(mu=1.0, bulk=5000.0) umat = fem.constitution.ThreeFieldVariation(neohooke) -Next, let's create a meshed cube for a Hood-Taylor element formulation. The family of Hood-Taylor elements have a pressure field which is one order lower than the displacement field. A Hood-Taylor Q2/P1 hexahedron element formulation is created, where a tri-quadratic continuous (Lagrange) 27-point per cell displacement formulation is used in combination with discontinuous (tetra) 4-point per cell formulations for the pressure and volume ratio fields. The mesh of the cube is converted to a tri-quadratic mesh for the displacement field. The tetra regions for the pressure and the volume ratio are created on a dual (disconnected) mesh for the generation of the discontinuous fields. +Next, let's create a meshed cube for a Hood-Taylor element formulation. The family of +Hood-Taylor elements have a pressure field which is one order lower than the +displacement field. A Hood-Taylor Q2/P1 hexahedron element formulation is created, +where a tri-quadratic continuous (Lagrange) 27-point per cell displacement formulation +is used in combination with discontinuous (tetra) 4-point per cell formulations for the +pressure and volume ratio fields. The mesh of the cube is converted to a tri-quadratic +mesh for the displacement field. The tetra regions for the pressure and the volume ratio +are created on a dual (disconnected) mesh for the generation of the discontinuous +fields. .. pyvista-plot:: :context: @@ -38,14 +52,18 @@ Next, let's create a meshed cube for a Hood-Taylor element formulation. The fami field = fem.FieldContainer(fields=[displacement, pressure, volumeratio]) solid = fem.SolidBody(umat=umat, field=field) -Boundary conditions are enforced on the displacement field. For the pre-defined loadcases like the clamped uniaxial compression, the boundaries are automatically applied on the first field. +Boundary conditions are enforced on the displacement field. For the pre-defined +loadcases like the clamped uniaxial compression, the boundaries are automatically +applied on the first field. .. pyvista-plot:: :context: boundaries, loadcase = fem.dof.uniaxial(field, clamped=True) -The Step and Job definitions are identical to ones used with single field formulations. The deformed cube is finally visualized by PyVista. The cell-based means of the maximum principal values of the logarithmic strain tensor are shown. +The Step and Job definitions are identical to ones used with single field formulations. +The deformed cube is finally visualized by PyVista. The cell-based means of the maximum +principal values of the logarithmic strain tensor are shown. .. pyvista-plot:: :context: diff --git a/docs/howto/solid.rst b/docs/howto/solid.rst index 0f347140..a125d52f 100644 --- a/docs/howto/solid.rst +++ b/docs/howto/solid.rst @@ -1,12 +1,16 @@ Solid Bodies as Items for Assembly ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -The mechanics submodule contains classes for the generation of solid bodies. Solid body objects are supported as items of a :class:`~felupe.Step` and the :func:`~felupe.newtonrhapson` procedure. +The mechanics submodule contains classes for the generation of solid bodies. Solid body +objects are supported as items of a :class:`~felupe.Step` and the +:func:`~felupe.newtonrhapson` function. Solid Body ---------- -The generation of internal force vectors and stiffness matrices of solid bodies are provided as assembly-methods of a :class:`~felupe.SolidBody` or a :class:`~felupe.SolidBodyNearlyIncompressible`. +The generation of internal force vectors and stiffness matrices of solid bodies are +provided as assembly-methods of a :class:`~felupe.SolidBody` or a +:class:`~felupe.SolidBodyNearlyIncompressible`. .. pyvista-plot:: :context: @@ -27,15 +31,24 @@ The generation of internal force vectors and stiffness matrices of solid bodies stiffness_matrix = body.assemble.matrix(field, parallel=False) -During assembly, several results are stored, e.g. the gradient of the strain energy density function per unit undeformed volume w.r.t. the deformation gradient (first Piola-Kirchhoff stress tensor). Other results are the deformation gradient or the fourth-order elasticity tensor associated to the first Piola-Kirchhoff stress tensor. +During assembly, several results are stored, e.g. the gradient of the strain energy +density function per unit undeformed volume w.r.t. the deformation gradient (first +Piola-Kirchhoff stress tensor). Other results are the deformation gradient or the +fourth-order elasticity tensor associated to the first Piola-Kirchhoff stress tensor. .. math:: - \boldsymbol{F} &= \frac{\partial \boldsymbol{x}}{\partial \boldsymbol{X}} + \boldsymbol{F} &= \frac{ + \partial \boldsymbol{x}}{\partial \boldsymbol{X} + } - \boldsymbol{P} &= \frac{\partial \psi(\boldsymbol{C}(\boldsymbol{F}))}{\partial \boldsymbol{F}} + \boldsymbol{P} &= \frac{ + \partial \psi(\boldsymbol{C}(\boldsymbol{F}))}{\partial \boldsymbol{F} + } - \mathbb{A} &= \frac{\partial^2 \psi(\boldsymbol{C}(\boldsymbol{F}))}{\partial \boldsymbol{F}\ \partial \boldsymbol{F}} + \mathbb{A} &= \frac{ + \partial^2 \psi(\boldsymbol{C}(\boldsymbol{F})) + }{\partial \boldsymbol{F}\ \partial \boldsymbol{F}} .. pyvista-plot:: :context: @@ -45,13 +58,19 @@ During assembly, several results are stored, e.g. the gradient of the strain ene A = body.results.elasticity -The Cauchy stress tensor, as well as the gradient and the hessian of the strain energy density function per unit undeformed volume are obtained by evaluate-methods of the solid body. +The Cauchy stress tensor, as well as the gradient and the hessian of the strain energy +density function per unit undeformed volume are obtained by evaluate-methods of the +solid body. .. math:: - \boldsymbol{P} &= \frac{\partial \psi(\boldsymbol{C}(\boldsymbol{F}))}{\partial \boldsymbol{F}} + \boldsymbol{P} &= \frac{ + \partial \psi(\boldsymbol{C}(\boldsymbol{F}))}{\partial \boldsymbol{F} + } - \mathbb{A} &= \frac{\partial^2 \psi(\boldsymbol{C}(\boldsymbol{F}))}{\partial \boldsymbol{F}\ \partial \boldsymbol{F}} + \mathbb{A} &= \frac{ + \partial^2 \psi(\boldsymbol{C}(\boldsymbol{F})) + }{\partial \boldsymbol{F}\ \partial \boldsymbol{F}} \boldsymbol{\sigma} &= \frac{1}{J} \boldsymbol{P} \boldsymbol{F}^T @@ -67,7 +86,8 @@ The Cauchy stress tensor, as well as the gradient and the hessian of the strain Body Force (Gravity) on a Solid Body ------------------------------------ -The generation of external force vectors or stiffness matrices of body forces acting on solid bodies are provided as assembly-methods of a :class:`~felupe.SolidBodyForce`. +The generation of external force vectors or stiffness matrices of body forces acting on +solid bodies are provided as assembly-methods of a :class:`~felupe.SolidBodyForce`. .. math:: @@ -86,11 +106,13 @@ The generation of external force vectors or stiffness matrices of body forces ac Pressure Boundary on a Solid Body --------------------------------- -The generation of force vectors or stiffness matrices of pressure boundaries on solid bodies are provided as assembly-methods of a :class:`~felupe.SolidBodyPressure`. +The generation of force vectors or stiffness matrices of pressure boundaries on solid +bodies are provided as assembly-methods of a :class:`~felupe.SolidBodyPressure`. .. math:: - \delta W_{ext} = \int_{\partial V} \delta \boldsymbol{u} \cdot (-p) \ J \boldsymbol{F}^{-T} \ d\boldsymbol{A} + \delta W_{ext} = \int_{\partial V} + \delta \boldsymbol{u} \cdot (-p) \ J \boldsymbol{F}^{-T} \ d\boldsymbol{A} .. pyvista-plot:: @@ -113,7 +135,8 @@ The generation of force vectors or stiffness matrices of pressure boundaries on ) -For axisymmetric problems the boundary region has to be created with the attribute ``ensure_3d=True``. +For axisymmetric problems the boundary region has to be created with the attribute +``ensure_3d=True``. .. pyvista-plot:: :context: diff --git a/docs/howto/umat.rst b/docs/howto/umat.rst index 36bffce7..b5b1fed4 100644 --- a/docs/howto/umat.rst +++ b/docs/howto/umat.rst @@ -1,7 +1,10 @@ Small-Strain based Materials ---------------------------- -A user material (``umat``) based on the incremental small-strain tensor, e.g. suitable for linear elastic-plastic material formulations, is provided by :class:`~felupe.MaterialStrain`. A user-defined function must be created with the arguments and must return: +A user material (``umat``) based on the incremental small-strain tensor, e.g. suitable +for linear elastic-plastic material formulations, is provided by +:class:`~felupe.MaterialStrain`. A user-defined function must be created with the +arguments and must return: +----------+---------------+---------------------------------------+ | **Kind** | **Symbol** | **Description** | @@ -26,7 +29,9 @@ A user material (``umat``) based on the incremental small-strain tensor, e.g. su def material(dε, εn, σn, ζn, **kwargs): return dσdε, σ, ζ -This function is further added as the ``material`` argument of :class:`~felupe.MaterialStrain`. If the material makes use of state variables, the shapes of these internal state variables must be provided. +This function is further added as the ``material`` argument of +:class:`~felupe.MaterialStrain`. If the material makes use of state variables, the +shapes of these internal state variables must be provided. .. code-block:: python @@ -34,7 +39,8 @@ This function is further added as the ``material`` argument of :class:`~felupe.M umat = fem.MaterialStrain(material=material, statevars=(0,), **kwargs) -FElupe contains two reference small-strain user materials, one for linear elastic materials and another one for linear elastic-plastic materials with isotropic hardening: +FElupe contains two reference small-strain user materials, one for linear elastic +materials and another one for linear elastic-plastic materials with isotropic hardening: * :func:`~felupe.linear_elastic` * :func:`~felupe.linear_elastic_plastic_isotropic_hardening` diff --git a/docs/howto/umat_hyperelasticity.rst b/docs/howto/umat_hyperelasticity.rst index af5eab1c..af605c75 100644 --- a/docs/howto/umat_hyperelasticity.rst +++ b/docs/howto/umat_hyperelasticity.rst @@ -26,11 +26,13 @@ This function is further added as the ``fun`` argument of :class:`~felupe.Hypere umat = fem.Hyperelastic(fun=strain_energy_function, **kwargs) -FElupe contains several reference implementations of hyperelastic user material formulations, like +FElupe contains several reference implementations of hyperelastic user material +formulations, like -* :func:`~felupe.neo_hooke`, -* :func:`~felupe.mooney_rivlin`, -* :func:`~felupe.yeoh` or -* :func:`~felupe.ogden`. +* :func:`~felupe.constitution.tensortrax.models.hyperelastic.neo_hooke`, +* :func:`~felupe.constitution.tensortrax.models.hyperelastic.mooney_rivlin`, +* :func:`~felupe.constitution.tensortrax.models.hyperelastic.yeoh` or +* :func:`~felupe.constitution.tensortrax.models.hyperelastic.ogden`. -A complete list of all available model formulations is available in the :ref:`hyperelasticity ` section of the API reference. +A complete list of all available model formulations is available in the +:ref:`hyperelasticity ` section of the API reference. From 2d829f1fc273b1c7942ed4ef628b7fe754b1ea1a Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Tue, 18 Feb 2025 16:00:10 +0100 Subject: [PATCH 44/47] Docs: Make more figures static --- src/felupe/mesh/_tools.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/felupe/mesh/_tools.py b/src/felupe/mesh/_tools.py index 5c6d1769..0159201a 100644 --- a/src/felupe/mesh/_tools.py +++ b/src/felupe/mesh/_tools.py @@ -235,7 +235,7 @@ def rotate(points, cells, cell_type, angle_deg, axis, center=None, mask=None): Rotate a rectangle in the xy-plane by 35 degree. .. pyvista-plot:: - :include-source: True + :force_static: >>> import felupe as fem >>> @@ -313,7 +313,7 @@ def revolve(points, cells, cell_type, n=11, phi=180, axis=0, expand_dim=True): Revolve a cylinder from a rectangle. .. pyvista-plot:: - :include-source: True + :force_static: >>> import felupe as fem >>> From 423808d526d3852dac81ed651c590a4e31f57568 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Mon, 3 Mar 2025 17:07:02 +0100 Subject: [PATCH 45/47] Enhance `mesh.interpolate_line(axis=None)` (#945) use a curve-progress variable if `axis=None` by default --- CHANGELOG.md | 1 + src/felupe/mesh/_interpolate.py | 97 ++++++++++++++++++++------------- tests/test_mesh.py | 20 ++++++- 3 files changed, 76 insertions(+), 42 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 002a9fbe..f70e2c63 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -62,6 +62,7 @@ All notable changes to this project will be documented in this file. The format - Change supported Python versions to 3.9 - 3.12. - Change the `dtype`-argument in `Region.astype(dtype)` from an optional to a required argument. - Make `tensortrax` an optional dependency (again). Now FElupe does only depend on NumPy and SciPy, all other extras are optional. +- Enhance `mesh.interpolate_line(mesh, xi, axis=None, ...)` which by default uses a curve-progress variable when `axis=None`. The curve progress evaluation points `xi` must be within `0 <= xi <= 1`. ### Fixed - Fix the number of points for non-disconnected dual meshes. This reduces the assembled (sparse) vector- and matrix-shapes, which are defined on mixed-fields. diff --git a/src/felupe/mesh/_interpolate.py b/src/felupe/mesh/_interpolate.py index f902fa13..ef0c3a5b 100644 --- a/src/felupe/mesh/_interpolate.py +++ b/src/felupe/mesh/_interpolate.py @@ -19,7 +19,7 @@ import numpy as np -def interpolate_line(mesh, xi, axis, Interpolator=None, **kwargs): +def interpolate_line(mesh, xi, axis=None, Interpolator=None, **kwargs): r"""Return an interpolated line mesh from an existing line mesh with provided interpolation points on a given axis. @@ -28,13 +28,15 @@ def interpolate_line(mesh, xi, axis, Interpolator=None, **kwargs): mesh : Mesh A :class:`~felupe.Mesh` with cell-type ``line``. xi : ndarray - Points at which to interpolate data. - axis : int - The axis of the points at which the data will be interpolated. + Points at which to interpolate data. If axis is None, a normalized curve- + progress must be provided (:math:`0 \le \xi \le 1`). + axis : int or None, optional + The axis of the points at which the data will be interpolated. If None, a + normalized curve-progress is used. Default is None. Interpolator : callable or None, optional An interpolator class (default is :class:`scipy.interpolate.PchipInterpolator`). **kwargs : dict, optional - Optional keyword arguments are passed to the interpolator. + Optional keyword arguments are passed to the Interpolator. Returns ------- @@ -51,53 +53,70 @@ def interpolate_line(mesh, xi, axis, Interpolator=None, **kwargs): >>> import felupe as fem >>> import numpy as np + >>> from scipy.interpolate import CubicSpline >>> - >>> mesh = fem.mesh.Line(n=5).expand(n=1) + >>> mesh = fem.mesh.Line(b=1.0, n=5).expand(n=1) >>> t = mesh.x.copy() - >>> mesh.points[:, 0] = np.sin(np.pi / 2 * t) - >>> mesh.points[:, 1] = np.cos(np.pi / 2 * t) + >>> mesh.points[:, 0] = np.sin(2 * np.pi * t) + >>> mesh.points[:, 1] = np.cos(2 * np.pi * t) >>> - >>> mesh.plot(style="points", color="black").show() - - .. pyvista-plot:: - :context: - :force_static: + >>> mesh_new = fem.mesh.interpolate_line( + ... mesh, xi=np.linspace(0, 1), Interpolator=CubicSpline, bc_type="periodic" + ... ) + >>> + >>> mesh_new.plot( + ... plotter=mesh.plot(style="points", color="red", point_size=15), + ... color="black", + ... ).show() - >>> mesh_new = fem.mesh.interpolate_line(mesh, xi=np.linspace(0, 1), axis=1) - >>> mesh_new.plot(style="points", color="black").show() """ - # line connectivity - # concatenate the first point of each cell and the last point of the last cell - line = np.concatenate([mesh.cells[:, :-1].ravel(), mesh.cells[-1, -1:]]) + if Interpolator is None: + from scipy.interpolate import PchipInterpolator as Interpolator - # independent spline variable - points = mesh.points[line, axis] - ascending = np.argsort(points) + if axis is None: # progress-based interpolation - # dependent spline variable(s) - mask = np.ones(mesh.dim, dtype=bool) - mask[axis] = False + distances = np.linalg.norm(np.diff(mesh.points, axis=0), axis=1) + progress = np.insert(np.cumsum(distances), 0, 0) + progress /= progress.max() - axes = np.arange(mesh.dim) - values = mesh.points[line.reshape(-1, 1), axes[mask]] + spline = Interpolator(progress, mesh.points, **kwargs) - # spline interpolator - if Interpolator is None: - from scipy.interpolate import PchipInterpolator as Interpolator + points_new = spline(xi) + cells_new = np.repeat(np.arange(len(xi)), 2)[1:-1].reshape(-1, 2) + + mesh_new = type(mesh)(points_new, cells_new, cell_type="line") + mesh_new.points_derivative = spline.derivative()(xi) + + else: # independent- and dependent-variables + + # line connectivity + # concatenate the first point of each cell and the last point of the last cell + line = np.concatenate([mesh.cells[:, :-1].ravel(), mesh.cells[-1, -1:]]) + + # independent spline variable + points = mesh.points[line, axis] + ascending = np.argsort(points) + + # dependent spline variable(s) + mask = np.ones(mesh.dim, dtype=bool) + mask[axis] = False + + axes = np.arange(mesh.dim) + values = mesh.points[line.reshape(-1, 1), axes[mask]] - # create a spline - spline = Interpolator(points[ascending], values[ascending], **kwargs) + # create a spline + spline = Interpolator(points[ascending], values[ascending], **kwargs) - # evaluation points for the independent spline variable - points_new = np.zeros((len(xi), mesh.dim)) - points_new[:, axis] = xi - points_new[:, axes[mask]] = spline(xi) + # evaluation points for the independent spline variable + points_new = np.zeros((len(xi), mesh.dim)) + points_new[:, axis] = xi + points_new[:, axes[mask]] = spline(xi) - cells_new = np.repeat(np.arange(len(xi)), 2)[1:-1].reshape(-1, 2) + cells_new = np.repeat(np.arange(len(xi)), 2)[1:-1].reshape(-1, 2) - mesh_new = type(mesh)(points_new, cells_new, cell_type="line") - mesh_new.points_derivative = np.zeros_like(mesh_new.points) - mesh_new.points_derivative[:, axes[mask]] = spline.derivative()(xi) + mesh_new = type(mesh)(points_new, cells_new, cell_type="line") + mesh_new.points_derivative = np.zeros_like(mesh_new.points) + mesh_new.points_derivative[:, axes[mask]] = spline.derivative()(xi) return mesh_new diff --git a/tests/test_mesh.py b/tests/test_mesh.py index 835cb286..c2eaab99 100644 --- a/tests/test_mesh.py +++ b/tests/test_mesh.py @@ -1,11 +1,11 @@ # -*- coding: utf-8 -*- """ - _______ _______ ___ __ __ _______ _______ + _______ _______ ___ __ __ _______ _______ | || || | | | | || || | | ___|| ___|| | | | | || _ || ___| -| |___ | |___ | | | |_| || |_| || |___ +| |___ | |___ | | | |_| || |_| || |___ | ___|| ___|| |___ | || ___|| ___| -| | | |___ | || || | | |___ +| | | |___ | || || | | |___ |___| |_______||_______||_______||___| |_______| This file is part of felupe. @@ -568,6 +568,7 @@ def test_expand(): def test_interpolate_line(): import felupe as fem + from scipy.interpolate import CubicSpline mesh = fem.mesh.Line(n=5).expand(n=1).expand(n=1) t = mesh.x.copy() @@ -581,6 +582,19 @@ def test_interpolate_line(): assert len(mesh_new.points_derivative) == len(xi) assert np.allclose(mesh_new.points_derivative[:, 1], 0) + mesh = fem.mesh.Line(n=5).expand(n=1).expand(n=1) + t = mesh.x.copy() + mesh.points[:, 0] = np.sin(2 * np.pi * t) + mesh.points[:, 1] = np.cos(2 * np.pi * t) + + xi = np.linspace(0, 1, 101) + mesh_new = fem.mesh.interpolate_line( + mesh, xi=xi, Interpolator=CubicSpline, bc_type="periodic" + ) + + assert mesh_new.npoints == len(xi) + assert len(mesh_new.points_derivative) == len(xi) + if __name__ == "__main__": test_meshes() From f8cf2173e907223e8516b1fa8ac4113fc796fadd Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Mon, 3 Mar 2025 17:08:09 +0100 Subject: [PATCH 46/47] Update CHANGELOG.md --- CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index f70e2c63..6de1ee52 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -28,6 +28,7 @@ All notable changes to this project will be documented in this file. The format - Link all field-values to the values of the first field if no other field is given in `FieldContainer.link()`. - Change the default arguments from `block=True` to `block=None` in `SolidBody.assemble.vector(block=None)` and `SolidBody.assemble.matrix(block=None)` and define `block` on creation, `SolidBody(..., block=True)` instead. - Integrate and assemble `None`-items in the `fun`-argument of `IntegralForm`. On integration, `None` is returned. `None` will be assembled to an emtpy sparse matrix. +- Enhance `mesh.interpolate_line(mesh, xi, axis=None, ...)` which by default uses a curve-progress variable when `axis=None`. The curve progress evaluation points `xi` must be within `0 <= xi <= 1`. ### Fixed - Fix `Boundary(..., mode="and")` by ignoring any undefined axis. @@ -62,7 +63,6 @@ All notable changes to this project will be documented in this file. The format - Change supported Python versions to 3.9 - 3.12. - Change the `dtype`-argument in `Region.astype(dtype)` from an optional to a required argument. - Make `tensortrax` an optional dependency (again). Now FElupe does only depend on NumPy and SciPy, all other extras are optional. -- Enhance `mesh.interpolate_line(mesh, xi, axis=None, ...)` which by default uses a curve-progress variable when `axis=None`. The curve progress evaluation points `xi` must be within `0 <= xi <= 1`. ### Fixed - Fix the number of points for non-disconnected dual meshes. This reduces the assembled (sparse) vector- and matrix-shapes, which are defined on mixed-fields. From eb7fb90ab168c1d916ef9e201b8933fc17c31300 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Tue, 4 Mar 2025 16:45:45 +0100 Subject: [PATCH 47/47] Set version tag to 9.2.0 (#946) * Add Python 3.13 * set version tag to 9.2.0 --- .github/workflows/upload-codecov.yml | 4 ++-- CHANGELOG.md | 2 ++ pyproject.toml | 1 + src/felupe/__about__.py | 2 +- 4 files changed, 6 insertions(+), 3 deletions(-) diff --git a/.github/workflows/upload-codecov.yml b/.github/workflows/upload-codecov.yml index 87c0570a..31a8439d 100644 --- a/.github/workflows/upload-codecov.yml +++ b/.github/workflows/upload-codecov.yml @@ -14,7 +14,7 @@ jobs: runs-on: ubuntu-latest strategy: matrix: - python-version: ["3.9", "3.12"] + python-version: ["3.9", "3.13"] steps: - uses: actions/setup-python@v5 with: @@ -25,7 +25,7 @@ jobs: pip install tox tox -- --cov felupe --cov-report xml --cov-report term - name: Upload coverage to Codecov - if: ${{ matrix.python-version == '3.12' }} + if: ${{ matrix.python-version == '3.13' }} uses: codecov/codecov-action@v4 with: token: ${{ secrets.CODECOV_TOKEN }} diff --git a/CHANGELOG.md b/CHANGELOG.md index 6de1ee52..e0931bd1 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.2.0] - 2025-03-04 + ### Added - Add `SolidBody.assemble.mass(density=1.0)` and `SolidBodyNearlyIncompressible.assemble.mass(density=1.0)` to assemble the mass matrix. - Add `SolidBody.evaluate.stress(field)` to evaluate the (first Piola-Kirchhoff) stress tensor (engineering stress in linear elasticity). diff --git a/pyproject.toml b/pyproject.toml index 5f96496a..112f393e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -42,6 +42,7 @@ classifiers = [ "Programming Language :: Python :: 3.10", "Programming Language :: Python :: 3.11", "Programming Language :: Python :: 3.12", + "Programming Language :: Python :: 3.13", "Topic :: Scientific/Engineering", "Topic :: Scientific/Engineering :: Mathematics", "Topic :: Utilities" diff --git a/src/felupe/__about__.py b/src/felupe/__about__.py index 982a458b..62ad3ee9 100644 --- a/src/felupe/__about__.py +++ b/src/felupe/__about__.py @@ -1 +1 @@ -__version__ = "9.2.0-dev" +__version__ = "9.2.0"