From 2aaf58f072929d43e1659c99239ad6ed774ce0ea Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Sun, 29 Jun 2025 21:46:35 +0200 Subject: [PATCH 1/4] Update __about__.py --- src/felupe/__about__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/felupe/__about__.py b/src/felupe/__about__.py index 1f28d710..62365d8c 100644 --- a/src/felupe/__about__.py +++ b/src/felupe/__about__.py @@ -1 +1 @@ -__version__ = "9.3.0" +__version__ = "9.4.0-dev" From a77fe0d778762003bce3365907592224306650c2 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Tue, 1 Jul 2025 17:37:31 +0200 Subject: [PATCH 2/4] Enhance docstrings for `GaussLobatto` and `Mesh.add_points()` --- src/felupe/mesh/_mesh.py | 2 +- src/felupe/quadrature/_gauss_lobatto.py | 99 ++++++++++++++++++++++++- 2 files changed, 99 insertions(+), 2 deletions(-) diff --git a/src/felupe/mesh/_mesh.py b/src/felupe/mesh/_mesh.py index 0ccb4ee9..d9b111e3 100644 --- a/src/felupe/mesh/_mesh.py +++ b/src/felupe/mesh/_mesh.py @@ -324,7 +324,7 @@ def imshow(self, *args, ax=None, **kwargs): return ax def add_points(self, points): - "Add additional points to the mesh." + "Add additional points to the mesh in-place." self.update(points=np.vstack([self.points, points])) def clear_points_without_cells(self): diff --git a/src/felupe/quadrature/_gauss_lobatto.py b/src/felupe/quadrature/_gauss_lobatto.py index c1223922..c9235e6a 100644 --- a/src/felupe/quadrature/_gauss_lobatto.py +++ b/src/felupe/quadrature/_gauss_lobatto.py @@ -101,7 +101,6 @@ class GaussLobatto(Scheme): Notes ----- - The approximation is given by .. math:: @@ -110,6 +109,61 @@ class GaussLobatto(Scheme): with quadrature points :math:`x_q` and corresponding weights :math:`w_q`. + Examples + -------- + .. note:: + + The quadrature points are not weighted in the plots because the end-points would + be almost invisible for the Gauss-Lobatto quadrature rule. + + .. pyvista-plot:: + :force_static: + + >>> import felupe as fem + >>> + >>> element = fem.Line() + >>> quadrature = fem.GaussLobatto(order=2, dim=1) + >>> quadrature.plot( + ... plotter=element.plot(add_point_labels=False, show_points=False), + ... weighted=False, + ... ).show() + + .. pyvista-plot:: + :force_static: + + >>> import felupe as fem + >>> + >>> element = fem.QuadraticQuad() + >>> quadrature = fem.GaussLobatto(order=2, dim=2) + >>> quadrature.plot( + ... plotter=element.plot(add_point_labels=False, show_points=False), + ... weighted=False, + ... ).show() + + .. pyvista-plot:: + :force_static: + + >>> import felupe as fem + >>> + >>> element = fem.QuadraticHexahedron() + >>> quadrature = fem.GaussLobatto(order=2, dim=3) + >>> quadrature.plot( + ... plotter=element.plot(add_point_labels=False, show_points=False), + ... weighted=False, + ... ).show() + + .. pyvista-plot:: + :force_static: + + >>> import felupe as fem + >>> + >>> element = fem.ArbitraryOrderLagrangeElement(order=5, dim=2) + >>> quadrature = fem.GaussLobatto(order=5, dim=2) + >>> quadrature.plot( + ... plotter=element.plot(add_point_labels=False, show_points=False), + ... weighted=False, + ... ).show() + """ def __init__(self, order: int, dim: int): @@ -154,6 +208,49 @@ class GaussLobattoBoundary(GaussLobatto): with quadrature points :math:`x_q` and corresponding weights :math:`w_q`. + Examples + -------- + .. note:: + + The quadrature points are not weighted in the plots because the end-points would + be almost invisible for the Gauss-Lobatto quadrature rule. + + .. pyvista-plot:: + :force_static: + + >>> import felupe as fem + >>> + >>> element = fem.QuadraticQuad() + >>> quadrature = fem.GaussLobattoBoundary(order=2, dim=2) + >>> quadrature.plot( + ... plotter=element.plot(add_point_labels=False, show_points=False), + ... weighted=False, + ... ).show() + + .. pyvista-plot:: + :force_static: + + >>> import felupe as fem + >>> + >>> element = fem.QuadraticHexahedron() + >>> quadrature = fem.GaussLobattoBoundary(order=2, dim=3) + >>> quadrature.plot( + ... plotter=element.plot(add_point_labels=False, show_points=False), + ... weighted=False, + ... ).show() + + .. pyvista-plot:: + :force_static: + + >>> import felupe as fem + >>> + >>> element = fem.ArbitraryOrderLagrangeElement(order=5, dim=3) + >>> quadrature = fem.GaussLobattoBoundary(order=5, dim=3) + >>> quadrature.plot( + ... plotter=element.plot(add_point_labels=False, show_points=False), + ... weighted=False, + ... ).show() + """ def __init__(self, order: int, dim: int): From e9b8da4cdb56f00e98c599ab89f5e4dacdfb74ea Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Sat, 12 Jul 2025 23:03:27 +0200 Subject: [PATCH 3/4] Update _linear_elastic.py --- src/felupe/constitution/small_strain/models/_linear_elastic.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/felupe/constitution/small_strain/models/_linear_elastic.py b/src/felupe/constitution/small_strain/models/_linear_elastic.py index dc7353f9..65583d7f 100644 --- a/src/felupe/constitution/small_strain/models/_linear_elastic.py +++ b/src/felupe/constitution/small_strain/models/_linear_elastic.py @@ -56,7 +56,7 @@ def linear_elastic(dε, εn, σn, ζn, λ, μ, **kwargs): :math:`\boldsymbol{\varepsilon} = \boldsymbol{\varepsilon}_n + \Delta\boldsymbol{\varepsilon}`. 3. Evaluation of the stress :math:`\boldsymbol{\sigma}` and the algorithmic - consistent tangent modulus :math:`\mathbb{C}` (=``dσdε``). + consistent tangent modulus :math:`\mathbb{C}` (=**dσdε**). .. math:: From 19aad0358d95a6e6d1e8f11601a091970a232c6e Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Mon, 14 Jul 2025 14:54:41 +0200 Subject: [PATCH 4/4] Remove `SolidBody._area_change` (#970) and add comments to the code --- src/felupe/mechanics/_solidbody.py | 54 ++++++++++++++++++++++++------ 1 file changed, 44 insertions(+), 10 deletions(-) diff --git a/src/felupe/mechanics/_solidbody.py b/src/felupe/mechanics/_solidbody.py index 4fd0419f..c50e374c 100644 --- a/src/felupe/mechanics/_solidbody.py +++ b/src/felupe/mechanics/_solidbody.py @@ -22,7 +22,6 @@ import numpy as np from ..assembly import IntegralForm -from ..constitution import AreaChange from ..math import det, dot, transpose from ..view import ViewSolid from ._helpers import Assemble, Evaluate, Results @@ -174,8 +173,7 @@ class SolidBody(Solid): Notes ----- - The total potential energy of internal forces is given in Eq. - :eq:`solidbody` + The total potential energy of internal forces is given in Eq. :eq:`solidbody` .. math:: :label: solidbody @@ -269,11 +267,20 @@ def __init__( self.results.kinematics = self._extract(self.field) if statevars is not None: + # state variables are located in the results self.results.statevars = statevars + else: + # init the shape of the state variables statevars_shape = (0,) + + # take the shape of the last item of the list if the provided umat has an + # x-attribute if hasattr(umat, "x"): statevars_shape = umat.x[-1].shape + + # create and fill the array of state variables with zeros + # state variables are located in the results self.results.statevars = np.zeros( ( *statevars_shape, @@ -296,10 +303,6 @@ def __init__( kirchhoff_stress=self._kirchhoff_stress, ) - self._area_change = AreaChange() - - self._form = IntegralForm - def _vector( self, field=None, @@ -322,13 +325,18 @@ def _vector( if apply is None: apply = self.apply + # evaluate the (first Piola-Kirchhoff) stress tensor and store it in the results self.results.stress = self._gradient(field, args=args, kwargs=kwargs) - self.results.force = self._form( + + # assemble the internal force vector + # optionally, use only the first n items for mixed-field formulations + self.results.force = IntegralForm( fun=self.results.stress[slice(items)], v=self.field, dV=self.field.region.dV, ).assemble(parallel=parallel, block=block) + # apply a callback on the assembled internal force vector if apply is not None: self.results.force = apply(self.results.force) @@ -356,29 +364,43 @@ def _matrix( if apply is None: apply = self.apply + # evaluate the fourth-order elasticity tensor and store it in the results + # (associated to the first Piola-Kirchhoff stress tensor, i.e. the partial + # derivative of the first Piola-Kirchhoff stress tensor w.r.t. the deformation + # gradient tensor) self.results.elasticity = self._hessian(field, args=args, kwargs=kwargs) - form = self._form( + # assemble the (sparse) tangent stiffness matrix + # optionally, use only the first n items for mixed-field formulations + form = IntegralForm( fun=self.results.elasticity[slice(items)], v=self.field, u=self.field, dV=self.field.region.dV, ) + # in a first step, integrate the weak-form and store the stiffness values + # (this dense array is only allocated once and will be re-used in all following + # evaluations) self.results.stiffness_values = form.integrate( parallel=parallel, out=self.results.stiffness_values ) + # finally, the dense array of stiffness-values is assembled into the sparse + # matrix self.results.stiffness = form.assemble( values=self.results.stiffness_values, block=block ) + # apply a callback on the assembled tangent stiffness matrix if apply is not None: self.results.stiffness = apply(self.results.stiffness) return self.results.stiffness def _extract(self, field): + "Evaluate and return the kinematics (the deformation gradient tensor)." + self.field = field self.results.kinematics = self.field.extract(out=self.results.kinematics) @@ -388,6 +410,7 @@ def _gradient(self, field=None, args=(), kwargs=None): if kwargs is None: kwargs = {} + # update the deformation gradient if field is not None: self.field = field self.results.kinematics = self._extract(self.field) @@ -395,11 +418,17 @@ def _gradient(self, field=None, args=(), kwargs=None): if "out" in inspect.signature(self.umat.gradient).parameters: kwargs["out"] = self.results.gradient + # evaluate the partial derivative of the strain energy density w.r.t. the + # deformation gradient gradient = self.umat.gradient( [*self.results.kinematics, self.results.statevars], *args, **kwargs ) + # store the first Piola-Kirchhoff stress tensor in the results self.results.gradient = gradient[0] + # store all gradients and the temporary state variables in the results + # the state variables are only updated to the public results-attribute if the + # solution is valid self.results.stress, self.results._statevars = gradient[:-1], gradient[-1] return self.results.stress @@ -415,9 +444,14 @@ def _hessian(self, field=None, args=(), kwargs=None): if "out" in inspect.signature(self.umat.hessian).parameters: kwargs["out"] = self.results.hessian + # evaluate the second partial derivative of the strain energy density w.r.t. the + # deformation gradient self.results.elasticity = self.umat.hessian( [*self.results.kinematics, self.results.statevars], *args, **kwargs ) + + # store the partial derivative of the first Piola-Kirchhoff stress tensor w.r.t. + # the deformation gradient in the results self.results.hessian = self.results.elasticity[0] return self.results.elasticity @@ -458,7 +492,7 @@ def _mass(self, density=None): field = self.field[0].as_container() dim = field[0].dim - form = self._form( + form = IntegralForm( fun=[density * np.eye(dim).reshape(dim, dim, 1, 1)], v=field, u=field,