Added 3rd- 2nd-order tensor contraction

src/TensorAlgebra/Operations.jl

@@ -650,10 +650,38 @@ The operation follows the **index contraction pattern**, where addition is perfo
   Meta.parse("TensorValue{$D²,$D²}($str)")
 end
 
+
+"""
+    contraction_IJK_KL(A::TensorValue{D,D*D}, B::TensorValue{D})::TensorValue{D,D*D}
+
+Performs a tensor contraction between third- and second-order tensors (represented as a `D × D²` matrix in flattened index notation).
+The operation follows the **index contraction pattern**, where addition is performed for repeated indices.
+"""
+@inline @generated function contraction_IJK_KL(A::TensorValue{D,D²}, B::TensorValue{D,D}) where {D, D²}
+  @assert D*D == D² "Third-order tensor size mismatch"
+  str = ""
+  for l in 1:D
+    for j in 1:D
+      for i in 1:D
+        for k in 1:D
+          a = _flat_idx(i,j,k,D)
+          b = _flat_idx(k,l,D)
+          str *= "+A[$a]*B[$b]"
+        end
+        str *= ","
+      end
+    end
+  end
+  Meta.parse("TensorValue{$D,$D²}($str)")
+end
+
+
 Gridap.TensorValues.dot(A::TensorValue{2,2}, B::TensorValue{4,4}) = contraction_IP_PJKL(A,B)
 Gridap.TensorValues.dot(A::TensorValue{3,3}, B::TensorValue{9,9}) = contraction_IP_PJKL(A,B)
 Gridap.TensorValues.dot(A::TensorValue{2,4}, B::TensorValue{2,4}) = contraction_IJK_KLP(A,B)
 Gridap.TensorValues.dot(A::TensorValue{3,9}, B::TensorValue{3,9}) = contraction_IJK_KLP(A,B)
+Gridap.TensorValues.dot(A::TensorValue{2,4}, B::TensorValue{2,2}) = contraction_IJK_KL(A,B)
+Gridap.TensorValues.dot(A::TensorValue{3,9}, B::TensorValue{3,3}) = contraction_IJK_KL(A,B)
 Gridap.TensorValues.dot(H::TensorValue{2,4}, V::VectorValue{2}) = H ⊙₁₂₃³ V
 Gridap.TensorValues.dot(H::TensorValue{3,9}, V::VectorValue{3}) = H ⊙₁₂₃³ V
 Gridap.TensorValues.dot(V::VectorValue{2}, H::TensorValue{2,4}) = V ⊙₁¹²³ H

test/TestTensorAlgebra/TensorAlgebraTests.jl

@@ -185,6 +185,23 @@ end
   end
 
   @test V · A == reference_I_IJK(V, A)
+
+  function reference_IJK_KL(A::TensorValue{3,9}, B::TensorValue{3,3})
+    D = size(A, 1)
+    C = zeros(Float64, D, D, D)
+    for i in 1:D, j in 1:D, l in 1:D
+      s = zero(Float64)
+      for k in 1:D
+        s += A[_flat_idx(i, j, k, D)] * B[k, l]
+      end
+      C[i, j, l] = s
+    end
+    TensorValue{D,D*D}(C...)
+  end
+
+  A = TensorValue{3,9}(digits3(3)...)
+  B = TensorValue{3,3}(digits2(3)...)
+  @test A · B == reference_IJK_KL(A, B)
 end