diff --git a/src/Exports.jl b/src/Exports.jl index 6670e98..895a6e3 100644 --- a/src/Exports.jl +++ b/src/Exports.jl @@ -17,7 +17,7 @@ end @publish TensorAlgebra (⊗₁₃²⁴) @publish TensorAlgebra (⊗₁₂³⁴) @publish TensorAlgebra (⊗₁₂₃⁴) -@publish TensorAlgebra (⊗₁₂₃₄²⁴) +@publish TensorAlgebra (⊙₁₂₃₄²⁴) @publish TensorAlgebra logreg @publish TensorAlgebra Box @publish TensorAlgebra Ellipsoid diff --git a/src/TensorAlgebra/Operations.jl b/src/TensorAlgebra/Operations.jl index bafdf12..f0ab2e4 100644 --- a/src/TensorAlgebra/Operations.jl +++ b/src/TensorAlgebra/Operations.jl @@ -19,6 +19,35 @@ end end +""" + transpose_IJK_KIJ(A::TensorValue{D,D²})::TensorValue{D,D²} + +Transpose of a third-order tensor in the sense `x·A·y·z = y·Aᵀ·z·x` for all +`x,y,z`, i.e. the cyclic index shift + + Aᵀ[i,j,k] = A[k,i,j] + +This operation has order 3, not 2. `transpose_IJK_KIJ` applied twice +is not the identity. +""" +@inline @generated function transpose_IJK_KIJ(A::TensorValue{D,D²}) where {D, D²} + @assert D*D == D² "Third-order tensor sizes mismatch" + str = "" + for k in 1:D + for j in 1:D + for i in 1:D + a = _flat_idx(k,i,j,D) + str *= "A[$a]," + end + end + end + Meta.parse("TensorValue{$D,$D²}($str)") +end + +Base.transpose(A::TensorValue{2,4}) = transpose_IJK_KIJ(A) +Base.transpose(A::TensorValue{3,9}) = transpose_IJK_KIJ(A) + + function Gridap.TensorValues.outer(A::TensorValue{D,D}, B::TensorValue{D,D}) where {D} return (A ⊗₁₂³⁴ B) end @@ -592,7 +621,7 @@ The operation follows the **index contraction pattern**, where addition is perfo Meta.parse("TensorValue{D}($str)") end -(⊗₁₂₃₄²⁴) = contraction_IJKL_JL +(⊙₁₂₃₄²⁴) = contraction_IJKL_JL """ diff --git a/src/TensorAlgebra/TensorAlgebra.jl b/src/TensorAlgebra/TensorAlgebra.jl index c34add2..dfc7839 100644 --- a/src/TensorAlgebra/TensorAlgebra.jl +++ b/src/TensorAlgebra/TensorAlgebra.jl @@ -16,6 +16,7 @@ export (⊗₁₂³⁴) export (⊗₁₃²⁴) export (⊗₁₄²³) export (⊗₁²) +export (⊙₁₂₃₄²⁴) export ×ᵢ⁴ export IIsym export I3 diff --git a/test/TestTensorAlgebra/TensorAlgebraTests.jl b/test/TestTensorAlgebra/TensorAlgebraTests.jl index 5b101c5..7c9c95c 100644 --- a/test/TestTensorAlgebra/TensorAlgebraTests.jl +++ b/test/TestTensorAlgebra/TensorAlgebraTests.jl @@ -4,6 +4,12 @@ using HyperFEM.TensorAlgebra using Test +digits1(D) = [Float64(i) for i in 1:D] +digits2(D) = [Float64(10i + j) for i in 1:D, j in 1:D] +digits3(D) = [Float64(100i + 10j + k) for i in 1:D, j in 1:D, k in 1:D] +digits4(D) = [Float64(1000i + 100j + 10k + l) for i in 1:D, j in 1:D, k in 1:D, l in 1:D] + + @testset "Flat indexing" begin A = rand(3,3) B = rand(3,3,3) @@ -20,6 +26,20 @@ using Test end +@testset "transpose" begin + function reference_transpose_IJK_KIJ(A::TensorValue{3,9}) + D = size(A, 1) + C = zeros(Float64, D, D, D) + for i in 1:D, j in 1:D, k in 1:D + C[i, j, k] = A[_flat_idx(k, i, j, D)] + end + TensorValue{D,D*D}(C...) + end + A = TensorValue{3,9}(digits3(3)...) + @test A' == reference_transpose_IJK_KIJ(A) +end + + @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 F = one(∇u) + ∇u @@ -117,11 +137,6 @@ end @test contraction_IP_PJKL(A,H) == TensorValue(7.0, 10.0, 15.0, 22.0, 23.0, 34.0, 31.0, 46.0, 39.0, 58.0, 47.0, 70.0, 55.0, 82.0, 63.0, 94.0) @test contraction_IP_JPKL(A,H) == TensorValue(10.0, 14.0, 14.0, 20.0, 26.0, 38.0, 30.0, 44.0, 42.0, 62.0, 46.0, 68.0, 58.0, 86.0, 62.0, 92.0) - digits1(D) = [Float64(i) for i in 1:D] - digits2(D) = [Float64(10i + j) for i in 1:D, j in 1:D] - digits3(D) = [Float64(100i + 10j + k) for i in 1:D, j in 1:D, k in 1:D] - digits4(D) = [Float64(1000i + 100j + 10k + l) for i in 1:D, j in 1:D, k in 1:D, l in 1:D] - function reference_IJK_KLP(A::TensorValue{3,9}, B::TensorValue{3,9}) D = size(A, 1) C = zeros(Float64, D, D, D, D)