Tensor products, contractions and their identities

The products of Tensor algebra, checked rather than asserted: the index formulas, the eight algebraic identities, the two order-4 identities, and the isotropic projection they make possible.

The operators, once:

symbol	function	index formula
⊗	otimes
⊗ˢ	sotimes	symmetrized over the last index of the left and the first of the right
⊠	otimesu
⊠ˢ	otimesul
⋅	dot	one contracted index
⊡	dcontract	two contracted indices (pair-wise convention)
⊙	qcontract	four contracted indices

using TensND
using LinearAlgebra
using SymPy
using Tensors
using Random

Random.seed!(20260804)

Random.TaskLocalRNG()

Symmetrized tensor product

     :

a = Tens(Vec{3}((i,) -> symbols("a$i", real = true)))
b = Tens(Vec{3}((i,) -> symbols("b$i", real = true)))
a ⊗ b

a ⊗ˢ b

Their difference is the antisymmetric part, which removes:

tsimplify(get_array(a ⊗ b) - get_array(a ⊗ˢ b))

The two order-4 identities

   is the identity of all order-2 tensors;    only of the symmetric ones.

𝟏 = Matrix(1.0I, 3, 3)
𝕀 = get_array(tens_Id4(Val(3), Val(Float64)))

norm(𝕀 - 𝟏 ⊠ˢ 𝟏)

0.0

On a non-symmetric argument they differ: returns the symmetric part, the tensor itself.

