Merge pull request #327 from adtzlr/fix-solid-incompr-axi

Fix incompressible solid: Wrong volume evaluation for axisymmetric fields

CHANGELOG.md

@@ -3,6 +3,11 @@ All notable changes to this project will be documented in this file. The format
 
 ## [Unreleased]
 
+## [5.3.1] - 2022-11-03
+
+### Fixed
+- Fix volume evaluation of (nearly) incompressible solids for axisymmetric fields.
+
 ## [5.3.0] - 2022-11-03
 
 ### Added

felupe/mechanics/_solidbody_incompressible.py

@@ -28,6 +28,7 @@
 import numpy as np
 
 from .._assembly import IntegralFormMixed
+from .._field import FieldAxisymmetric
 from ..constitution import AreaChange
 from ..math import dot, transpose, det, dya, ddot
 from ._helpers import Assemble, Evaluate, Results
@@ -60,7 +61,12 @@ def h(self, parallel=False, jit=False):
 
     def v(self):
         "Cell volumes of the deformed configuration."
-        return (det(self.F[0]) * self.field.region.dV).sum(0)
+        dV = self.field.region.dV
+        if isinstance(self.field[0], FieldAxisymmetric):
+            R = self.field[0].radius
+            dA = self.field.region.dV
+            dV = 2 * np.pi * R * dA
+        return (det(self.F[0]) * dV).sum(0)
 
 
 class SolidBodyNearlyIncompressible:
@@ -103,7 +109,13 @@ def __init__(self, umat, field, bulk, state=None):
         self._form = IntegralFormMixed
 
         # volume of undeformed configuration
-        self.V = self.field.region.dV.sum(0)
+        if isinstance(self.field[0], FieldAxisymmetric):
+            R = self.field[0].radius
+            dA = self.field.region.dV
+            dV = 2 * np.pi * R * dA
+        else:
+            dV = self.field.region.dV
+        self.V = dV.sum(0)
 
         self.results = Results(stress=True, elasticity=True)
 

pyproject.toml

@@ -1,6 +1,6 @@
 [project]
 name = "felupe"
-version = "5.4.0"
+version = "5.3.1"
 description = "Finite Element Analysis"
 readme = "README.md"
 requires-python = ">=3.6"

setup.cfg

@@ -1,6 +1,6 @@
 [metadata]
 name = felupe
-version = 5.4.0
+version = 5.3.1
 author = Andreas Dutzler
 author_email = [email protected]
 description = Finite Element Analysis

tests/test_mechanics.py

@@ -130,9 +130,9 @@ def pre(dim, bulk=2):
     return umat, fe.FieldContainer([u])
 
 
-def pre_axi():
+def pre_axi(bulk=2):
 
-    umat = fe.NeoHooke(mu=1, bulk=2)
+    umat = fe.NeoHooke(mu=1, bulk=bulk)
 
     m = fe.Rectangle(n=3)
     r = fe.RegionQuad(m)
@@ -267,6 +267,54 @@ def test_solidbody_incompressible():
 
 def test_solidbody_axi():
 
+    umat, u = pre_axi(bulk=None)
+    b = fe.SolidBodyNearlyIncompressible(umat=umat, field=u, bulk=5000)
+    b = fe.SolidBodyNearlyIncompressible(
+        umat=umat, field=u, bulk=5000, state=fe.StateNearlyIncompressible(u)
+    )
+
+    for parallel in [False, True]:
+        for jit in [False, True]:
+
+            kwargs = {"parallel": parallel, "jit": jit}
+
+            r1 = b.assemble.vector(u, **kwargs)
+            assert r1.shape == (18, 1)
+
+            r2 = b.assemble.vector(**kwargs)
+            assert np.allclose(r1.toarray(), r2.toarray())
+
+            K1 = b.assemble.matrix(u, **kwargs)
+            assert K1.shape == (18, 18)
+
+            K2 = b.assemble.matrix(**kwargs)
+            assert np.allclose(K1.toarray(), K2.toarray())
+
+            P1 = b.results.stress
+            P2 = b.evaluate.gradient()
+            P2 = b.evaluate.gradient(u)
+            assert np.allclose(P1, P2)
+
+            A1 = b.results.elasticity
+            A2 = b.evaluate.hessian()
+            A2 = b.evaluate.hessian(u)
+            assert np.allclose(A1, A2)
+
+            F1 = b.results.kinematics
+            F2 = b._extract(u)
+            assert np.allclose(F1, F2)
+
+            s1 = b.evaluate.cauchy_stress()
+            s2 = b.evaluate.cauchy_stress(u)
+            assert np.allclose(s1, s2)
+
+            t1 = b.evaluate.kirchhoff_stress()
+            t2 = b.evaluate.kirchhoff_stress(u)
+            assert np.allclose(t1, t2)
+
+
+def test_solidbody_axi_incompressible():
+
     umat, u = pre_axi()
     b = fe.SolidBody(umat=umat, field=u)