Add constitution.lame_converter_orthotropic() (#904)

docs/felupe/constitution/tools.rst

@@ -29,6 +29,7 @@ This page contains tools and helpers for constitutive material formulations.
 .. autosummary::
 
    constitution.lame_converter
+   constitution.lame_converter_orthotropic
 
 **Detailed API Reference**
 
@@ -68,3 +69,5 @@ This page contains tools and helpers for constitutive material formulations.
    :inherited-members:
 
 .. autofunction:: felupe.constitution.lame_converter
+
+.. autofunction:: felupe.constitution.lame_converter_orthotropic

src/felupe/constitution/__init__.py

@@ -12,6 +12,7 @@
     LinearElasticPlaneStress,
     LinearElasticTensorNotation,
     lame_converter,
+    lame_converter_orthotropic,
 )
 from .poisson import Laplace
 from .small_strain import MaterialStrain
@@ -32,6 +33,7 @@
     "LinearElasticPlaneStress",
     "LinearElasticTensorNotation",
     "lame_converter",
+    "lame_converter_orthotropic",
     "OgdenRoxburgh",
     "LinearElasticPlasticIsotropicHardening",
     "Material",

src/felupe/constitution/linear_elasticity/__init__.py

@@ -1,4 +1,4 @@
-from ._lame_converter import lame_converter
+from ._lame_converter import lame_converter, lame_converter_orthotropic
 from ._linear_elastic import (
     LinearElastic,
     LinearElasticPlaneStrain,
@@ -10,6 +10,7 @@
 
 __all__ = [
     "lame_converter",
+    "lame_converter_orthotropic",
     "LinearElastic",
     "LinearElasticLargeStrain",
     "LinearElasticOrthotropic",

src/felupe/constitution/linear_elasticity/_lame_converter.py

@@ -15,6 +15,7 @@
 You should have received a copy of the GNU General Public License
 along with FElupe.  If not, see <http://www.gnu.org/licenses/>.
 """
+from numpy.linalg import inv
 
 
 def lame_converter(E, nu):
@@ -50,3 +51,73 @@ def lame_converter(E, nu):
     mu = E / (2 * (1 + nu))
 
     return lmbda, mu
+
+
+def lame_converter_orthotropic(E, nu, G):
+    r"""Convert elastic orthotropic material parameters to Lamé parameter.
+
+    Parameters
+    ----------
+    E : list of float
+        List of the three elastic moduli :math:`E_1, E_2, E_3`.
+    nu : list of float
+        List of three poisson ratios :math:`\nu_{12}, \nu_{23}, \nu_{31}`.
+    G : list of float
+        List of three shear moduli :math:`G_{12}, G_{23}, G_{31}`.
+
+    Returns
+    -------
+    lmbda : list of float
+        List of six (upper triangle) first Lamé parameters
+        :math:`\lambda_11, \lambda_12, \lambda_13, \lambda_22, \lambda_23, \lambda_33`.
+    mu : list of float
+        List of the three second Lamé parameters :math:`\mu_1,\mu_2, \mu3`.
+
+    """
+
+    # unpack orthotropic elastic material parameters
+    E1, E2, E3 = E
+    ν12, ν23, ν31 = nu
+    G12, G23, G31 = G
+
+    # orthotropic symmetry
+    ν21 = ν12 * E2 / E1
+    ν32 = ν23 * E3 / E2
+    ν13 = ν31 * E1 / E3
+
+    C = inv(
+        [
+            [1 / E1, -ν21 / E2, -ν31 / E3, 0, 0, 0],
+            [-ν12 / E1, 1 / E2, -ν32 / E3, 0, 0, 0],
+            [-ν13 / E1, -ν23 / E2, 1 / E3, 0, 0, 0],
+            [0, 0, 0, 1 / G12, 0, 0],
+            [0, 0, 0, 0, 1 / G23, 0],
+            [0, 0, 0, 0, 0, 1 / G31],
+        ]
+    )
+
+    # take the components from this matrix
+    #
+    # [λ11 + 2 * μ1, λ12, λ13, 0, 0, 0]
+    # [λ12, λ22 + 2 * μ2, λ23, 0, 0, 0]
+    # [λ13, λ23, λ33 + 2 * μ3, 0, 0, 0]
+    # [0, 0, 0, (μ1 + μ2) / 2, 0, 0]
+    # [0, 0, 0, 0, (μ2 + μ3) / 2, 0]
+    # [0, 0, 0, 0, 0, (μ1 + μ3) / 2]
+
+    λ12 = C[0, 1]
+    λ23 = C[1, 2]
+    λ13 = C[0, 2]
+
+    μ1 = C[3, 3] - C[4, 4] + C[5, 5]
+    μ2 = C[3, 3] + C[4, 4] - C[5, 5]
+    μ3 = -C[3, 3] + C[4, 4] + C[5, 5]
+
+    λ11 = C[0, 0] - 2 * μ1
+    λ22 = C[1, 1] - 2 * μ2
+    λ33 = C[2, 2] - 2 * μ3
+
+    lmbda = [λ11, λ12, λ13, λ22, λ23, λ33]
+    mu = [μ1, μ2, μ3]
+
+    return lmbda, mu

src/felupe/constitution/linear_elasticity/_linear_elastic_orthotropic.py

@@ -27,24 +27,12 @@ class LinearElasticOrthotropic(ConstitutiveMaterial):
 
     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.
+    E : float
+        Young's modulus (E1, E2, E3).
+    nu : float
+        Poisson ratio (nu12, nu23, n31).
+    G : float
+        Shear modulus (G12, G23, G31).
 
     Examples
     --------
@@ -54,8 +42,8 @@ class LinearElasticOrthotropic(ConstitutiveMaterial):
         >>> 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
-        >>> )
+        ...     E=[1, 1, 1], nu=[0.3, 0.3, 0.3], G=[0.4, 0.4, 0.4]
+        ... )
         >>> ax = umat.plot()
 
     ..  pyvista-plot::
@@ -71,30 +59,12 @@ class LinearElasticOrthotropic(ConstitutiveMaterial):
 
     """
 
-    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,
-        }
+    def __init__(self, E, nu, G):
+        self.E = E
+        self.nu = nu
+        self.G = G
+
+        self.kwargs = {"E": self.E, "nu": self.nu, "G": self.G}
 
         # aliases for gradient and hessian
         self.stress = self.gradient
@@ -149,17 +119,9 @@ def hessian(self, x=None, shape=(1, 1), dtype=None):
 
         """
 
-        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
+        E1, E2, E3 = self.E
+        nu12, nu23, nu31 = self.nu
+        G12, G23, G31 = self.G
 
         if x is not None:
             dtype = x[0].dtype
@@ -168,7 +130,7 @@ def hessian(self, x=None, shape=(1, 1), dtype=None):
 
         nu21 = nu12 * E2 / E1
         nu32 = nu23 * E3 / E2
-        nu31 = nu13 * E3 / E1
+        nu13 = nu31 * E1 / E3
 
         J = 1 / (1 - nu12 * nu21 - nu23 * nu32 - nu31 * nu13 - 2 * nu21 * nu32 * nu13)
 
@@ -192,7 +154,7 @@ def hessian(self, x=None, shape=(1, 1), dtype=None):
             [2, 0, 2, 0],
             [0, 0, 2, 2],
             [2, 2, 0, 0],
-        ] = G13
+        ] = G31
 
         elast[
             [1, 2, 1, 2],

src/felupe/constitution/tensortrax/models/hyperelastic/_saint_venant_kirchhoff_orthotropic.py

@@ -30,9 +30,10 @@ def saint_venant_kirchhoff_orthotropic(C, mu, lmbda, r1, r2, r3):
     C : tensortrax.Tensor
         Right Cauchy-Green deformation tensor.
     mu : list of float
-        List of the three shear moduli.
+        List of the three second Lamé parameters :math:`\mu_1,\mu_2, \mu3`.
     lmbda : list of float
-        The six (upper triangle) orthotropic moduli.
+        List of six (upper triangle) first Lamé parameters
+        :math:`\lambda_11, \lambda_12, \lambda_13, \lambda_22, \lambda_23, \lambda_33`.
     r1 : list of float
         First normal vector of planes of symmetry.
     r2 : list of float
@@ -65,16 +66,21 @@ def saint_venant_kirchhoff_orthotropic(C, mu, lmbda, r1, r2, r3):
         >>> import felupe.constitution.tensortrax as mat
         >>> import felupe as fem
         >>>
-        >>> r = fem.math.rotation_matrix(30, axis=2)
+        >>> r = fem.math.rotation_matrix(0, axis=2)
+        >>> lmbda, mu = fem.constitution.lame_converter_orthotropic(
+        ...     E=[10, 10, 10],
+        ...     nu=[0.3, 0.3, 0.3],
+        ...     G=[1, 1, 1],
+        ... )
         >>> umat = mat.Hyperelastic(
         ...     mat.models.hyperelastic.saint_venant_kirchhoff_orthotropic,
-        ...     mu=[1, 1, 1],
-        ...     lmbda=[20, 20, 20, 20, 20, 20],
+        ...     mu=mu,
+        ...     lmbda=lmbda,
         ...     r1=r[:, 0],
         ...     r2=r[:, 1],
         ...     r3=r[:, 2],
         ... )
-        >>> ax = umat.plot(incompressible=False)
+        >>> ax = umat.plot(ux=fem.math.linsteps([1, 1.1], num=10), ps=None, bx=None)
 
     ..  pyvista-plot::
         :include-source: False

tests/test_constitution.py

@@ -136,22 +136,22 @@ def test_linear():
 def test_linear_orthotropic():
     r, F = pre(sym=False, add_identity=True)
 
+    lmbda, mu = fem.constitution.lame_converter_orthotropic(
+        E=[6919, 271, 450],
+        nu=[0.388, 0.278, 0.375],
+        G=[262, 34, 354],
+    )
+
     for Material in [
         (fem.constitution.LinearElasticOrthotropic, {}),
     ]:
-        LinearElastic, kwargs = Material
+        LinearElasticOrtho, 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,
+        le = LinearElasticOrtho(
+            E=[6919, 271, 450],
+            nu=[0.388, 0.278, 0.375],
+            G=[262, 34, 354],
             **kwargs,
         )