m = rand(3, 3)
norm((𝟏 ⊠ 𝟏) ⊡ m - m), norm((𝟏 ⊠ˢ 𝟏) ⊡ m - m), norm((𝟏 ⊠ˢ 𝟏) ⊡ m - (m + m') / 2)

(0.0, 0.6530437126839718, 0.0)

The identities of the box algebra

All five hold to machine precision on random arguments.

A, B, C, D = rand(3, 3), rand(3, 3), rand(3, 3), rand(3, 3)

ids = [
    "(a⊠b):(c⊠d) = (a·c)⊠(b·d)" => norm((A ⊠ B) ⊡ (C ⊠ D) - (A * C) ⊠ (B * D)),
    "(a⊠b):c     = a·c·ᵗb" => norm((A ⊠ B) ⊡ C - A * C * B'),
    "(a⊠ˢb):c    = (a·c·ᵗb + a·ᵗc·ᵗb)/2" => norm((A ⊠ˢ B) ⊡ C - (A * C * B' + A * C' * B') / 2),
    "(a⊗b):(c⊗d) = (b:c) a⊗d" => norm((A ⊗ B) ⊡ (C ⊗ D) - sum(B .* C) * (A ⊗ D)),
    "(a⊠b):(a⁻¹⊠b⁻¹) = 1⊠1" => norm((A ⊠ B) ⊡ (inv(A) ⊠ inv(B)) - 𝟏 ⊠ 𝟏),
]
for (name, residual) in ids
    println(rpad(name, 38), " residual = ", residual)
end

(a⊠b):(c⊠d) = (a·c)⊠(b·d)              residual = 8.874191536472157e-16
(a⊠b):c     = a·c·ᵗb                   residual = 3.3306690738754696e-16
(a⊠ˢb):c    = (a·c·ᵗb + a·ᵗc·ᵗb)/2     residual = 5.978733960281817e-16
(a⊗b):(c⊗d) = (b:c) a⊗d                residual = 1.667235859327271e-15
(a⊠b):(a⁻¹⊠b⁻¹) = 1⊠1                  residual = 4.483933624673779e-15

The one that fails

There is no termwise inverse for . Equality requires and to be proportional — commuting is not enough, and taking   does not help either.

chk(x, y) = norm((x ⊠ˢ y) ⊡ (inv(x) ⊠ˢ inv(y)) - 𝟏 ⊠ˢ 𝟏)
D1, D2 = diagm(rand(3)), diagm(rand(3))     # diagonal ⟹ they commute

for (name, r) in [
        "b = a          " => chk(A, A),
        "b = 3a         " => chk(A, 3A),
        "a, b commuting " => chk(D1, D2),
        "b = 1          " => chk(A, 𝟏),
    ]
    println(name, " residual = ", round(r, sigdigits = 4))
end

b = a           residual = 7.499e-15
b = 3a          residual = 6.066e-15
a, b commuting  residual = 9.694
b = 1           residual = 7.128

Only proportionality works. This is why inversion is implemented per symmetry class rather than by one generic formula.

Quadruple contraction is the Frobenius product

and are complementary orthogonal projectors, and their norms are the dimensions of the subspaces they project onto.

𝕁 = tens_J4(Val(3), Val(Float64))
𝕂 = tens_K4(Val(3), Val(Float64))

(𝕁 ⊙ 𝕁, 𝕂 ⊙ 𝕂, 𝕁 ⊙ 𝕂)

(1.0, 5.0, 0.0)

   (one spherical direction),    (five deviatoric ones), and they are orthogonal. In dimension the second is   :

[(d, tens_K4(Val(d), Val(Float64)) ⊙ tens_K4(Val(d), Val(Float64)), d * (d + 1) ÷ 2 - 1) for d in 2:3]

2-element Vector{Tuple{Int64, Float64, Int64}}:
 (2, 2.0, 2)
 (3, 5.0, 5)

Isotropic projection

The closest isotropic tensor for the Frobenius distance, obtained by dividing each scalar product by the corresponding norm:

On a genuinely isotropic input it is exact — here a stiffness with bulk modulus and shear modulus :

k, μ = symbols("k μ", positive = true)
𝕀s, 𝕁s, 𝕂s = iso_projectors(Val(3), Val(Sym))
ℂ = 3k * 𝕁s + 2μ * 𝕂s

(tsimplify(ℂ ⊙ 𝕁s / 3), tsimplify(ℂ ⊙ 𝕂s / 10))

(k, μ)

recovering and exactly.

On an anisotropic input it is a genuine approximation. Take an orthotropic stiffness and isotropize it:

t = TensOrtho(10.0, 8.0, 9.0, 3.0, 2.0, 4.0, 2.5, 3.0, 1.5, CanonicalBasis{3, Float64}())
Ct = get_array(t)
kiso = (Ct ⊙ get_array(𝕁)) / 3
μiso = (Ct ⊙ get_array(𝕂)) / 10
(kiso, μiso)

(5.0, 2.6)

Biso, d, drel = proj_tens(:ISO, Ct)
println("closest isotropic tensor : ", Biso)
println("relative distance        : ", round(drel, sigdigits = 4))

closest isotropic tensor : (15.0) 𝕁 + (5.2) 𝕂
relative distance        : 0.1852

The projection agrees with the closed form:

norm(get_array(Biso) - (3kiso * get_array(𝕁) + 2μiso * get_array(𝕂)))

1.0877919644084146e-15

Isotropization does not commute with inversion

Projecting a stiffness and projecting its compliance give different isotropic materials, because the Euclidean distance is not invariant under inversion.

Siso, _, _ = proj_tens(:ISO, get_array(inv(t)))
lhs = get_array(inv(Biso))          # ISO(ℂ) then invert
rhs = get_array(Siso)               # invert then ISO

println("‖inv(ISO(ℂ)) − ISO(inv(ℂ))‖ = ", round(norm(lhs - rhs), sigdigits = 4))
println("relative                    = ", round(norm(lhs - rhs) / norm(rhs), sigdigits = 4))

‖inv(ISO(ℂ)) − ISO(inv(ℂ))‖ = 0.04566
relative                    = 0.09506

A reported isotropic estimate must therefore always say which of the stiffness or the compliance was projected. Distances that repair this are discussed in Isotropic tensors.

---

This page was generated using Literate.jl.