adtzlr · adtzlr · Jul 9, 2023 · Jul 9, 2023 · Jul 9, 2023 · Jul 9, 2023
diff --git a/src/hyperelastic/__about__.py b/src/hyperelastic/__about__.py
@@ -1 +1 @@
-__version__ = "0.4.0"
+__version__ = "0.5.0"
diff --git a/src/hyperelastic/math/_voigt.py b/src/hyperelastic/math/_voigt.py
@@ -221,40 +221,157 @@ def astensor(A, mode=2):
         return A[i.reshape(3, 3, 3, 3), j.reshape(3, 3, 3, 3)]
 
 
-def tril_from_triu(A, dim=6, inplace=True):
+def tril_from_triu(A, dim=6, out=None):
     "Copy upper triangle values to lower triangle values of a nxn tensor inplace."
-    B = A
-    if not inplace:
-        B = A.copy()
+
+    if out is None:
+        out = A.copy()
 
     i, j = np.triu_indices(dim, 1)
-    B[j, i] = A[i, j]
-    return B
+    out[j, i] = A[i, j]
+    return out
 
 
-def triu_from_tril(A, dim=6, inplace=True):
+def triu_from_tril(A, dim=6, out=None):
     "Copy lower triangle values to upper triangle values of a nxn tensor inplace."
-    B = A
-    if not inplace:
-        B = A.copy()
+
+    if out is None:
+        out = A.copy()
 
     i, j = np.tril_indices(dim, -1)
-    B[j, i] = A[i, j]
-    return B
+    out[j, i] = A[i, j]
+    return out
 
 
 def trace(A):
-    "Return the sum of the diagonal values of a 3x3 tensor."
+    r"""The trace of a symmetric second-order tensor in reduced vector storage as the
+    sum of the main diagonal.
+
+    Parameters
+    ----------
+    A : np.ndarray
+        Symmetric second-order tensor in reduced vector storage.
+
+    Returns
+    -------
+    np.ndarray
+        Trace of a symmetric second-order tensor in reduced vector storage.
+
+    Notes
+    -----
+
+    ..  math::
+
+        \text{tr}\left(\boldsymbol{C}\right) &= \boldsymbol{C} : \boldsymbol{I}
+            = C_{11} + C_{22} + C_{33}
+
+        C_{kk} &= C_{ij} : \delta_{ij}
+
+    Examples
+    --------
+    >>> from hyperelastic.math import asvoigt, trace
+    >>> import numpy as np
+
+    >>> A = np.array([1.0, 1.3, 1.5, 1.3, 1.1, 1.4, 1.5, 1.4, 1.2]).reshape(3, 3)
+
+    >>> trA = trace(asvoigt(A))
+    >>> trA
+    3.3
+
+    >>> np.allclose(trA, np.trace(A))
+    True
+
+    """
+
     return np.sum(A[:3], axis=0)
 
 
 def transpose(A):
-    "Return the major transpose of a fourth order tensor in Voigt-storage."
-    return np.einsum("ij...->ji...", A)
+    r"""The major-transpose of a symmetric fourth-order tensor in reduced vector
+    storage.
+
+    Parameters
+    ----------
+    A : np.ndarray
+        Symmetric fourth-order tensor in reduced vector storage.
+
+    Returns
+    -------
+    np.ndarray
+        Major-transpose of a symmetric fourth-order tensor in reduced vector storage.
+
+    Notes
+    -----
+
+    ..  math::
+
+        \mathbb{A}_{ijkl}^T = \mathbb{A}_{klij}
+
+
+    Examples
+    --------
+    >>> from hyperelastic.math import asvoigt, dya, transpose
+    >>> import numpy as np
+
+    >>> A = np.array([1.0, 1.3, 1.5, 1.3, 1.1, 1.4, 1.5, 1.4, 1.2]).reshape(3, 3)
+    >>> B = np.array([1.0, 1.3, 1.5, 1.3, 1.1, 1.4, 1.5, 1.4, 1.2])[::-1].reshape(3, 3)
+
+    >>> C = dya(asvoigt(A), asvoigt(B))
+    >>> C
+    array([[1.2 , 1.1 , 1.  , 1.4 , 1.3 , 1.5 ],
+           [1.32, 1.21, 1.1 , 1.54, 1.43, 1.65],
+           [1.44, 1.32, 1.2 , 1.68, 1.56, 1.8 ],
+           [1.56, 1.43, 1.3 , 1.82, 1.69, 1.95],
+           [1.68, 1.54, 1.4 , 1.96, 1.82, 2.1 ],
+           [1.8 , 1.65, 1.5 , 2.1 , 1.95, 2.25]])
+
+    >>> transpose(C)
+    array([[1.2 , 1.32, 1.44, 1.56, 1.68, 1.8 ],
+           [1.1 , 1.21, 1.32, 1.43, 1.54, 1.65],
+           [1.  , 1.1 , 1.2 , 1.3 , 1.4 , 1.5 ],
+           [1.4 , 1.54, 1.68, 1.82, 1.96, 2.1 ],
+           [1.3 , 1.43, 1.56, 1.69, 1.82, 1.95],
+           [1.5 , 1.65, 1.8 , 1.95, 2.1 , 2.25]])
+
+    >>> D = astensor(C, mode=4)
+    >>> E = astensor(transpose(C), mode=4)
+    >>> np.allclose(D, np.einsum("klij", E))
+    True
+
+    """
+
+    return np.swapaxes(A, 0, 1)
 
 
 def det(A):
