Added third-order tensor transpose

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

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
 
 
 """

src/TensorAlgebra/TensorAlgebra.jl

@@ -16,6 +16,7 @@ export (⊗₁₂³⁴)
 export (⊗₁₃²⁴)
 export (⊗₁₄²³)
 export (⊗₁²)
+export (⊙₁₂₃₄²⁴)
 export ×ᵢ⁴
 export IIsym
 export I3

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)