Update Kinematics class and add new tests

Author: jpsferreira

hyper_surrogate/kinematics.py

@@ -68,41 +68,42 @@ def left_cauchy_green(f: np.ndarray) -> Any:
         # and 'jk' are the indices for the second matrix (transposed to 'kj' for multiplication).
         return np.einsum("nij,nkj->nik", f, f)
 
-    def strain_tensor(self, f: np.ndarray) -> Any:
+    @staticmethod
+    def strain_tensor(f: np.ndarray) -> Any:
         """
         Compute the strain tensor.
 
         Args:
-            f (np.ndarray): The deformation gradient.
+            f (np.ndarray): The deformation gradients. batched with shape (N, 3, 3).
 
         Returns:
-            np.ndarray: The strain tensor.
+            np.ndarray: The strain tensors. batched with shape (N, 3, 3).
         """
-        return 0.5 * (f.T @ f - np.eye(3))
+        return 0.5 * (np.einsum("nji,njk->nik", f, f) - np.eye(3))
 
     def stretch_tensor(self, f: np.ndarray) -> Any:
         """
         Compute the stretch tensor.
 
         Args:
-            f (np.ndarray): The deformation gradient.
+            f (np.ndarray): The deformation gradients. batched with shape (N, 3, 3).
 
         Returns:
-            np.ndarray: The stretch tensor.
+            np.ndarray: The stretch tensor. batched with shape (N, 3, 3).
         """
-        return np.sqrt(self.right_cauchy_green(f))
+        return np.sqrt(np.einsum("nji,njk->nik", f, f))
 
-    def rotation_tensor(self, f: np.ndarray) -> np.ndarray:
+    def rotation_tensor(self, f: np.ndarray) -> Any:
         """
-        Compute the rotation tensor.
+        Compute the rotation tensors.
 
         Args:
-            f (np.ndarray): The deformation gradient.
+            f (np.ndarray): The deformation gradients. batched with shape (N, 3, 3).
 
         Returns:
-            np.ndarray: The rotation tensor.
+            np.ndarray: The rotation tensors. batched with shape (N, 3, 3).
         """
-        return f @ np.linalg.inv(self.stretch_tensor(f))
+        return np.einsum("nij,njk->nik", f, np.linalg.inv(f))
 
     def principal_stretches(self, f: np.ndarray) -> np.ndarray:
         """

hyper_surrogate/symbolic.py

@@ -25,7 +25,7 @@ def _c_tensor(self) -> sym.Matrix:
         return sym.Matrix(3, 3, lambda i, j: sym.Symbol(f"C_{i+1}{j+1}"))
 
     # multuply c_tensor by itself
-    def _c_tensor_squared(self) -> sym.Matrix:
+    def _c_tensor_squared(self) -> Any:
         """
         Compute the square of the c_tensor.
 
@@ -68,7 +68,7 @@ def invariant3(self) -> Any:
         """
         return self.c_tensor.det()
 
-    def pk2_tensor(self, sef: sym.Expr) -> sym.MutableDenseNDimArray:
+    def pk2_tensor(self, sef: sym.Expr) -> sym.Matrix:
         """
         Compute the pk2 tensor.
 
@@ -78,9 +78,9 @@ def pk2_tensor(self, sef: sym.Expr) -> sym.MutableDenseNDimArray:
         Returns:
             sym.Matrix: The pk2 tensor.
         """
-        return sym.MutableDenseNDimArray([[sym.diff(sef, self.c_tensor[i, j]) for j in range(3)] for i in range(3)])
+        return sym.Matrix([[sym.diff(sef, self.c_tensor[i, j]) for j in range(3)] for i in range(3)])
 
-    def cmat_tensor(self, pk2: sym.Matrix) -> sym.MutableDenseNDimArray:
+    def cmat_tensor(self, pk2: sym.Matrix) -> sym.ImmutableDenseNDimArray:
         """
         Compute the cmat tensor.
 
@@ -90,27 +90,51 @@ def cmat_tensor(self, pk2: sym.Matrix) -> sym.MutableDenseNDimArray:
         Returns:
             sym.MutableDenseNDimArray: The cmat tensor.
         """
-        return sym.MutableDenseNDimArray(
+        return sym.ImmutableDenseNDimArray(
             [
                 [[[sym.diff(pk2[i, j], self.c_tensor[k, ll]) for ll in range(3)] for k in range(3)] for j in range(3)]
                 for i in range(3)
             ]
         )
 
-    def substitute(self, symbolic_tensor: sym.Expr, numerical_c_tensor: np.ndarray) -> Any:
+    def substitute(self, symbolic_tensor: sym.Expr, numerical_c_tensor: np.ndarray, *args: dict) -> Any:
         """
         Automatically substitute numerical values from a given 3x3 numerical matrix into c_tensor.
 
         Args:
             symbolic_tensor (sym.Matrix): A symbolic tensor to substitute numerical values into.
             numerical_matrix (np.ndarray): A 3x3 numerical matrix to substitute into c_tensor.
