Add LinearElasticOrthotropic (#840)

* Add `LinearElasticOrthotropic`

* Update _linear_elastic_orthotropic.py

* Update test_constitution.py

* Update CHANGELOG.md

* Update test_constitution.py

* Update _linear_elastic_orthotropic.py

* Add `LinearElasticOrthotropic` to the tests

CHANGELOG.md

@@ -15,6 +15,7 @@ All notable changes to this project will be documented in this file. The format
 - Add an optional `mechanics.Assemble(..., multiplier=None)` argument which is used in external items like  `SolidBodyForce`, `SolidBodyGravity` and `PointLoad` and is applied in `newtonrhapson(items, ...)`.
 - Add a new submodule `view` which contains the `View...` classes like `ViewSolid` or `ViewField`, previously located at `tools._plot`.
 - Add the `grad`- and `hess`-arguments to the `reload()`- and `copy()`-methods of a `Region`, i.e. `Region.reload(grad=None, hess=None)`.
+- Add `LinearElasticOrthotropic`.
 
 ### Changed
 - Change the internal initialization of `field = Field(region, values=1, dtype=None)` values from `field.values = np.ones(shape) * values` to `field = np.full(shape, fill_value=values, dtype=dtype)`. This enforces `field = Field(region, values=1)` to return the gradient array with data-type `int` which was of type `float` before.

docs/felupe/constitution/core.rst

@@ -22,6 +22,7 @@ This page contains the core (hard-coded) constitutive material model formulation
    LinearElasticPlaneStress
    constitution.LinearElasticTensorNotation
    LinearElasticLargeStrain
+   LinearElasticOrthotropic
 
 **Plasticity**
 
@@ -93,6 +94,11 @@ This page contains the core (hard-coded) constitutive material model formulation
 
 .. autofunction:: felupe.linear_elastic_plastic_isotropic_hardening
 
+.. autoclass:: felupe.LinearElasticOrthotropic
+   :members:
+   :undoc-members:
+   :inherited-members:
+
 .. autoclass:: felupe.MaterialStrain
    :members:
    :undoc-members:

src/felupe/__init__.py

@@ -23,6 +23,7 @@
     Laplace,
     LinearElastic,
     LinearElasticLargeStrain,
+    LinearElasticOrthotropic,
     LinearElasticPlaneStrain,
     LinearElasticPlaneStress,
     LinearElasticPlasticIsotropicHardening,
@@ -171,6 +172,7 @@
     "Laplace",
     "LinearElastic",
     "LinearElasticLargeStrain",
+    "LinearElasticOrthotropic",
     "LinearElasticPlaneStrain",
     "LinearElasticPlaneStress",
     "LinearElasticPlasticIsotropicHardening",

src/felupe/constitution/__init__.py

@@ -33,6 +33,7 @@
 from .linear_elasticity import (
     LinearElastic,
     LinearElasticLargeStrain,
+    LinearElasticOrthotropic,
     LinearElasticPlaneStrain,
     LinearElasticPlaneStress,
     LinearElasticTensorNotation,
@@ -70,6 +71,7 @@
     "Laplace",
     "LinearElastic",
     "LinearElasticLargeStrain",
+    "LinearElasticOrthotropic",
     "LinearElasticPlaneStrain",
     "LinearElasticPlaneStress",
     "LinearElasticTensorNotation",

src/felupe/constitution/linear_elasticity/__init__.py

@@ -6,11 +6,13 @@
     LinearElasticTensorNotation,
 )
 from ._linear_elastic_large_strain import LinearElasticLargeStrain
+from ._linear_elastic_orthotropic import LinearElasticOrthotropic
 
 __all__ = [
     "lame_converter",
     "LinearElastic",
     "LinearElasticLargeStrain",
+    "LinearElasticOrthotropic",
     "LinearElasticPlaneStrain",
     "LinearElasticPlaneStress",
     "LinearElasticTensorNotation",

src/felupe/constitution/linear_elasticity/_linear_elastic_orthotropic.py

@@ -0,0 +1,204 @@
+# -*- coding: utf-8 -*-
+"""
+This file is part of FElupe.
+
+FElupe is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+FElupe is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with FElupe.  If not, see <http://www.gnu.org/licenses/>.
+"""
+
+import numpy as np
+
+from ...math import ddot, identity, transpose
+from .._base import ConstitutiveMaterial
+
+
+class LinearElasticOrthotropic(ConstitutiveMaterial):
+    r"""Orthotropic linear-elastic material formulation.
+
+    Parameters
+    ----------
+    E1 : float
+        Young's modulus.
+    E2 : float
+        Young's modulus.
+    E3 : float
+        Young's modulus.
+    nu12 : float
+        Poisson ratio.
+    nu23 : float
+        Poisson ratio.
+    nu13 : float
+        Poisson ratio.
+    G12 : float
+        Shear modulus.
+    G23 : float
+        Shear modulus.
+    G13 : float
+        Shear modulus.
+
+    Examples
+    --------
+    ..  pyvista-plot::
+        :context:
+
+        >>> import felupe as fem
+        >>>
+        >>> umat = fem.LinearElasticOrthotropic(
+        >>>     E1=1, E2=1, E3=1, nu12=0.3, nu23=0.3, nu13=0.3, G12=0.4, G23=0.4, G13=0.4
+        >>> )
+        >>> ax = umat.plot()
+
+    ..  pyvista-plot::
+        :include-source: False
+        :context:
+        :force_static:
+
+        >>> import pyvista as pv
+        >>>
+        >>> fig = ax.get_figure()
+        >>> chart = pv.ChartMPL(fig)
+        >>> chart.show()
+
+    """
+
+    def __init__(self, E1, E2, E3, nu12, nu23, nu13, G12, G23, G13):
+        self.E1 = E1
+        self.E2 = E1
+        self.E3 = E3
+
+        self.nu12 = nu12
+        self.nu23 = nu23
+        self.nu13 = nu13
+
+        self.G12 = G12
+        self.G23 = G23
+        self.G13 = G13
+
+        self.kwargs = {
+            "E1": self.E1,
+            "E2": self.E2,
+            "E3": self.E3,
+            "nu12": self.nu12,
+            "nu23": self.nu23,
+            "nu13": self.nu13,
+            "G12": self.G12,
+            "G23": self.G23,
+            "G13": self.G13,
+        }
+
+        # aliases for gradient and hessian
+        self.stress = self.gradient
+        self.elasticity = self.hessian
+
+        # initial variables for calling
+        # ``self.gradient(self.x)`` and ``self.hessian(self.x)``
+        self.x = [np.eye(3), np.zeros(0)]
+
+    def gradient(self, x):
+        r"""Evaluate the stress tensor (as a function of the deformation gradient).
+
+        Parameters
+        ----------
+        x : list of ndarray
+            List with Deformation gradient :math:`\boldsymbol{F}` (3x3) as first item.
+
+        Returns
+        -------
+        ndarray of shape (3, 3, ...)
+            Stress tensor
+
+        """
+
+        F, statevars = x[0], x[-1]
+
+        # convert the deformation gradient to strain
+        H = F - identity(F)
+        strain = (H + transpose(H)) / 2
+
+        # init stress
+        elast = self.hessian(x=x)[0]
+
+        return [ddot(elast, strain, mode=(4, 2)), statevars]
+
+    def hessian(self, x=None, shape=(1, 1), dtype=None):
+        r"""Evaluate the elasticity tensor. The Deformation gradient is only
+        used for the shape of the trailing axes.
+
+        Parameters
+        ----------
+        x : list of ndarray, optional
+            List with Deformation gradient :math:`\boldsymbol{F}` (3x3) as first item
+            (default is None).
+        shape : tuple of int, optional
+            Tuple with shape of the trailing axes (default is (1, 1)).
+
+        Returns
+        -------
+        ndarray of shape (3, 3, 3, 3, ...)
+            elasticity tensor
+
+        """
+
+        E1 = self.E1
+        E2 = self.E2
+        E3 = self.E3
+
+        nu12 = self.nu12
+        nu23 = self.nu23
+        nu13 = self.nu13
+
+        G12 = self.G12
+        G23 = self.G23
+        G13 = self.G13
+
+        if x is not None:
+            dtype = x[0].dtype
+
+        elast = np.zeros((3, 3, 3, 3, *shape), dtype=dtype)
+
+        nu21 = nu12 * E2 / E1
+        nu32 = nu23 * E3 / E2
+        nu31 = nu13 * E3 / E1
+
+        J = 1 / (1 - nu12 * nu21 - nu23 * nu32 - nu31 * nu13 - 2 * nu21 * nu32 * nu13)
+
+        elast[0, 0, 0, 0] = E1 * (1 - nu23 * nu32) * J
+        elast[1, 1, 1, 1] = E2 * (1 - nu13 * nu31) * J
+        elast[2, 2, 2, 2] = E3 * (1 - nu12 * nu21) * J
+
+        elast[0, 0, 1, 1] = elast[1, 1, 0, 0] = E1 * (nu21 + nu31 * nu23) * J
+        elast[0, 0, 2, 2] = elast[2, 2, 0, 0] = E1 * (nu31 + nu21 * nu32) * J
+        elast[1, 1, 2, 2] = elast[2, 2, 1, 1] = E2 * (nu32 + nu12 * nu31) * J
+
+        elast[
+            [0, 1, 0, 1],
+            [1, 0, 1, 0],
+            [0, 0, 1, 1],
+            [1, 1, 0, 0],
+        ] = G12
+
+        elast[
+            [0, 2, 0, 2],
+            [2, 0, 2, 0],
+            [0, 0, 2, 2],
+            [2, 2, 0, 0],
+        ] = G13
+
+        elast[
+            [1, 2, 1, 2],
+            [2, 1, 2, 1],
+            [1, 1, 2, 2],
+            [2, 2, 1, 1],
+        ] = G23
+
+        return [elast]

tests/test_constitution.py

@@ -150,6 +150,39 @@ def test_linear():
     assert np.allclose(*check_dsde)
 
 
+def test_linear_orthotropic():
+    r, F = pre(sym=False, add_identity=True)
+
+    for Material in [
+        (fem.constitution.LinearElasticOrthotropic, {}),
+    ]:
+        LinearElastic, kwargs = Material
+
+        # doi.org/10.2478/ace-2018-0027 (pine wood)
+        le = LinearElastic(
+            E1=6919,
+            E2=271,
+            E3=450,
+            nu12=0.388,
+            nu23=0.278,
+            nu13=0.375,
+            G12=262,
+            G23=34,
+            G13=354,
+            **kwargs,
+        )
+
+        stress = le.gradient(F)[:-1]
+        dsde = le.hessian(F)
+
+        assert le.elasticity()[0].shape[-2:] == (1, 1)
+
+        assert stress[0].shape == (3, 3, *F[0].shape[-2:])
+        assert dsde[0].shape == (3, 3, 3, 3, 1, 1)
+
+        assert np.allclose(stress, 0)
+
+
 def test_linear_planestress():
     r, F = pre(sym=False, add_identity=True)
     F = [F[0][:2][:, :2]]
@@ -708,6 +741,7 @@ def test_laplace():
 if __name__ == "__main__":
     test_nh()
     test_linear()
+    test_linear_orthotropic()
     test_linear_planestress()
     test_linear_planestrain()
     test_kinematics()