Merge pull request #623 from adtzlr/add-composite-material

Add `Volumetric`, `CompositeMaterial` and `&`-operator for `ConstitutiveMaterial`

CHANGELOG.md

@@ -10,6 +10,8 @@ All notable changes to this project will be documented in this file. The format
 - Add `ViewMaterialIncompressible(umat)` to view force-stretch curves for incompressible uniaxial tension/compression, planar shear and equi-biaxial tension.
 - Add a base class for constitutive materials with methods `ConstitutiveMaterial.view(incompressible=False)`, `ConstitutiveMaterial.plot(incompressible=False)` and `ConstitutiveMaterial.screenshot(incompressible=False)`.
 - Add a dict-attribute with material parameters to all built-in materials, e.g. `NeoHooke.kwargs = {"mu": self.mu, "bulk": self.bulk}`.
+- Add `umat = CompositeMaterial(material, other_material)`.
+- Add `&`-operator to combine constitutive materials `umat = material & other_material`. Note that only the first material must contain state variables.
 
 ### Changed
 - Don't disconnect the dual mesh by default for regions `RegionQuadraticTriangle` and `RegionQuadraticTetra` in `FieldsMixed`.

docs/felupe/constitution.rst

@@ -64,6 +64,12 @@ There are many different pre-defined constitutive material formulations availabl
    ViewMaterial
    ViewMaterialIncompressible
 
+**Merge Constitutive Materials**
+
+.. autosummary::
+
+   CompositeMaterial
+
 **Linear-Elasticity**
 
 .. autosummary::
@@ -150,6 +156,11 @@ There are many different pre-defined constitutive material formulations availabl
    :undoc-members:
    :inherited-members:
 
+.. autoclass:: felupe.CompositeMaterial
+   :members:
+   :undoc-members:
+   :inherited-members:
+
 .. autoclass:: felupe.NeoHooke
    :members:
    :undoc-members:

src/felupe/__init__.py

@@ -28,10 +28,12 @@
     MaterialStrain,
     NeoHooke,
     NeoHookeCompressible,
+    Volumetric,
     OgdenRoxburgh,
     ThreeFieldVariation,
     ViewMaterial,
     ViewMaterialIncompressible,
+    CompositeMaterial,
     VolumeChange,
     arruda_boyce,
     extended_tube,

src/felupe/constitution/__init__.py

@@ -1,6 +1,6 @@
 from ._kinematics import AreaChange, LineChange, VolumeChange
 from ._mixed import ThreeFieldVariation
-from ._models_hyperelasticity import NeoHooke, NeoHookeCompressible
+from ._models_hyperelasticity import NeoHooke, NeoHookeCompressible, Volumetric
 from ._models_hyperelasticity_ad import (
     arruda_boyce,
     extended_tube,
@@ -23,7 +23,12 @@
 )
 from ._models_linear_elasticity_large_strain import LinearElasticLargeStrain
 from ._models_pseudo_elasticity import OgdenRoxburgh
-from ._preview import ConstitutiveMaterial, ViewMaterial, ViewMaterialIncompressible
+from ._preview import (
+    ConstitutiveMaterial,
+    ViewMaterial,
+    ViewMaterialIncompressible,
+    CompositeMaterial,
+)
 from ._user_materials import (
     LinearElasticPlasticIsotropicHardening,
     Material,
@@ -70,4 +75,5 @@
     "ViewMaterial",
     "ViewMaterialIncompressible",
     "ConstitutiveMaterial",
+    "CompositeMaterial",
 ]

src/felupe/constitution/_models_hyperelasticity.py

@@ -188,7 +188,11 @@ def __init__(self, mu=None, bulk=None, parallel=False):
         self.mu = mu
         self.bulk = bulk
 
-        self.kwargs = {"mu": self.mu}
+        self.kwargs = {}
+
+        if self.mu is not None:
+            self.kwargs["mu"] = self.mu
+
         if self.bulk is not None:
             self.kwargs["bulk"] = self.bulk
 
@@ -225,8 +229,10 @@ def function(self, x, mu=None, bulk=None):
 
         J = det(F)
         C = dot(transpose(F), F, parallel=self.parallel)
+        W = np.zeros_like(J)
 