+            args (dict): Additional substitution dictionaries.
 
         Returns:
-            sym.Matrix: The c_tensor with numerical values substituted.
+            sym.Matrix: The symbolic_tensor with numerical values substituted.
 
         Raises:
             ValueError: If numerical_tensor is not a 3x3 matrix.
         """
-        return symbolic_tensor.subs(
-            {self.c_tensor[i, j]: np.array(numerical_c_tensor)[i, j] for i in range(3) for j in range(3)}
-        )
+        if not isinstance(numerical_c_tensor, np.ndarray) or numerical_c_tensor.shape != (3, 3):
+            raise ValueError("c_tensor.shape")
+
+        # Start with substitutions for c_tensor elements
+        substitutions = {self.c_tensor[i, j]: numerical_c_tensor[i, j] for i in range(3) for j in range(3)}
+        # Merge additional substitution dictionaries from *args
+        substitutions.update(*args)
+        return symbolic_tensor.subs(substitutions)
+
+    def substitute_iterator(self, symbolic_tensor: sym.Expr, numerical_c_tensors: np.ndarray, *args: dict) -> Any:
+        """
+        Automatically substitute numerical values from a given 3x3 numerical matrix into c_tensor.
+
+        Args:
+            symbolic_tensor (sym.Matrix): A symbolic tensor to substitute numerical values into.
+            numerical_matrix (np.ndarray): A 3x3 numerical matrix to substitute into c_tensor.
+            args (dict): Additional substitution dictionaries.
+
+        Returns:
+            sym.Matrix: The symbolic_tensor with numerical values substituted.
+
+        Raises:
+            ValueError: If numerical_tensor is not a 3x3 matrix.
+        """
+        for numerical_c_tensor in numerical_c_tensors:
+            yield self.substitute(symbolic_tensor, numerical_c_tensor, *args)

tests/test_kinematics.py

@@ -2,9 +2,14 @@
 import pytest
 
 from hyper_surrogate.deformation_gradient import DeformationGradientGenerator as FGen
-from hyper_surrogate.kinematics import Kinematics as K
+from hyper_surrogate.kinematics import Kinematics
 
-SIZE = 2
+SIZE = 1
+
+
+@pytest.fixture
+def K():
+    return Kinematics()
 
 
 @pytest.fixture
@@ -18,6 +23,27 @@ def right_cauchys(def_gradients):
     return np.array([np.matmul(f, f.T) for f in def_gradients])
 
 
-def test_right_cauchys(def_gradients):
+@pytest.fixture
+def left_cauchys(def_gradients):
+    return np.array([np.matmul(f.T, f) for f in def_gradients])
+
+
+def test_right_cauchys(def_gradients, K):
     right_cauchys = K.right_cauchy_green(def_gradients)
     assert np.allclose(right_cauchys, np.array([np.matmul(f.T, f) for f in def_gradients]))
+
+
+def test_left_cauchys(def_gradients, K):
+    left_cauchys = K.left_cauchy_green(def_gradients)
+    assert np.allclose(left_cauchys, np.array([np.matmul(f, f.T) for f in def_gradients]))
+
+
+def test_strain_tensor(def_gradients, K):
+    strain_tensors = K.strain_tensor(def_gradients)
+    assert np.allclose(strain_tensors, 0.5 * (np.array([f.T @ f for f in def_gradients]) - np.eye(3)))
+
+
+# def test_stretch_tensor(def_gradients, K):
+#     stretch_tensors = K.stretch_tensor(def_gradients)
+#     logging.info(stretch_tensors)
+# assert np.allclose(stretch_tensors, np.sqrt(np.array([f.T @ f for f in def_gradients])))

tests/test_symbolic_rules.py

@@ -1,5 +1,3 @@
-import logging
-
 import numpy as np
 import pytest
 import sympy as sym
@@ -11,41 +9,59 @@
 SIZE = 2
 
 
+# numeric fixtures
 @pytest.fixture
-def def_gradients():
+def def_gradients() -> np.ndarray:
     return FGen(seed=42, size=SIZE).generate()
 
 
 @pytest.fixture
-def handler():
+def right_cauchys(def_gradients) -> np.ndarray:
+    return K.right_cauchy_green(def_gradients)
+
+
+# symbolic fixtures
+@pytest.fixture
+def handler() -> SymbolicHandler:
     return SymbolicHandler()
 
 
+# mooney_rivlin
 @pytest.fixture
-def sef():
-    symbolic_handler = SymbolicHandler()
-    return (symbolic_handler.invariant1 - 3) * sym.Symbol("C10") + (symbolic_handler.invariant2 - 3) * sym.Symbol("C01")
+def sef(handler) -> sym.Expr:
+    return (handler.invariant1 - 3) * sym.Symbol("C10") + (handler.invariant2 - 3) * sym.Symbol("C01")
 
 
-# kinematics testing
 @pytest.fixture
-def right_cauchys(def_gradients):
-    return np.array([np.matmul(f, f.T) for f in def_gradients])
+def sef_args() -> dict:
+    return {"C10": 1, "C01": 1}
+
+
+@pytest.fixture
+def pk2(handler, sef) -> sym.Matrix:
+    return handler.pk2_tensor(sef)
+
+
+@pytest.fixture
+def cmat(handler, pk2) -> sym.ImmutableDenseNDimArray:
+    return handler.cmat_tensor(pk2)
+
 
+# testing
 
-def test_symbolic_pk2_cmat(handler, sef):
+
+def test_symbolic_pk2_cmat(pk2, cmat):
     # derivative of sef in order to c_tensor
-    pk2 = handler.pk2_tensor(sef)
-    cmat = handler.cmat_tensor(pk2)
+
     # assert instance
-    assert isinstance(pk2, sym.MutableDenseNDimArray)
-    assert isinstance(cmat, sym.MutableDenseNDimArray)
+    assert isinstance(pk2, sym.Matrix)
+    assert isinstance(cmat, sym.ImmutableDenseNDimArray)
     # assert shape
     assert pk2.shape == (3, 3)
     assert cmat.shape == (3, 3, 3, 3)
 
 
-def test_symbolic_subs(handler):
+def test_symbolic_subs_in_c(handler):
     # substitute the c_tensor with a 3x3 matrix of ones
     c_tensor = handler.c_tensor
     # dummy c values
@@ -56,34 +72,29 @@ def test_symbolic_subs(handler):
     assert all(c_tensor_subs[i, j] == 1 for i in range(3) for j in range(3))
 
 
-def test_symbolic_subs_in_batch(handler, def_gradients):
-    #
-    c_tensor = handler.c_tensor
-    # get first c from def_gradients
-    c = K.right_cauchy_green(def_gradients)
-    # subs c_tensor with c values
-    c_all = [handler.substitute(c_tensor, cc) for cc in c]
-    logging.info(c_all)
-
-
-# voigt notation
-# def test_voigt_notation(handler):
-#     c_tensor = handler.c_tensor
-#     c_voigt = sym.Matrix([c_tensor[0, 0], c_tensor[1, 1], c_tensor[2, 2], c_tensor[0, 1], c_tensor[0, 2], c_tensor[1, 2]])
-
-# invariant3 = symbolic_handler.invariant3_symbolic()
-
-
-# def test_derivative(def_gradients,c_tensor_symbolic):
-#     logging.info(def_gradients)
-
-#     C_10_sym=sym.Symbol("C_10_sym")
-#     #get symbols from c_tensor_symbolic
-#     C11,C12,C13,C21,C22,C23,C31,C32,C33=c_tensor_symbolic.args
-#     I3=C11*C22*C33 - C11*C23*C32 - C12*C21*C33 + C12*C23*C31 + C13*C21*C32 - C13*C22*C31
-#     trace=C11+C22+C33
-#     I1=trace*(I3)**(-1/3)
-
-#     #Symbolic Strain Energy Function
-#     SEF=(I1-3)*C_10_sym
-#     logging.info(SEF)
+def test_symbolic_subs_in_pk2(handler, pk2, right_cauchys, sef_args):
+    # right_cauchys # (N, 3, 3)
+    # for each c_tensor in c, substitute the pk2 tensor with c values and material parameters values.
+    assert all(
+        isinstance(subs, sym.Matrix)
+        and subs.shape == (3, 3)
+        and all(isinstance(subs[i, j], sym.Expr) for i in range(3) for j in range(3))
+        for subs in handler.substitute_iterator(pk2, right_cauchys, sef_args)
+    )
+
+
+def test_symbolic_subs_in_cmat(handler, cmat, right_cauchys, sef_args):
+    # right_cauchys # (N, 3, 3)
+    # for each c_tensor in c, substitute the cmat tensor with c values and material parameters values.
+    assert all(
+        isinstance(subs, sym.ImmutableDenseNDimArray)
+        and subs.shape == (3, 3, 3, 3)
+        and all(
+            isinstance(subs[i, j, k, ll], sym.Expr)
+            for i in range(3)
+            for j in range(3)
+            for k in range(3)
+            for ll in range(3)
+        )
+        for subs in handler.substitute_iterator(cmat, right_cauchys, sef_args)
+    )