-    "The determinant of a symmetric 3x3 tensor in Voigt-storage (rule of Sarrus)."
+    """The determinant of a symmetric second-order tensor in reduced vector storage.
+
+    Parameters
+    ----------
+    A : np.ndarray
+        Symmetric second-order tensor in reduced vector storage.
+
+    Returns
+    -------
+    np.ndarray
+        Determinant of a symmetric second-order tensor in reduced vector storage.
+
+    Examples
+    --------
+    >>> from hyperelastic.math import asvoigt, det
+    >>> import numpy as np
+
+    >>> A = np.array([1.0, 1.3, 1.5, 1.3, 1.1, 1.4, 1.5, 1.4, 1.2]).reshape(3, 3)
+    >>> B = asvoigt(A)
+
+    >>> det(B)
+    0.3169999999999993
+
+    >>> np.allclose(det(B), np.linalg.det(A))
+    True
+
+    """
+
     return (
         A[0] * A[1] * A[2]
         + 2 * A[3] * A[4] * A[5]
@@ -264,9 +381,48 @@ def det(A):
     )
 
 
-def inv(A, determinant=None):
-    """The inverse of a symmetric 3x3 tensor in Voigt-storage with optional provided
-    determinant."""
+def inv(A, determinant=None, out=None):
+    r"""The inverse of a symmetric second-order tensor in reduced vector storage.
+
+    Parameters
+    ----------
+    A : np.ndarray
+        Symmetric second-order tensor in reduced vector storage.
+    determinant : np.ndarray or None, optional
+        The determinant of the symmetric second-order tensor (default is None).
+    out : np.ndarray or None, optional
+        A location into which the result is stored. If provided, it must have a shape
+        that the inputs broadcast to. If not provided or None, a freshly-allocated array
+        is returned (default is None).
+
+    Returns
+    -------
+    np.ndarray
+        Inverse of a symmetric second-order tensor in reduced vector storage.
+
+    Notes
+    -----
+
+    ..  math::
+
+        \boldsymbol{C} \boldsymbol{C}^{-1} = \boldsymbol{I}
+
+
+    Examples
+    --------
+    >>> from hyperelastic.math import asvoigt, inv
+    >>> import numpy as np
+
+    >>> A = np.array([1.0, 1.3, 1.5, 1.3, 1.1, 1.4, 1.5, 1.4, 1.2]).reshape(3, 3)
+    >>> B = asvoigt(A)
+
+    >>> inv(B)
+    array([-2.01892744, -3.31230284, -1.86119874,  1.70347003,  1.73501577, 0.5362776 ])
+
+    >>> np.allclose(dot(B, inv(B)), np.eye(3))
+    True
+
+    """
 
     detAinvA = np.zeros_like(A)
 
@@ -275,15 +431,20 @@ def inv(A, determinant=None):
     else:
         detA = determinant
 
-    detAinvA[0] = A[1] * A[2] - A[4] ** 2
-    detAinvA[1] = A[0] * A[2] - A[5] ** 2
-    detAinvA[2] = A[0] * A[1] - A[3] ** 2
+    Ai = A[[1, 0, 0, 4, 3, 3]]
+    Aj = A[[2, 2, 1, 5, 5, 4]]
+    Ak = A[[4, 5, 3, 3, 0, 1]]
+    Al = A[[4, 5, 3, 2, 4, 5]]
+
+    if out is None:
+        out = Ai
+
+    Aij = np.multiply(Ai, Aj, out=Ai)
+    Akl = np.multiply(Ak, Al, out=Ak)
 
-    detAinvA[3] = A[4] * A[5] - A[3] * A[2]
-    detAinvA[4] = A[3] * A[5] - A[0] * A[4]
-    detAinvA[5] = A[3] * A[4] - A[1] * A[5]
+    detAinvA = np.subtract(Aij, Akl, out=out)
 
-    return detAinvA / detA
+    return np.divide(*np.broadcast_arrays(detAinvA, detA), out=out)
 
 
 def dot(A, B, mode=(2, 2)):
@@ -482,7 +643,7 @@ def ddot(A, B, mode=(2, 2)):
         return np.einsum("i...,ij...,i->j...", A, B, weights)
 
 
-def dya(A, B):
+def dya(A, B, out=None):
     r"""The dyadic product of two symmetric second-order tensors in reduced vector
     storage.
 
@@ -492,6 +653,10 @@ def dya(A, B):
         First symmetric second-order tensor in reduced vector storage.
     B : np.ndarray
         Second symmetric second-order tensor in reduced vector storage.
+    out : np.ndarray or None, optional
+        A location into which the result is stored. If provided, it must have a shape
+        that the inputs broadcast to. If not provided or None, a freshly-allocated array
+        is returned (default is None).
 
     Returns
     -------
@@ -528,7 +693,12 @@ def dya(A, B):
     True
 
     """
