Merge pull request #61 from adtzlr/add-planestrain-planestress

Add plane-strain plane-stress-incompressible

README.md

@@ -64,7 +64,7 @@ A = Mat.hessian([defgrad])[0]
 ```
 
 ## Template classes for hyperelasticity
-matADi provides several simple template classes suitable for simple hyperelastic materials. Some common isotropic hyperelastic material formulations are located in `matadi.models` (see list below). These strain energy functions have to be passed as the `fun` argument into an instance of `MaterialHyperelastic`. Usage is exactly the same as described above. To convert a hyperelastic material based on the deformation gradient into a mixed three-field formulation suitable for nearly-incompressible behavior (*displacements*, *pressure* and *volume ratio*) an instance of a `MaterialHyperelastic` class has to be passed to `ThreeFieldVariation`.
+matADi provides several simple template classes suitable for simple hyperelastic materials. Some isotropic hyperelastic material formulations are located in `matadi.models` (see list below). These strain energy functions have to be passed as the `fun` argument into an instance of `MaterialHyperelastic`. Usage is exactly the same as described above. To convert a hyperelastic material based on the deformation gradient into a mixed three-field formulation suitable for nearly-incompressible behavior (*displacements*, *pressure* and *volume ratio*) an instance of a `MaterialHyperelastic` class has to be passed to `ThreeFieldVariation`. For *plane strain* or *plane stress* use `MaterialHyperelasticPlaneStrain` or `MaterialHyperelasticPlaneStressIncompressible` instead of `MaterialHyperelastic`.
 
 ```python
 

matadi/__about__.py

@@ -1 +1 @@
-__version__ = "0.0.21"
+__version__ = "0.0.22"

matadi/__init__.py

@@ -15,6 +15,8 @@
     ThreeFieldVariation,
     MaterialHyperelastic,
     MaterialComposite,
+    MaterialHyperelasticPlaneStrain,
+    MaterialHyperelasticPlaneStressIncompressible,
 )
 from ._lab import Lab
 

matadi/_templates.py

@@ -3,7 +3,7 @@
 import numpy as np
 
 from . import Material, Variable
-from .math import det
+from .math import det, horzcat, vertcat, zeros
 
 
 class ThreeFieldVariation:
@@ -39,6 +39,33 @@ def _fun_wrapper(self, x, **kwargs):
         return self.fun(x[0], **kwargs)
 
 
+class MaterialHyperelasticPlaneStrain:
+    def __init__(self, fun, **kwargs):
+        F = Variable("F", 2, 2)
+        self.x = [F]
+        self.fun = fun
+        self.kwargs = kwargs
+        self.W = Material(self.x, self._fun_wrapper, kwargs=self.kwargs)
+        self.function = self.W.function
+        self.gradient = self.W.gradient
+        self.hessian = self.W.hessian
+
+    def _fun_wrapper(self, x, **kwargs):
+        F = horzcat(vertcat(x[0], zeros(1, 2)), zeros(3, 1))
+        F[2, 2] = 1  # fixed thickness ratio `h / H = 1`
+        return self.fun(F, **kwargs)
+
+
+class MaterialHyperelasticPlaneStressIncompressible(MaterialHyperelasticPlaneStrain):
+    def __init__(self, fun, **kwargs):
+        super().__init__(fun, **kwargs)
+
+    def _fun_wrapper(self, x, **kwargs):
+        F = horzcat(vertcat(x[0], zeros(1, 2)), zeros(3, 1))
+        F[2, 2] = 1 / det(x[0])  # thickness ratio `h / H = 1 / (a / A)`
+        return self.fun(F, **kwargs)
+
+
 class MaterialComposite:
     def __init__(self, materials):
         "Composite Material as a sum of a list of materials."

matadi/math.py

@@ -52,7 +52,6 @@
     logic_not,
     floor,
     ceil,
-
     SX,
     DM,
     MX,

pyproject.toml

@@ -1,6 +1,6 @@
 [project]
 name = "matadi"
-version = "0.0.21"
+version = "0.0.22"
 description = "Material Definition with Automatic Differentiation"
 readme = "README.md"
 requires-python = ">=3.6"

setup.cfg

@@ -1,6 +1,6 @@
 [metadata]
 name = matadi
-version = 0.0.21
+version = 0.0.22
 author = Andreas Dutzler
 author_email = [email protected]
 description = Material Definition with Automatic Differentiation

tests/test_plane.py

@@ -0,0 +1,69 @@
+import numpy as np
+
+from matadi import (
+    MaterialHyperelasticPlaneStrain,
+    MaterialHyperelasticPlaneStressIncompressible,
+)
+from matadi.models import neo_hooke
+
+
+def pre():
+
+    FF = np.random.rand(2, 2, 8, 1000)
+    for a in range(2):
+        FF[a, a] += 1
+
+    return FF
+
+
+def test_plane_strain():
+
+    # data
+    FF = pre()
+
+    # init Material
+    W = MaterialHyperelasticPlaneStrain(
+        fun=neo_hooke,
+        C10=0.5,
+    )
+
+    W0 = W.function([FF])
+    dW = W.gradient([FF])
+    DW = W.hessian([FF])
+
+    assert W0[0].shape == (8, 1000)
+    assert dW[0].shape == (2, 2, 8, 1000)
+    assert DW[0].shape == (2, 2, 2, 2, 8, 1000)
+
+    assert W0[0].shape == (8, 1000)
+    assert dW[0].shape == (2, 2, 8, 1000)
+    assert DW[0].shape == (2, 2, 2, 2, 8, 1000)
+
+
+def test_plane_stress():
+
+    # data
+    FF = pre()
+
+    # init Material
+    W = MaterialHyperelasticPlaneStressIncompressible(
+        fun=neo_hooke,
+        C10=0.5,
+    )
+
+    W0 = W.function([FF])
+    dW = W.gradient([FF])
+    DW = W.hessian([FF])
+
+    assert W0[0].shape == (8, 1000)
+    assert dW[0].shape == (2, 2, 8, 1000)
+    assert DW[0].shape == (2, 2, 2, 2, 8, 1000)
+
+    assert W0[0].shape == (8, 1000)
+    assert dW[0].shape == (2, 2, 8, 1000)
+    assert DW[0].shape == (2, 2, 2, 2, 8, 1000)
+
+
+if __name__ == "__main__":
+    test_plane_strain()
+    test_plane_stress()