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,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
 

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,15 +89,15 @@
 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
 # plotted with matplotlib.
 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-",
 )

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})

src/felupe/mechanics/_helpers.py

@@ -36,14 +36,21 @@ 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
 
         if cauchy_stress is not 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."

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,

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

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,20 +516,18 @@ 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()
 
     field = fem.FieldContainer([fem.Field(region, dim=2)])
     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():

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