diff --git a/src/hyperelastic/__about__.py b/src/hyperelastic/__about__.py index 6a9beea..3d18726 100644 --- 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 index 7a11b99..82d8a99 100644 --- 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 index 11fb455..467d704 100644 --- 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)) diff --git a/tests/test_sympy.py b/tests/test_sympy.py index 135d649..06cc5d7 100644 --- a/tests/test_sympy.py +++ b/tests/test_sympy.py @@ -1,5 +1,4 @@ import felupe as fem -import numpy as np import sympy as sym import hyperelastic as hel