From e560154167f46f31254a1fead0877a7bb8db371e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3?= Date: Tue, 16 Sep 2025 17:53:47 +0200 Subject: [PATCH 1/6] fix docstrings --- src/TensorAlgebra/TensorAlgebra.jl | 172 ++++++++++++++++++----------- 1 file changed, 107 insertions(+), 65 deletions(-) diff --git a/src/TensorAlgebra/TensorAlgebra.jl b/src/TensorAlgebra/TensorAlgebra.jl index c495533..35c5bd7 100644 --- a/src/TensorAlgebra/TensorAlgebra.jl +++ b/src/TensorAlgebra/TensorAlgebra.jl @@ -48,15 +48,9 @@ export Ellipsoid """ - sqrt(A::TensorValue{3})::TensorValue{3} + sqrt(A::TensorValue{3})::TensorValue{3} - Compute the square root of a 3x3 matrix by means of eigen decomposition. - - # Arguments - - `A::TensorValue{3}`: the matrix to calculate the square root - - # Returns - - `::TensorValue{3}`: the squared root matrix +Compute the square root of a 3x3 matrix by means of eigen decomposition. """ function sqrt(A::TensorValue{3}) λ, Q = eigen(get_array(A)) # TODO: the get_array must be removed as long as it is supported after https://github.com/gridap/Gridap.jl/pull/1157 @@ -66,27 +60,61 @@ end """ - cof(A::TensorValue)::TensorValue - - Calculate the cofactor of a matrix. - - # Arguments - - `A::TensorValue`: the matrix to calculate. + cof(A::TensorValue)::TensorValue - # Returns - - `TensorValue`: the cofactor matrix. +Calculate the cofactor of a matrix. """ function cof(A::TensorValue) return det(A)*inv(A') end +"""Return the linear index of a N-dimensional tensor""" _flat_idx(i::Int, j::Int, N::Int) = i + N*(j-1) -_flat_idx(i::Int, j::Int, k::Int, l::Int, N::Int) = _flat_idx(_flat_idx(i,j,N), _flat_idx(k,l,N), N) +_flat_idx(i::Int, j::Int, k::Int, l::Int, N::Int) = _flat_idx(_flat_idx(i,j,N), _flat_idx(k,l,N), N*N) + +"""Return the cartesian indices of an N-dimensional second-order tensor""" _full_idx2(α::Int, N::Int) = ((α-1)%N+1 ,(α-1)÷N+1) + +"""Return the cartesian indices of an N-dimensional fourth-order tensor""" _full_idx4(α::Int, β::Int, N::Int) = (_full_idx2(α,N)..., _full_idx2(β,N)...) _full_idx4(α::Int, N::Int) = _full_idx4(_full_idx2(α,N*N)...,N) + +# ===================== +# Identity matrix +# ===================== + +const I_(N) = TensorValue{N,N,Float64}(ntuple(α -> begin + i,j = _full_idx2(α,N) + i==j ? 1.0 : 0.0 +end,N*N)) + +""" + I2::TensorValue{2} +Identity matrix 2D""" +const I2 = I_(2) + +""" + I3::TensorValue{3} +Identity matrix 3D""" +const I3 = I_(3) + +""" + I4::TensorValue{4} +Identity fourth-order tensor 2D""" +const I4 = I_(4) + +""" + I9::TensorValue{9} +Identity fourth-order tensor 3D""" +const I9 = I_(9) + + +# ===================== +# Delta Kronecker +# ===================== + function _Kroneckerδδ(δδ::Function, N::Int) TensorValue{N*N,N*N,Float64}(ntuple(α -> begin i, j, k, l = _full_idx4(α,N) @@ -95,12 +123,40 @@ function _Kroneckerδδ(δδ::Function, N::Int) N*N*N*N)) end +""" + δᵢⱼδₖₗ2D::TensorValue{4} + +Delta Kronecker outer product 2D""" const δᵢⱼδₖₗ2D = _Kroneckerδδ((i,j,k,l) -> i==j && k==l, 2) + +""" + δᵢₖδⱼₗ2D::TensorValue{4} + +Delta Kronecker outer product 2D""" const δᵢₖδⱼₗ2D = _Kroneckerδδ((i,j,k,l) -> i==k && j==l, 2) + +""" + δᵢₗδⱼₖ2D::TensorValue{4} + +Delta Kronecker outer product 2D""" const δᵢₗδⱼₖ2D = _Kroneckerδδ((i,j,k,l) -> i==l && j==k, 2) +""" + δᵢⱼδₖₗ3D::TensorValue{9} + +Delta Kronecker outer product 3D""" const δᵢⱼδₖₗ3D = _Kroneckerδδ((i,j,k,l) -> i==j && k==l, 3) + +""" + δᵢₖδⱼₗ3D::TensorValue{9} + +Delta Kronecker outer product 3D""" const δᵢₖδⱼₗ3D = _Kroneckerδδ((i,j,k,l) -> i==k && j==l, 3) + +""" + δᵢₗδⱼₖ3D::TensorValue{9} + +Delta Kronecker outer product 3D""" const δᵢₗδⱼₖ3D = _Kroneckerδδ((i,j,k,l) -> i==l && j==k, 3) @@ -119,9 +175,9 @@ end """ - **`⊗₁²(A::VectorValue{D}, B::VectorValue{D})::TensorValue{D,D}`** + ⊗₁²(A::VectorValue{D}, B::VectorValue{D})::TensorValue{D,D} - Outer product of two first-order tensors (vectors), returning a second-order tensor (matrix). +Outer product of two first-order tensors (vectors), returning a second-order tensor (matrix). """ @generated function (⊗₁²)(A::VectorValue{D,Float64}, B::VectorValue{D,Float64}) where {D} str = "" @@ -135,10 +191,10 @@ end """ - **`⊗₁₃²⁴(A::TensorValue{D}, B::TensorValue{D})::TensorValue{D*D}`** + ⊗₁₃²⁴(A::TensorValue{D}, B::TensorValue{D})::TensorValue{D*D} - Outer product of two second-order tensors (matrices), returning a fourth-order tensor - represented in a `D² x D²` flattened matrix using combined indices. +Outer product of two second-order tensors (matrices), returning a fourth-order tensor +represented in a `D² x D²` flattened matrix using combined indices. """ @generated function (⊗₁₂³⁴)(A::TensorValue{D,D,Float64}, B::TensorValue{D,D,Float64}) where {D} str = "" @@ -152,10 +208,10 @@ end """ - **`⊗₁₃²⁴(A::TensorValue{D}, B::TensorValue{D})::TensorValue{D*D}`** + ⊗₁₃²⁴(A::TensorValue{D}, B::TensorValue{D})::TensorValue{D*D} - Outer product of two second-order tensors (matrices), returning a fourth-order tensor - represented in a `D² x D²` flattened matrix using combined indices. +Outer product of two second-order tensors (matrices), returning a fourth-order tensor +represented in a `D² x D²` flattened matrix using combined indices. """ @generated function (⊗₁₃²⁴)(A::TensorValue{D}, B::TensorValue{D}) where D str = "" @@ -173,10 +229,10 @@ end """ - **`⊗₁₄²³(A::TensorValue{D}, B::TensorValue{D})::TensorValue{D*D}`** + ⊗₁₄²³(A::TensorValue{D}, B::TensorValue{D})::TensorValue{D*D} - Outer product of two second-order tensors (matrices), returning a fourth-order tensor - represented in a `D² x D²` flattened matrix using combined indices. +Outer product of two second-order tensors (matrices), returning a fourth-order tensor +represented in a `D² x D²` flattened matrix using combined indices. """ @generated function (⊗₁₄²³)(A::TensorValue{D}, B::TensorValue{D}) where D str = "" @@ -194,10 +250,10 @@ end """ - **`⊗₁²³(A::VectorValue{D}, B::TensorValue{D})::TensorValue{D,D*D}`** + ⊗₁²³(A::VectorValue{D}, B::TensorValue{D})::TensorValue{D,D*D} - Outer product of a first-order and second-order tensors (vector and matrix), - returning a third-order tensor represented in a `D x D²` flattened matrix using combined indices. +Outer product of a first-order and second-order tensors (vector and matrix), +returning a third-order tensor represented in a `D x D²` flattened matrix using combined indices. """ @generated function (⊗₁²³)(V::VectorValue{D,Float64}, A::TensorValue{D,D,Float64}) where {D} str = "" @@ -211,10 +267,10 @@ end """ - **`⊗₁²³(A::TensorValue{D}, B::VectorValue{D})::TensorValue{D,D*D}`** + ⊗₁²³(A::TensorValue{D}, B::VectorValue{D})::TensorValue{D,D*D} - Outer product of a second-order and first-order tensors (matrix and vector), - returning a third-order tensor represented in a `D x D²` flattened matrix using combined indices. +Outer product of a second-order and first-order tensors (matrix and vector), +returning a third-order tensor represented in a `D x D²` flattened matrix using combined indices. """ @generated function (⊗₁₂³)(A::TensorValue{D,D,Float64}, V::VectorValue{D,Float64}) where {D} str = "" @@ -228,10 +284,10 @@ end """ - **`⊗₁²³(A::TensorValue{D}, B::TensorValue{D})::TensorValue{D,D*D}`** + ⊗₁²³(A::TensorValue{D}, B::TensorValue{D})::TensorValue{D,D*D} - Outer product of a second-order and first-order tensors (matrix and vector), - returning a third-order tensor represented in a `D x D²` flattened matrix using combined indices. +Outer product of a second-order and first-order tensors (matrix and vector), +returning a third-order tensor represented in a `D x D²` flattened matrix using combined indices. """ @generated function (⊗₁₃²)(A::TensorValue{D}, V::VectorValue{D}) where D str = "" @@ -551,18 +607,6 @@ end end -# Identity matrix -const I_(N) = TensorValue{N,N,Float64}(ntuple(α -> begin - i,j = _full_idx2(α,N) - i==j ? 1.0 : 0.0 -end,N*N)) - -const I2 = I_(2) -const I3 = I_(3) -const I4 = I_(4) -const I9 = I_(9) - - # Jacobian regularization function logreg(J; Threshold=0.01) if J >= Threshold @@ -582,11 +626,11 @@ end """ - **`contraction_IP_PJKL(A::TensorValue{D}, H::TensorValue{D*D})::TensorValue{D*D}`** + contraction_IP_PJKL(A::TensorValue{D}, H::TensorValue{D*D})::TensorValue{D*D} - Performs a tensor contraction between a second-order tensor (of size `D × D`) - and a fourth-order tensor (represented as a `D² × D²` matrix in flattened index notation). - The operation follows the **index contraction pattern**, where addition is performed for repeated indices. +Performs a tensor contraction between a second-order tensor (of size `D × D`) +and a fourth-order tensor (represented as a `D² × D²` matrix in flattened index notation). +The operation follows the **index contraction pattern**, where addition is performed for repeated indices. """ @generated function contraction_IP_PJKL(A::TensorValue{D}, H::TensorValue{D²}) where {D, D²} @assert D*D == D² "Second and Fourth-order tensors size mismatch" @@ -596,25 +640,24 @@ end for j in 1:D for i in 1:D for p in 1:D - a = _flat_idx(p,j,D) - b = _flat_idx(k,l,D) - str *= "+A[$i,$p]*H[$a,$b]" + a = _flat_idx(p,j,k,l,D) + str *= "+A[$i,$p]*H[$a]" end str *= "," end end end end - Meta.parse("TensorValue{D²,D², Float64}($str)") + Meta.parse("TensorValue{D²}($str)") end """ - **`contraction_IP_JPKL(A::TensorValue{D}, H::TensorValue{D*D})::TensorValue{D*D}`** + contraction_IP_JPKL(A::TensorValue{D}, H::TensorValue{D*D})::TensorValue{D*D} - Performs a tensor contraction between a second-order tensor (of size `D × D`) - and a fourth-order tensor (represented as a `D² × D²` matrix in flattened index notation). - The operation follows the **index contraction pattern**, where addition is performed for repeated indices. +Performs a tensor contraction between a second-order tensor (of size `D × D`) +and a fourth-order tensor (represented as a `D² × D²` matrix in flattened index notation). +The operation follows the **index contraction pattern**, where addition is performed for repeated indices. """ @generated function contraction_IP_JPKL(A::TensorValue{D}, H::TensorValue{D²}) where {D, D²} @assert D*D == D² "Second and Fourth-order tensors size mismatch" @@ -624,16 +667,15 @@ end for j in 1:D for i in 1:D for p in 1:D - a = _flat_idx(j,p,D) - b = _flat_idx(k,l,D) - str *= "+A[$i,$p]*H[$a,$b]" + a = _flat_idx(j,p,k,l,D) + str *= "+A[$i,$p]*H[$a]" end str *= "," end end end end - Meta.parse("TensorValue{D²,D², Float64}($str)") + Meta.parse("TensorValue{D²}($str)") end From 7d6973211bcfc9c6d581f41311c833e1c15d5753 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3?= Date: Tue, 16 Sep 2025 18:10:10 +0200 Subject: [PATCH 2/6] changed inner to loops --- src/TensorAlgebra/TensorAlgebra.jl | 114 ++++++++++++++---------- test/TestTensorAlgebra/TensorAlgebra.jl | 19 ++-- 2 files changed, 76 insertions(+), 57 deletions(-) diff --git a/src/TensorAlgebra/TensorAlgebra.jl b/src/TensorAlgebra/TensorAlgebra.jl index 35c5bd7..6e31ba9 100644 --- a/src/TensorAlgebra/TensorAlgebra.jl +++ b/src/TensorAlgebra/TensorAlgebra.jl @@ -536,53 +536,73 @@ function Gridap.TensorValues.outer(A::SVector, B::SVector) return get_array(VectorValue(A) ⊗ VectorValue(B)) end -function Gridap.TensorValues.inner(Ten1::TensorValue{9,9,Float64}, Ten2::TensorValue{3,3,Float64}) - TensorValue(Ten1[1] * Ten2[1] + Ten1[10] * Ten2[2] + Ten1[19] * Ten2[3] + Ten1[28] * Ten2[4] + Ten1[37] * Ten2[5] + Ten1[46] * Ten2[6] + Ten1[55] * Ten2[7] + Ten1[64] * Ten2[8] + Ten1[73] * Ten2[9], - Ten1[2] * Ten2[1] + Ten1[11] * Ten2[2] + Ten1[20] * Ten2[3] + Ten1[29] * Ten2[4] + Ten1[38] * Ten2[5] + Ten1[47] * Ten2[6] + Ten1[56] * Ten2[7] + Ten1[65] * Ten2[8] + Ten1[74] * Ten2[9], - Ten1[3] * Ten2[1] + Ten1[12] * Ten2[2] + Ten1[21] * Ten2[3] + Ten1[30] * Ten2[4] + Ten1[39] * Ten2[5] + Ten1[48] * Ten2[6] + Ten1[57] * Ten2[7] + Ten1[66] * Ten2[8] + Ten1[75] * Ten2[9], - Ten1[4] * Ten2[1] + Ten1[13] * Ten2[2] + Ten1[22] * Ten2[3] + Ten1[31] * Ten2[4] + Ten1[40] * Ten2[5] + Ten1[49] * Ten2[6] + Ten1[58] * Ten2[7] + Ten1[67] * Ten2[8] + Ten1[76] * Ten2[9], - Ten1[5] * Ten2[1] + Ten1[14] * Ten2[2] + Ten1[23] * Ten2[3] + Ten1[32] * Ten2[4] + Ten1[41] * Ten2[5] + Ten1[50] * Ten2[6] + Ten1[59] * Ten2[7] + Ten1[68] * Ten2[8] + Ten1[77] * Ten2[9], - Ten1[6] * Ten2[1] + Ten1[15] * Ten2[2] + Ten1[24] * Ten2[3] + Ten1[33] * Ten2[4] + Ten1[42] * Ten2[5] + Ten1[51] * Ten2[6] + Ten1[60] * Ten2[7] + Ten1[69] * Ten2[8] + Ten1[78] * Ten2[9], - Ten1[7] * Ten2[1] + Ten1[16] * Ten2[2] + Ten1[25] * Ten2[3] + Ten1[34] * Ten2[4] + Ten1[43] * Ten2[5] + Ten1[52] * Ten2[6] + Ten1[61] * Ten2[7] + Ten1[70] * Ten2[8] + Ten1[79] * Ten2[9], - Ten1[8] * Ten2[1] + Ten1[17] * Ten2[2] + Ten1[26] * Ten2[3] + Ten1[35] * Ten2[4] + Ten1[44] * Ten2[5] + Ten1[53] * Ten2[6] + Ten1[62] * Ten2[7] + Ten1[71] * Ten2[8] + Ten1[80] * Ten2[9], - Ten1[9] * Ten2[1] + Ten1[18] * Ten2[2] + Ten1[27] * Ten2[3] + Ten1[36] * Ten2[4] + Ten1[45] * Ten2[5] + Ten1[54] * Ten2[6] + Ten1[63] * Ten2[7] + Ten1[72] * Ten2[8] + Ten1[81] * Ten2[9]) -end - -function Gridap.TensorValues.inner(Ten1::TensorValue{3,9,Float64}, Ten2::TensorValue{3,3,Float64}) - VectorValue(Ten1[1] * Ten2[1] + Ten1[4] * Ten2[2] + Ten1[7] * Ten2[3] + Ten1[10] * Ten2[4] + Ten1[13] * Ten2[5] + Ten1[16] * Ten2[6] + Ten1[19] * Ten2[7] + Ten1[22] * Ten2[8] + Ten1[25] * Ten2[9], - Ten1[2] * Ten2[1] + Ten1[5] * Ten2[2] + Ten1[8] * Ten2[3] + Ten1[11] * Ten2[4] + Ten1[14] * Ten2[5] + Ten1[17] * Ten2[6] + Ten1[20] * Ten2[7] + Ten1[23] * Ten2[8] + Ten1[26] * Ten2[9], - Ten1[3] * Ten2[1] + Ten1[6] * Ten2[2] + Ten1[9] * Ten2[3] + Ten1[12] * Ten2[4] + Ten1[15] * Ten2[5] + Ten1[18] * Ten2[6] + Ten1[21] * Ten2[7] + Ten1[24] * Ten2[8] + Ten1[27] * Ten2[9]) -end - -function Gridap.TensorValues.inner(Ten1::TensorValue{2,4,Float64}, Ten2::TensorValue{2,2,Float64}) - VectorValue(Ten1[1] * Ten2[1] + Ten1[3] * Ten2[2] + Ten1[5] * Ten2[3] + Ten1[7] * Ten2[4], - Ten1[2] * Ten2[1] + Ten1[4] * Ten2[2] + Ten1[6] * Ten2[3] + Ten1[8] * Ten2[4]) -end - -function Gridap.TensorValues.inner(Ten1::TensorValue{3,9,Float64}, Ten2::VectorValue{3,Float64}) - TensorValue(Ten1[1] * Ten2[1] + Ten1[10] * Ten2[2] + Ten1[19] * Ten2[3], - Ten1[2] * Ten2[1] + Ten1[11] * Ten2[2] + Ten1[20] * Ten2[3], - Ten1[3] * Ten2[1] + Ten1[12] * Ten2[2] + Ten1[21] * Ten2[3], - Ten1[4] * Ten2[1] + Ten1[13] * Ten2[2] + Ten1[22] * Ten2[3], - Ten1[5] * Ten2[1] + Ten1[14] * Ten2[2] + Ten1[23] * Ten2[3], - Ten1[6] * Ten2[1] + Ten1[15] * Ten2[2] + Ten1[24] * Ten2[3], - Ten1[7] * Ten2[1] + Ten1[16] * Ten2[2] + Ten1[25] * Ten2[3], - Ten1[8] * Ten2[1] + Ten1[17] * Ten2[2] + Ten1[26] * Ten2[3], - Ten1[9] * Ten2[1] + Ten1[18] * Ten2[2] + Ten1[27] * Ten2[3]) -end - -function Gridap.TensorValues.inner(Ten1::TensorValue{2,4,Float64}, Ten2::VectorValue{2,Float64}) - TensorValue(Ten1[1] * Ten2[1] + Ten1[5] * Ten2[2], - Ten1[2] * Ten2[1] + Ten1[6] * Ten2[2], - Ten1[3] * Ten2[1] + Ten1[7] * Ten2[2], - Ten1[4] * Ten2[1] + Ten1[8] * Ten2[2]) -end - -function Gridap.TensorValues.inner(Ten1::TensorValue{4,4,Float64}, Ten2::TensorValue{2,2,Float64}) - TensorValue(Ten1[1] * Ten2[1] + Ten1[5] * Ten2[2] + Ten1[9] * Ten2[3] + Ten1[13] * Ten2[4], - Ten1[2] * Ten2[1] + Ten1[6] * Ten2[2] + Ten1[10] * Ten2[3] + Ten1[14] * Ten2[4], - Ten1[3] * Ten2[1] + Ten1[7] * Ten2[2] + Ten1[11] * Ten2[3] + Ten1[15] * Ten2[4], - Ten1[4] * Ten2[1] + Ten1[8] * Ten2[2] + Ten1[12] * Ten2[3] + Ten1[16] * Ten2[4]) + +function str_inner_fourth_second(D) + str = "" + for j in 1:D + for i in 1:D + for l in 1:D + for k in 1:D + a = _flat_idx(i,j,k,l,D) + str *= "+H[$a]*A[$k,$l]" + end + end + str *= "," + end + end + "TensorValue{$D}($str)" +end + +function str_inner_third_second(D) + str = "" + for i in 1:D + for k in 1:D + for j in 1:D + a = _flat_idx(j,k,D) + str *= "+H[$i,$a]*A[$j,$k]" + end + end + str *= "," + end + "VectorValue{$D}($str)" +end + +function str_inner_third_first(D) + str = "" + for j in 1:D + for i in 1:D + for k in 1:D + a = _flat_idx(j,k,D) + str *= "+H[$i,$a]*V[$k]" + end + str *= "," + end + end + "TensorValue{$D}($str)" +end + +@generated function Gridap.TensorValues.inner(H::TensorValue{4,4}, A::TensorValue{2,2}) + Meta.parse(str_inner_fourth_second(2)) +end + +@generated function Gridap.TensorValues.inner(H::TensorValue{9,9}, A::TensorValue{3,3}) + Meta.parse(str_inner_fourth_second(3)) +end + +@generated function Gridap.TensorValues.inner(H::TensorValue{2,4}, A::TensorValue{2,2}) + Meta.parse(str_inner_third_second(2)) +end + +@generated function Gridap.TensorValues.inner(H::TensorValue{3,9}, A::TensorValue{3,3}) + Meta.parse(str_inner_third_second(3)) +end + +@generated function Gridap.TensorValues.inner(H::TensorValue{2,4}, V::VectorValue{2}) + Meta.parse(str_inner_third_first(2)) +end + +@generated function Gridap.TensorValues.inner(H::TensorValue{3,9}, V::VectorValue{3}) + Meta.parse(str_inner_third_first(3)) end function Gridap.TensorValues.inner(Vec1::VectorValue, Ten1::TensorValue) diff --git a/test/TestTensorAlgebra/TensorAlgebra.jl b/test/TestTensorAlgebra/TensorAlgebra.jl index d557fec..fe3e305 100644 --- a/test/TestTensorAlgebra/TensorAlgebra.jl +++ b/test/TestTensorAlgebra/TensorAlgebra.jl @@ -67,19 +67,18 @@ end end + @testset "inner" begin - A = TensorValue(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0) * 1e-3 - B = TensorValue(4.6, 2.1, 1.7, 3.2, 6.5, 1.4, 9.2, 8.0, 9.0) * 1e-3 - C = A ⊗ B - D = TensorValue([4.6 2.1 1.7 3.2 6.5 1.4 9.2 8.0 9.0; - 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0; - 5.3 2.0 3.1 1.9 5.4 9.8 0.4 8.8 3.1] * 1e-3) - E = VectorValue(1.0, 2.0, 3.0) * 1e-3 - @test norm(C ⊙ A) == 4.676298215469156e-6 - @test norm(D ⊙ E) == 0.00010313946868197451 - @test norm(D ⊙ A) == 0.0004509607632599537 + H = TensorValue(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0) + G = TensorValue{2,4}(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0) + A = TensorValue(1.0, 2.0, 3.0, 4.0) + V = VectorValue(5.0, 6.0) + @test inner(H,A) == TensorValue(90.0, 100.0, 110.0, 120.0) + @test inner(G,A) == VectorValue(50.0, 60.0) + @test inner(G,V) == TensorValue(35.0, 46.0, 57.0, 68.0) end + @testset "sum" begin A = TensorValue(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0) * 1e-3 B = TensorValue(4.6, 2.1, 1.7, 3.2, 6.5, 1.4, 9.2, 8.0, 9.0) * 1e-3 From 0f708bd0055c57c9032cc698229c3d5a848bb6c6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3?= Date: Tue, 16 Sep 2025 18:17:38 +0200 Subject: [PATCH 3/6] added more robust tests --- test/TestTensorAlgebra/TensorAlgebra.jl | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/test/TestTensorAlgebra/TensorAlgebra.jl b/test/TestTensorAlgebra/TensorAlgebra.jl index fe3e305..649dbf6 100644 --- a/test/TestTensorAlgebra/TensorAlgebra.jl +++ b/test/TestTensorAlgebra/TensorAlgebra.jl @@ -73,6 +73,10 @@ end G = TensorValue{2,4}(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0) A = TensorValue(1.0, 2.0, 3.0, 4.0) V = VectorValue(5.0, 6.0) + @test inner(H,H) == 1496.0 + @test inner(G,G) == 204.0 + @test inner(A,A) == 30.0 + @test inner(V,V) == 61.0 @test inner(H,A) == TensorValue(90.0, 100.0, 110.0, 120.0) @test inner(G,A) == VectorValue(50.0, 60.0) @test inner(G,V) == TensorValue(35.0, 46.0, 57.0, 68.0) @@ -80,9 +84,10 @@ end @testset "sum" begin - A = TensorValue(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0) * 1e-3 - B = TensorValue(4.6, 2.1, 1.7, 3.2, 6.5, 1.4, 9.2, 8.0, 9.0) * 1e-3 - @test norm(A + B) == 0.03393449572337859 + A = TensorValue(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0) + B = TensorValue(4.1, 5.2, 6.3, 7.4, 8.5, 9.6, 1.7, 2.8, 3.9) + @test A + B == TensorValue(5.1, 7.2, 9.3, 11.4, 13.5, 15.6, 8.7, 10.8, 12.9) + @test norm(A + B) ≈ 32.842807431765024 end From abb0e38c258f65d13a415488547ca60567b455cd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3?= Date: Tue, 16 Sep 2025 18:19:18 +0200 Subject: [PATCH 4/6] minor --- test/TestTensorAlgebra/TensorAlgebra.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/test/TestTensorAlgebra/TensorAlgebra.jl b/test/TestTensorAlgebra/TensorAlgebra.jl index 649dbf6..615ceb4 100644 --- a/test/TestTensorAlgebra/TensorAlgebra.jl +++ b/test/TestTensorAlgebra/TensorAlgebra.jl @@ -4,7 +4,6 @@ using HyperFEM.TensorAlgebra using Test using BenchmarkTools - @testset "Jacobian regularization" begin ∇u = TensorValue(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0) * 1e-3 From 004a26f16144142902bd498106b646b5580356db Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3?= Date: Wed, 17 Sep 2025 10:24:18 +0200 Subject: [PATCH 5/6] minor docstrings --- src/TensorAlgebra/TensorAlgebra.jl | 21 +++++++++++++++++---- 1 file changed, 17 insertions(+), 4 deletions(-) diff --git a/src/TensorAlgebra/TensorAlgebra.jl b/src/TensorAlgebra/TensorAlgebra.jl index 6e31ba9..f562e81 100644 --- a/src/TensorAlgebra/TensorAlgebra.jl +++ b/src/TensorAlgebra/TensorAlgebra.jl @@ -85,6 +85,7 @@ _full_idx4(α::Int, N::Int) = _full_idx4(_full_idx2(α,N*N)...,N) # Identity matrix # ===================== +"""The scaling N×N matrix""" const I_(N) = TensorValue{N,N,Float64}(ntuple(α -> begin i,j = _full_idx2(α,N) i==j ? 1.0 : 0.0 @@ -92,22 +93,30 @@ end,N*N)) """ I2::TensorValue{2} -Identity matrix 2D""" + +Identity matrix 2D +""" const I2 = I_(2) """ I3::TensorValue{3} -Identity matrix 3D""" + +Identity matrix 3D +""" const I3 = I_(3) """ I4::TensorValue{4} -Identity fourth-order tensor 2D""" + +Identity fourth-order tensor 2D +""" const I4 = I_(4) """ I9::TensorValue{9} -Identity fourth-order tensor 3D""" + +Identity fourth-order tensor 3D +""" const I9 = I_(9) @@ -115,6 +124,10 @@ const I9 = I_(9) # Delta Kronecker # ===================== +""" + _Kroneckerδδ(δδ::Function, N::Int)::TensorValue{N*N,N*N,Float64} + +Delta Kronecker outer product according to the `δδ(i,j,k,l)` function""" function _Kroneckerδδ(δδ::Function, N::Int) TensorValue{N*N,N*N,Float64}(ntuple(α -> begin i, j, k, l = _full_idx4(α,N) From c29ee9b1fed0021d67ce019ce885c1e319646fe9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3?= Date: Wed, 17 Sep 2025 10:31:57 +0200 Subject: [PATCH 6/6] minor update tests --- test/TestTensorAlgebra/TensorAlgebra.jl | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/test/TestTensorAlgebra/TensorAlgebra.jl b/test/TestTensorAlgebra/TensorAlgebra.jl index 615ceb4..ea3e8a5 100644 --- a/test/TestTensorAlgebra/TensorAlgebra.jl +++ b/test/TestTensorAlgebra/TensorAlgebra.jl @@ -19,15 +19,15 @@ end B = TensorValue(5.0, 6.0, 7.0, 8.0) u = VectorValue(1.0, 2.0) v = VectorValue(3.0, 4.0) - @test u ⊗ v == TensorValue(3.0, 6.0, 4.0, 8.0) - @test u ⊗₁² v == TensorValue(3.0, 6.0, 4.0, 8.0) + @test u ⊗ v == TensorValue(3.0, 6.0, 4.0, 8.0) + @test u ⊗₁² v == TensorValue(3.0, 6.0, 4.0, 8.0) @test A ⊗ B == TensorValue(5.0, 10.0, 15.0, 20.0, 6.0, 12.0, 18.0, 24.0, 7.0, 14.0, 21.0, 28.0, 8.0, 16.0, 24.0, 32.0) @test A ⊗₁₂³⁴ B == TensorValue(5.0, 10.0, 15.0, 20.0, 6.0, 12.0, 18.0, 24.0, 7.0, 14.0, 21.0, 28.0, 8.0, 16.0, 24.0, 32.0) @test A ⊗₁₃²⁴ B == TensorValue(5.0, 10.0, 6.0, 12.0, 15.0, 20.0, 18.0, 24.0, 7.0, 14.0, 8.0, 16.0, 21.0, 28.0, 24.0, 32.0) @test A ⊗₁₄²³ B == TensorValue(5.0, 10.0, 6.0, 12.0, 7.0, 14.0, 8.0, 16.0, 15.0, 20.0, 18.0, 24.0, 21.0, 28.0, 24.0, 32.0) - @test u ⊗₁²³ A == TensorValue{2,4}(1.0, 2.0, 2.0, 4.0, 3.0, 6.0, 4.0, 8.0) - @test A ⊗₁₂³ u == TensorValue{2,4}(1.0, 2.0, 3.0, 4.0, 2.0, 4.0, 6.0, 8.0) - @test A ⊗₁₃² u == TensorValue{2,4}(1.0, 2.0, 2.0, 4.0, 3.0, 4.0, 6.0, 8.0) + @test u ⊗₁²³ A == TensorValue{2,4}(1.0, 2.0, 2.0, 4.0, 3.0, 6.0, 4.0, 8.0) + @test A ⊗₁₂³ u == TensorValue{2,4}(1.0, 2.0, 3.0, 4.0, 2.0, 4.0, 6.0, 8.0) + @test A ⊗₁₃² u == TensorValue{2,4}(1.0, 2.0, 2.0, 4.0, 3.0, 4.0, 6.0, 8.0) end @@ -79,6 +79,9 @@ end @test inner(H,A) == TensorValue(90.0, 100.0, 110.0, 120.0) @test inner(G,A) == VectorValue(50.0, 60.0) @test inner(G,V) == TensorValue(35.0, 46.0, 57.0, 68.0) + @test H ⊙ A == TensorValue(90.0, 100.0, 110.0, 120.0) + @test G ⊙ A == VectorValue(50.0, 60.0) + @test G ⊙ V == TensorValue(35.0, 46.0, 57.0, 68.0) end