-        W = mu / 2 * (J ** (-2 / 3) * trace(C) - 3)
+        if mu is not None:
+            W += mu / 2 * (J ** (-2 / 3) * trace(C) - 3)
 
         if bulk is not None:
             W += bulk * (J - 1) ** 2 / 2
@@ -257,9 +263,11 @@ def gradient(self, x, mu=None, bulk=None):
 
         J = det(F)
         iFT = transpose(inv(F, J))
+        P = np.zeros_like(F)
 
-        # "physical"-deviatoric (not math-deviatoric!) part of P
-        P = mu * (F - ddot(F, F, parallel=self.parallel) / 3 * iFT) * J ** (-2 / 3)
+        if mu is not None:
+            # "physical"-deviatoric (not math-deviatoric!) part of P
+            P += mu * (F - ddot(F, F, parallel=self.parallel) / 3 * iFT) * J ** (-2 / 3)
 
         if bulk is not None:
             # "physical"-volumetric (not math-volumetric!) part of P
@@ -291,17 +299,20 @@ def hessian(self, x, mu=None, bulk=None):
 
         J = det(F)
         iFT = transpose(inv(F, J))
-        eye = identity(F)
-
-        # "physical"-deviatoric (not math-deviatoric!) part of A4
-        FF = ddot(F, F, parallel=self.parallel)
-        A4 = (mu * J ** (-2 / 3)) * (
-            cdya_ik(eye, eye, parallel=self.parallel)
-            - 2 / 3 * dya(F, iFT, parallel=self.parallel)
-            - 2 / 3 * dya(iFT, F, parallel=self.parallel)
-            + 2 / 9 * FF * dya(iFT, iFT, parallel=self.parallel)
-            + 1 / 3 * FF * cdya_il(iFT, iFT, parallel=self.parallel)
-        )
+        A4 = np.zeros((*F.shape[:2], *F.shape[:2], *F.shape[-2:]))
+
+        if mu is not None:
+            eye = identity(F)
+
+            # "physical"-deviatoric (not math-deviatoric!) part of A4
+            FF = ddot(F, F, parallel=self.parallel)
+            A4 += (mu * J ** (-2 / 3)) * (
+                cdya_ik(eye, eye, parallel=self.parallel)
+                - 2 / 3 * dya(F, iFT, parallel=self.parallel)
+                - 2 / 3 * dya(iFT, F, parallel=self.parallel)
+                + 2 / 9 * FF * dya(iFT, iFT, parallel=self.parallel)
+                + 1 / 3 * FF * cdya_il(iFT, iFT, parallel=self.parallel)
+            )
 
         if bulk is not None:
             p = bulk * (J - 1)
@@ -496,3 +507,10 @@ def hessian(self, x, mu=None, lmbda=None):
             A4 += lmbda * (dya(iFT, iFT, parallel=self.parallel) - np.log(J) * iFTiFT)
 
         return [A4]
+
+
+class Volumetric(NeoHooke):
+    "Neo-Hookean material formulation with deactivated shear modulus."
+
+    def __init__(self, bulk, parallel=False):
+        super().__init__(mu=None, bulk=bulk, parallel=parallel)

src/felupe/constitution/_models_linear_elasticity.py

@@ -584,7 +584,7 @@ def gradient(self, x, E=None, nu=None):
 
         return [stress, statevars]
 
