Merge pull request #650 from adtzlr/add-equivalent-von-mises

Add `math.equivalent_von_mises(A)`

CHANGELOG.md

@@ -9,6 +9,8 @@ All notable changes to this project will be documented in this file. The format
 - Add `NearlyIncompressible` as a simplified version of `ThreeFieldVariation`. A constitutive material formulation on the distortional part of a strain energy function in terms of the deformation gradient has to be provided, e.g. by `umat = NearlyIncompressible(NeoHooke(mu=1), bulk=5000)`.
 - Add optional kwargs to a job-callback `Job(callback=lambda stepnumber, substepnumber, substep, **kwargs: None, **kwargs)` and `CharacteristicCurve(callback=lambda stepnumber, substepnumber, substep, **kwargs: None, **kwargs)`.
 - Add `DiscreteGeometry` properties `x`, `y` and `z` to access the colums of the points-array.
+- Add new math-function `math.equivalent_von_mises(A)`.
+- Add the evaluation of the equivalent von Mises Cauchy stress as cell-data in `ViewSolid`, available as `Solid.plot("Equivalent of Cauchy Stress")`.
 
 ### Changed
 - Rename `Mesh.save()` to `Mesh.write()` and add `Mesh.save()` as an alias to `Mesh.write()`.

src/felupe/math/__init__.py

@@ -18,6 +18,7 @@
     eigh,
     eigvals,
     eigvalsh,
+    equivalent_von_mises,
     identity,
     inv,
     majortranspose,
@@ -54,6 +55,7 @@
     "eigh",
     "eigvals",
     "eigvalsh",
+    "equivalent_von_mises",
     "identity",
     "inv",
     "majortranspose",

src/felupe/math/_tensor.py

@@ -316,3 +316,27 @@ def reshape(A, shape, trailing_axes=2):
 def ravel(A, trailing_axes=2):
     ij, shape = np.split(A.shape, [-trailing_axes])
     return reshape(A, shape=np.prod(ij))
+
+
+def equivalent_von_mises(A):
+    r"""Return the Equivalent von Mises values of symmetric matrices.
+
+    ..  math::
+
+        \boldsymbol{A}_{v} = \sqrt{
+            \frac{3}{2} \text{dev}(\boldsymbol{A}) : \text{dev}(\boldsymbol{A})
+        }
+
+    Parameters
+    ----------
+    A : (M, M, ...) ndarray
+        Symmetric matrices for which the equivalent von Mises value will be computed.
+
+    Returns
+    -------
+    AvM : (...) ndarray
+        The equivalent von Mises values.
+    """
+
+    devA = dev(A)
+    return np.sqrt(3 / 2 * ddot(devA, devA))

src/felupe/mechanics/_solidbody.py

@@ -164,7 +164,7 @@ class SolidBody(Solid):
 
         \boldsymbol{K}_{\boldsymbol{u}\boldsymbol{u}} \cdot \delta \boldsymbol{u} +
             \boldsymbol{f}_\boldsymbol{u} = \boldsymbol{0}
-    
+
     ..  note::
         This class also supports ``umat`` with mixed-field formulations like
         :class:`~felupe.NearlyIncompressible` or :class:`~felupe.ThreeFieldVariation`.

src/felupe/tools/_plot.py

@@ -18,7 +18,7 @@
 
 import numpy as np
 
-from ..math import eigvalsh, tovoigt
+from ..math import eigvalsh, equivalent_von_mises, tovoigt
 from ..mechanics._job import (
     deformation_gradient,
     displacement,
@@ -169,6 +169,9 @@ def plot(
                     ]
                     data_label = data_label[20:]
 
+                if "Equivalent of " in data_label:
+                    component_labels_dict[1] = [""]
+
                 component_labels = np.arange(dim)
                 if dim in component_labels_dict.keys():
                     component_labels = component_labels_dict[dim]
@@ -448,6 +451,9 @@ def __init__(
             cell_data_from_solid[f"Principal Values of {stress_label}"] = (
                 eigvalsh(stress).mean(-2)[::-1].T
             )
+            cell_data_from_solid[f"Equivalent of {stress_label}"] = (
+                equivalent_von_mises(stress).mean(-2).T
+            )
 
         super().__init__(
             field=field,

tests/test_math.py

@@ -104,6 +104,7 @@ def test_math():
     assert fem.math.ddot(A, C, mode=(4, 2)).shape == (3, 3, 5, 7)
 
     assert fem.math.ddot(A, A, mode=(4, 4)).shape == (3, 3, 3, 3, 5, 7)
+    assert fem.math.equivalent_von_mises(C).shape == (5, 7)
 
     with pytest.raises(TypeError):
         fem.math.ddot(A, B, mode=(4, 3))

tests/test_plot.py

@@ -164,6 +164,7 @@ def test_model():
         mesh, field, solid = pre(n=3)
         view = fem.View(field, solid=solid)
         plotter = view.plot(
+            "Equivalent of Cauchy Stress",
             off_screen=True,
             notebook=True,
         )