-    return A.reshape(-1, 1, *A.shape[1:]) * B.reshape(1, -1, *B.shape[1:])
+
+    return np.multiply(
+        A.reshape(-1, 1, *A.shape[1:]),
+        B.reshape(1, -1, *B.shape[1:]),
+        out=out,
+    )
 
 
 def dev(A):
@@ -555,7 +725,7 @@ def dev(A):
     return A - trace(A) / 3 * eye(A)
 
 
-def cdya_ik(A, B):
+def cdya_ik(A, B, out=None):
     r"""The overlined-dyadic product of two symmetric second-order
     tensors in reduced vector storage, where the inner indices (the second index of the
     first tensor and the first index of the second tensor) are interchanged.
@@ -566,6 +736,10 @@ def cdya_ik(A, B):
         First symmetric second-order tensor in reduced vector storage.
     B : np.ndarray
         Second symmetric second-order tensor in reduced vector storage.
+    out : np.ndarray or None, optional
+        A location into which the result is stored. If provided, it must have a shape
+        that the inputs broadcast to. If not provided or None, a freshly-allocated array
+        is returned (default is None).
 
     Returns
     -------
@@ -602,10 +776,10 @@ def cdya_ik(A, B):
 
     """
 
-    return np.einsum("ijkl...->ikjl...", astensor(dya(A, B), mode=4))
+    return np.swapaxes(astensor(dya(A, B, out=out), mode=4), 1, 2)
 
 
-def cdya_il(A, B):
+def cdya_il(A, B, out=None):
     r"""The underlined-dyadic product of two symmetric second-order
     tensors in reduced vector storage, where the right indices of the two tensors are
     interchanged.
@@ -616,6 +790,10 @@ def cdya_il(A, B):
         First symmetric second-order tensor in reduced vector storage.
     B : np.ndarray
         Second symmetric second-order tensor in reduced vector storage.
+    out : np.ndarray or None, optional
+        A location into which the result is stored. If provided, it must have a shape
+        that the inputs broadcast to. If not provided or None, a freshly-allocated array
+        is returned (default is None).
 
     Returns
     -------
@@ -652,10 +830,10 @@ def cdya_il(A, B):
 
     """
 
-    return np.einsum("ijkl...->ilkj...", astensor(dya(A, B), mode=4))
+    return np.swapaxes(astensor(dya(A, B, out=out), mode=4), 1, 3)
 
 
-def cdya(A, B):
+def cdya(A, B, out=None):
     r"""The full-symmetric crossed-dyadic product of two symmetric second-order tensors
     in reduced vector storage.
 
@@ -665,6 +843,10 @@ def cdya(A, B):
         First symmetric second-order tensor in reduced vector storage.
     B : np.ndarray
         Second symmetric second-order tensor in reduced vector storage.
+    out : np.ndarray or None, optional
+        A location into which the result is stored. If provided, it must have a shape
+        that the inputs broadcast to. If not provided or None, a freshly-allocated array
+        is returned (default is None).
 
     Returns
     -------
@@ -740,7 +922,8 @@ def cdya(A, B):
     il = b[i, l]
     kj = b[k, j]
 
-    C = np.zeros((6, *np.broadcast_shapes(A.shape, B.shape)))
+    if out is None:
+        out = np.zeros((6, *np.broadcast_shapes(A.shape, B.shape)))
 
     A, B = np.broadcast_arrays(A, B)
     Aik, Bjl, Ail, Bkj = A[ik], B[jl], A[il], B[kj]
@@ -755,9 +938,9 @@ def cdya(A, B):
     else:
         values /= 2
 
-    C[ij, kl] = C[kl, ij] = values
+    out[ij, kl] = out[kl, ij] = values
 
-    return C
+    return out
 
 
 def eigh(A, fun=None):

diff --git a/tests/test_math.py b/tests/test_math.py
@@ -1,6 +1,5 @@
 import numpy as np
 
-import hyperelastic
 import hyperelastic.math as hm
 
 
@@ -14,8 +13,9 @@ def test_math():
     D6 = hm.asvoigt(D, mode=2)
     A6 = hm.asvoigt(A, mode=4)
 
-    assert np.allclose(hm.tril_from_triu(A6, inplace=False), A6)
-    assert np.allclose(hm.triu_from_tril(A6, inplace=False), A6)
+    assert np.allclose(hm.tril_from_triu(A6), A6)
+    assert np.allclose(hm.triu_from_tril(A6), A6)
+    assert np.allclose(hm.trace(hm.dev(C6)), 0)
 
     assert np.allclose(hm.inv(C6), hm.inv(C6, determinant=hm.det(C6)))
     assert np.allclose(hm.dot(C6, C6), np.einsum("ik...,kj...->ij...", C, C))