-    def hessian(self, x, E=None, nu=None, shape=(1, 1)):
+    def hessian(self, x=None, E=None, nu=None, shape=(1, 1)):
         """Evaluate the elasticity tensor from the deformation gradient.
 
         Parameters

src/felupe/constitution/_preview.py

@@ -80,7 +80,6 @@ def plot(self, ax=None, show_title=True, show_kwargs=True, **kwargs):
         if show_title:
             title = self.umat.__class__.__name__
             if hasattr(self.umat, "fun"):
-
                 fun = self.umat.fun
                 label = [fun.__class__.__name__]
                 if callable(fun):
@@ -594,3 +593,56 @@ def screenshot(self, filename="umat.png", incompressible=False, **kwargs):
         plt.close(fig)
 
         return ax
+
+    def __and__(self, other_material):
+        return CompositeMaterial(self, other_material)
+
+
+class CompositeMaterial(ConstitutiveMaterial):
+    """A composite material with two constitutive materials merged. State variables are
+    only considered for the first material.
+
+    Parameters
+    ----------
+    material : ConstitutiveMaterial
+        First constitutive material.
+    other_material : ConstitutiveMaterial
+        Second constitutive material.
+
+    Notes
+    -----
+    ..  warning::
+        Do not merge two constitutive materials with the same keys of material
+        parameters. In this case, the values of these material parameters are taken from
+        the first constitutive material.
+
+    Examples
+    --------
+    >>> import felupe as fem
+    >>>
+    >>> nh = fem.NeoHooke(mu=1.0)
+    >>> vol = fem.Volumetric(bulk=2.0)
+    >>> umat = nh & vol
+    >>> ax = umat.plot()
+
+    ..  image:: images/umat_composite.png
+        :width: 400px
+
+    """
+
+    def __init__(self, material, other_material):
+        self.materials = [material, other_material]
+        self.kwargs = {**other_material.kwargs, **material.kwargs}
+        self.x = material.x
+
+    def gradient(self, x, **kwargs):
+        gradients = [material.gradient(x, **kwargs) for material in self.materials]
+        nfields = len(x) - 1
+        P = [np.sum([grad[i] for grad in gradients], axis=0) for i in range(nfields)]
+        statevars_new = gradients[0][-1]
+        return [*P, statevars_new]
+
+    def hessian(self, x, **kwargs):
+        hessians = [material.hessian(x, **kwargs) for material in self.materials]
+        nfields = len(x) - 1
+        return [np.sum([hess[i] for hess in hessians], axis=0) for i in range(nfields)]

tests/test_constitution.py

@@ -338,6 +338,8 @@ def neo_hooke(C, mu=1):
                 filename=f"../docs/felupe/images/umat_{umat.fun.__name__}.png",
                 incompressible=True,
             )
+
+    ax = umat.plot()
 
 
 def test_umat_hyperelastic2():
@@ -369,7 +371,7 @@ def neo_hooke(F, mu=1):
     assert np.allclose(dsde, dsde2)
 
 
-def test_umat_viscoelastic():
+def test_umat_viscoelastic(savefig=False):
     r, x = pre(sym=False, add_identity=True, add_random=True)
     F = x[0]
 
@@ -396,13 +398,14 @@ def viscoelastic(C, Cin, mu, eta, dtime):
     umat = fem.Hyperelastic(
         fem.constitution.finite_strain_viscoelastic, nstatevars=6, **kwargs
     )
-    ax = umat.screenshot(
-        filename="../docs/felupe/images/umat_finite_strain_viscoelastic.png",
-        ux=fem.math.linsteps([1, 1.5, 1, 2, 1, 2.5, 1], num=15),
-        ps=None,
-        bx=None,
-        incompressible=True,
-    )
+    if savefig:
+        ax = umat.screenshot(
+            filename="../docs/felupe/images/umat_finite_strain_viscoelastic.png",
+            ux=fem.math.linsteps([1, 1.5, 1, 2, 1, 2.5, 1], num=15),
+            ps=None,
+            bx=None,
+            incompressible=True,
+        )
 
     s2, statevars_new = umat.gradient([F, statevars])
     dsde2 = umat.hessian([F, statevars])
@@ -492,6 +495,20 @@ def test_elpliso():
     dsde = umat.hessian([F, statevars])
 
 
+def test_composite():
+    r, x = pre(sym=False, add_identity=True)
+    F = x[0]
+
+    nh = fem.NeoHooke(mu=1.0)
+    vol = fem.Volumetric(bulk=2.0)
+    umat = nh & vol
+
+    ax = umat.plot()
+
+    dWdF, statevars_new = umat.gradient([F, None])
+    d2WdFdF, = umat.hessian([F, None])
+
+
 if __name__ == "__main__":
     test_nh()
     test_linear()
@@ -501,8 +518,9 @@ def test_elpliso():
     test_umat()
     test_umat_hyperelastic(savefig=True)
     test_umat_hyperelastic2()
-    test_umat_viscoelastic()
+    test_umat_viscoelastic(savefig=True)
     test_umat_viscoelastic2()
     test_umat_strain()
     test_umat_strain_plasticity()
     test_elpliso()
+    test_composite()