Tensors

A tensor, parametrized by an order and a dimension, is in general defined by

an array or a set of condensed parameters (e.g. isotropic tensors),
a basis,
a set of variances (covariant :cov or contravariant :cont) useful if the basis is not orthonormal.

In practice, the type of basis conditions the type of tensor (TensCanonical, TensRotated, TensOrthogonal, Tens or even TensISO in case of isotropic tensor).

julia> using TensND, SymPy, Tensors
julia> ℬ = Basis(Sym[0 1 1; 1 0 1; 1 1 0])3×3 Basis{3, Sym}:
 0  1  1
 1  0  1
 1  1  0
julia> V = Tens(Tensor{1,3}(i -> symbols("v$i", real = true)))3-element TensND.TensCanonical{1, 3, Sym{PyCall.PyObject}, Tensors.Vec{3, Sym{PyCall.PyObject}}}:
 v₁
 v₂
 v₃
julia> components(V, ℬ, (:cont,))3-element Vector{Sym{PyCall.PyObject}}:
 -v1/2 + v2/2 + v3/2
  v1/2 - v2/2 + v3/2
  v1/2 + v2/2 - v3/2
julia> components(V, ℬ, (:cov,))3-element Vector{Sym{PyCall.PyObject}}:
 v₂ + v₃
 v₁ + v₃
 v₁ + v₂
julia> ℬ̄ = normalize(ℬ)ERROR: UndefVarError: `normalize` not defined in `Main`
Suggestion: check for spelling errors or missing imports.
Hint: a global variable of this name also exists in LinearAlgebra.
julia> components(V, ℬ̄, (:cov,))ERROR: UndefVarError: `ℬ̄` not defined in `Main`
Suggestion: add an appropriate import or assignment. This global was declared but not assigned.
julia> T = Tens(Tensor{2,3}((i, j) -> symbols("t$i$j", real = true)))3×3 TensND.TensCanonical{2, 3, Sym{PyCall.PyObject}, Tensors.Tensor{2, 3, Sym{PyCall.PyObject}, 9}}:
 t₁₁  t₁₂  t₁₃
 t₂₁  t₂₂  t₂₃
 t₃₁  t₃₂  t₃₃
julia> components(T, ℬ, (:cov, :cov))3×3 Matrix{Sym{PyCall.PyObject}}:
 t₂₂ + t₂₃ + t₃₂ + t₃₃  t₂₁ + t₂₃ + t₃₁ + t₃₃  t₂₁ + t₂₂ + t₃₁ + t₃₂
 t₁₂ + t₁₃ + t₃₂ + t₃₃  t₁₁ + t₁₃ + t₃₁ + t₃₃  t₁₁ + t₁₂ + t₃₁ + t₃₂
 t₁₂ + t₁₃ + t₂₂ + t₂₃  t₁₁ + t₁₃ + t₂₁ + t₂₃  t₁₁ + t₁₂ + t₂₁ + t₂₂
julia> factor(simplify(components(T, ℬ, (:cont, :cov))))ERROR: MethodError: no method matching factor(::Matrix{Sym{PyCall.PyObject}})
The function `factor` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  factor()
   @ SymPyCore ~/.julia/packages/SymPyCore/l4Bbg/src/gen_methods.jl:280
  factor(::Sym, Any...; kwargs...)
   @ SymPy ~/.julia/packages/SymPyCore/l4Bbg/src/SymPy/gen_methods_sympy.jl:33

Special tensors are available

tensId2(::Val{dim} = Val(3), ::Val{T} = Val(Sym)) where {dim,T<:Number}: second-order identity (𝟏ᵢⱼ = δᵢⱼ = 1 if i=j otherwise 0)
tensId4(::Val{dim} = Val(3), ::Val{T} = Val(Sym)) where {dim,T<:Number}: fourth-order identity with minor symmetries (𝕀 = 𝟏 ⊠ˢ 𝟏 i.e. (𝕀)ᵢⱼₖₗ = (δᵢₖδⱼₗ+δᵢₗδⱼₖ)/2)
tensJ4(::Val{dim} = Val(3), ::Val{T} = Val(Sym)) where {dim,T<:Number}: fourth-order spherical projector (𝕁 = (𝟏 ⊗ 𝟏) / dim i.e. (𝕁)ᵢⱼₖₗ = δᵢⱼδₖₗ/dim)
tensK4(::Val{dim} = Val(3), ::Val{T} = Val(Sym)) where {dim,T<:Number}: fourth-order deviatoric projector (𝕂 = 𝕀 - 𝕁 i.e. (𝕂)ᵢⱼₖₗ = (δᵢₖδⱼₗ+δᵢₗδⱼₖ)/2 - δᵢⱼδₖₗ/dim)
ISO(::Val{dim} = Val(3), ::Val{T} = Val(Sym)) where {dim,T<:Number}: returns 𝕀, 𝕁, 𝕂

The useful tensor products are the following:

⊗ tensor product
⊗ˢ symmetrized tensor product
⊠ modified tensor product
⊠ˢ symmetrized modified tensor product
⋅ contracted product
⊡ double contracted product
⊙ quadruple contracted product

julia> 𝟏 = tensId2(3, Sym)3×3 TensISO{2, 3, Sym, 1}:
 1  ⋅  ⋅
 ⋅  1  ⋅
 ⋅  ⋅  1
julia> 𝕀, 𝕁, 𝕂 = ISO(3, Sym) ;
julia> 𝕀 == 𝟏 ⊠ˢ 𝟏true
julia> 𝕁 == (𝟏 ⊗ 𝟏)/3true
julia> a = Tens(Vec{3}((i,) -> symbols("a$i", real = true))) ;
julia> b = Tens(Vec{3}((i,) -> symbols("b$i", real = true))) ;
julia> a ⊗ b3×3 TensND.TensCanonical{2, 3, Sym{PyCall.PyObject}, Tensors.Tensor{2, 3, Sym{PyCall.PyObject}, 9}}:
 a₁⋅b₁  a₁⋅b₂  a₁⋅b₃
 a₂⋅b₁  a₂⋅b₂  a₂⋅b₃
 a₃⋅b₁  a₃⋅b₂  a₃⋅b₃
julia> a ⊗ˢ b3×3 TensND.TensCanonical{2, 3, Sym{PyCall.PyObject}, Tensors.SymmetricTensor{2, 3, Sym{PyCall.PyObject}, 6}}:
             a₁⋅b₁  a1*b2/2 + a2*b1/2  a1*b3/2 + a3*b1/2
 a1*b2/2 + a2*b1/2              a₂⋅b₂  a2*b3/2 + a3*b2/2
 a1*b3/2 + a3*b1/2  a2*b3/2 + a3*b2/2              a₃⋅b₃

The predefined spherical coordinate system init_spherical() provides the local orthonormal basis $(\mathbf{e}_\theta, \mathbf{e}_\varphi, \mathbf{e}_r)$ in terms of polar angle $\theta$ (from the $z$-axis) and azimuthal angle $\varphi$:

\[\mathbf{e}_\theta = \cos\theta\cos\varphi\,\mathbf{e}_1 + \cos\theta\sin\varphi\,\mathbf{e}_2 - \sin\theta\,\mathbf{e}_3\]

\[\mathbf{e}_\varphi = -\sin\varphi\,\mathbf{e}_1 + \cos\varphi\,\mathbf{e}_2\]

\[\mathbf{e}_r = \sin\theta\cos\varphi\,\mathbf{e}_1 + \sin\theta\sin\varphi\,\mathbf{e}_2 + \cos\theta\,\mathbf{e}_3\]

The rotation matrix $R = [\mathbf{e}_\theta \mid \mathbf{e}_\varphi \mid \mathbf{e}_r]$ encodes this change of basis. For any vector $\mathbf{A}$ in the canonical frame, change_tens(A, ℬˢ) returns its components in the spherical basis. The example below verifies that if $\mathbf{A} = R\,\mathbf{a}$, then expressing $\mathbf{A}$ in $\mathcal{B}^s$ recovers the original components $(a_1, a_2, a_3)$:

julia> (θ, ϕ, r), (𝐞ᶿ, 𝐞ᵠ, 𝐞ʳ), ℬˢ = init_spherical()((θ, ϕ, r), (Sym{PyCall.PyObject}[1, 0, 0], Sym{PyCall.PyObject}[0, 1, 0], Sym{PyCall.PyObject}[0, 0, 1]), Sym{PyCall.PyObject}[cos(θ)*cos(ϕ) -sin(ϕ) sin(θ)*cos(ϕ); sin(ϕ)*cos(θ) cos(ϕ) sin(θ)*sin(ϕ); -sin(θ) 0 cos(θ)])
julia> R = rot3(θ, ϕ)3×3 RotZYZ{Sym{PyObject}} with indices SOneTo(3)×SOneTo(3)(ϕ, θ, 0):
 cos(θ)⋅cos(ϕ)  -sin(ϕ)  sin(θ)⋅cos(ϕ)
 sin(ϕ)⋅cos(θ)   cos(ϕ)  sin(θ)⋅sin(ϕ)
       -sin(θ)        0         cos(θ)
julia> A = Tens(R * a)3-element TensND.TensCanonical{1, 3, Sym{PyCall.PyObject}, Tensors.Vec{3, Sym{PyCall.PyObject}}}:
 a₁⋅cos(θ)⋅cos(ϕ) - a₂⋅sin(ϕ) + a₃⋅sin(θ)⋅cos(ϕ)
 a₁⋅sin(ϕ)⋅cos(θ) + a₂⋅cos(ϕ) + a₃⋅sin(θ)⋅sin(ϕ)
                          -a₁⋅sin(θ) + a₃⋅cos(θ)
julia> simplify(change_tens(A, ℬˢ))3-element TensND.TensRotated{1, 3, Sym{PyCall.PyObject}, Tensors.Vec{3, Sym{PyCall.PyObject}}}:
 -(-a1*sin(θ) + a3*cos(θ))*sin(θ) + (a1*sin(ϕ)*cos(θ) + a2*cos(ϕ) + a3*sin(θ)*sin(ϕ))*sin(ϕ)*cos(θ) + (a1*cos(θ)*cos(ϕ) - a2*sin(ϕ) + a3*sin(θ)*cos(ϕ))*cos(θ)*cos(ϕ)
                                                  (a1*sin(ϕ)*cos(θ) + a2*cos(ϕ) + a3*sin(θ)*sin(ϕ))*cos(ϕ) - (a1*cos(θ)*cos(ϕ) - a2*sin(ϕ) + a3*sin(θ)*cos(ϕ))*sin(ϕ)
  (-a1*sin(θ) + a3*cos(θ))*cos(θ) + (a1*sin(ϕ)*cos(θ) + a2*cos(ϕ) + a3*sin(θ)*sin(ϕ))*sin(θ)*sin(ϕ) + (a1*cos(θ)*cos(ϕ) - a2*sin(ϕ) + a3*sin(θ)*cos(ϕ))*sin(θ)*cos(ϕ)

Isotropic tensors (TensISO)

Isotropic tensors are stored compactly: a second-order isotropic tensor $\lambda\mathbf{1}$ is parametrized by one scalar, while a fourth-order isotropic tensor $\alpha\mathbb{J} + \beta\mathbb{K}$ is parametrized by two scalars. All arithmetic operations ($+$, $-$, $\times$, $\mathbb{A}:\mathbb{B}$, $\mathbb{A}^{-1}$) exploit this compact form and remain in the TensISO type whenever possible.

The type predicates isISO, isTI, isOrtho allow querying the symmetry class of any tensor:

julia> using TensND, Tensors
julia> 𝟏 = tensId2(Val(3), Val(Float64))3×3 TensISO{2, 3, Float64, 1}:
 1.0   ⋅    ⋅
  ⋅   1.0   ⋅
  ⋅    ⋅   1.0
julia> 𝕀, 𝕁, 𝕂 = ISO(Val(3), Val(Float64)) ;
julia> isISO(𝕀)true
julia> isTI(𝕀)false
julia> isOrtho(𝕀)false

The compact display reflects the algebraic form directly:

julia> show(stdout, 𝕁 + 𝕂)   # prints "(1.0) 𝕁 + (1.0) 𝕂"(1.0) 𝕁 + (1.0) 𝕂
julia> show(stdout, 2.0 * 𝟏)  # prints "(2.0) 𝟏"(2.0) 𝟏

Transverse isotropy and orthotropy

TensWalpole

A transversely isotropic 4th-order tensor with symmetry axis $\mathbf{n}$ is decomposed in the Walpole basis:

\[L = \ell_1 W_1 + \ell_2 W_2 + \ell_3 W_3 + \ell_4 W_4 + \ell_5 W_5 + \ell_6 W_6\]

where $\mathbf{n}_n = \mathbf{n}\otimes\mathbf{n}$, $\mathbf{n}_T = \mathbf{1} - \mathbf{n}_n$ and

Tensor	Expression
$W_1$	$\mathbf{n}_n\otimes\mathbf{n}_n$
$W_2$	$(\mathbf{n}_T\otimes\mathbf{n}_T)/2$
$W_3$	$(\mathbf{n}_n\otimes\mathbf{n}_T)/\sqrt{2}$
$W_4$	$(\mathbf{n}_T\otimes\mathbf{n}_n)/\sqrt{2}$
$W_5$	$\mathbf{n}_T\,\overline{\boxtimes}^s\,\mathbf{n}_T - (\mathbf{n}_T\otimes\mathbf{n}_T)/2$
$W_6$	$\mathbf{n}_T\,\overline{\boxtimes}^s\,\mathbf{n}_n + \mathbf{n}_n\,\overline{\boxtimes}^s\,\mathbf{n}_T$

The double contraction follows the synthetic Walpole rule:

\[L\colon M \equiv \left(\begin{bmatrix}\ell_1 & \ell_3\\\ell_4 & \ell_2\end{bmatrix}\begin{bmatrix}m_1 & m_3\\m_4 & m_2\end{bmatrix},\; \ell_5 m_5,\; \ell_6 m_6\right)\]

For major-symmetric tensors ($\ell_3=\ell_4$), use N=5; for general tensors, N=6.

The show method displays the tensor in its compact Walpole form, including the symmetry axis:

julia> using TensND, Tensors
julia> n = 𝐞(3) ;
julia> W1, W2, W3, W4, W5, W6 = Walpole(n) ;
julia> L = TensWalpole(2., 1., 0.5, 0.3, 0.8, n)3×3×3×3 TensWalpole{Sym{PyCall.PyObject}, 5}:
[:, :, 1, 1] =
 0.650000000000000                  0             0
                 0  0.350000000000000             0
                 0                  0  0.25⋅√2

[:, :, 2, 1] =
                 0  0.150000000000000  0
 0.150000000000000                  0  0
                 0                  0  0

[:, :, 3, 1] =
                 0  0  0.400000000000000
                 0  0                  0
 0.400000000000000  0                  0

[:, :, 1, 2] =
                 0  0.150000000000000  0
 0.150000000000000                  0  0
                 0                  0  0

[:, :, 2, 2] =
 0.350000000000000                  0             0
                 0  0.650000000000000             0
                 0                  0  0.25⋅√2

[:, :, 3, 2] =
 0                  0                  0
 0                  0  0.400000000000000
 0  0.400000000000000                  0

[:, :, 1, 3] =
                 0  0  0.400000000000000
                 0  0                  0
 0.400000000000000  0                  0

[:, :, 2, 3] =
 0                  0                  0
 0                  0  0.400000000000000
 0  0.400000000000000                  0

[:, :, 3, 3] =
 0.25⋅√2             0                 0
            0  0.25⋅√2                 0
            0             0  2.00000000000000
julia> show(stdout, L)(2.00000000000000) W₁ˢ + (1.00000000000000) W₂ˢ + (0.500000000000000) W₃ˢ + (0.300000000000000) W₄ˢ + (0.800000000000000) W₅ˢ
  axis n = (0, 0, 1)
julia> maximum(abs.(getarray(L ⊡ inv(L)) - getarray(tensId4(Val(3), Val(Float64)))))0
julia> 𝕀, 𝕁, 𝕂 = ISO() ; L2 = fromISO(3𝕁 + 2𝕂, n)3×3×3×3 TensWalpole{Sym{PyCall.PyObject}, 5}:
[:, :, 1, 1] =
 7/3    0    0
   0  1/3    0
   0    0  1/3

[:, :, 2, 1] =
 0  1  0
 1  0  0
 0  0  0

[:, :, 3, 1] =
 0  0  1
 0  0  0
 1  0  0

[:, :, 1, 2] =
 0  1  0
 1  0  0
 0  0  0

[:, :, 2, 2] =
 1/3    0    0
   0  7/3    0
   0    0  1/3

[:, :, 3, 2] =
 0  0  0
 0  0  1
 0  1  0

[:, :, 1, 3] =
 0  0  1
 0  0  0
 1  0  0

[:, :, 2, 3] =
 0  0  0
 0  0  1
 0  1  0

[:, :, 3, 3] =
 1/3    0    0
   0  1/3    0
   0    0  7/3
julia> isTI(L)true
julia> isISO(L)false
julia> isOrtho(L)false

An isotropic tensor converted to TensWalpole via fromISO retains the isTI predicate, and symbolic manipulations via tsimplify, tsubs, tdiff, etc. preserve the TensWalpole type:

julia> using SymPy
julia> ℓ₁, ℓ₂, ℓ₃ = symbols("ℓ₁ ℓ₂ ℓ₃", real = true) ;
julia> ns = 𝐞(Val(3), Val(3), Val(Sym)) ;
julia> Ls = TensWalpole(ℓ₁, ℓ₂, ℓ₃, ℓ₁ + ℓ₂, ℓ₂ + ℓ₃, ns) ;
julia> Ls_simp = tsimplify(Ls) ;
julia> Ls_simp isa TensWalpoletrue

TensOrtho

An orthotropic 4th-order tensor in material frame $(\mathbf{e}_1,\mathbf{e}_2,\mathbf{e}_3)$ with $P_m = \mathbf{e}_m\otimes\mathbf{e}_m$ has 9 independent elastic constants:

\[\mathbb{C} = C_{11}P_1{\otimes}P_1 + C_{22}P_2{\otimes}P_2 + C_{33}P_3{\otimes}P_3 + C_{12}(P_1{\otimes}P_2+P_2{\otimes}P_1) + C_{13}(P_1{\otimes}P_3+P_3{\otimes}P_1) + C_{23}(P_2{\otimes}P_3+P_3{\otimes}P_2) + 2C_{44}(P_2\,\overline{\boxtimes}^s P_3) + 2C_{55}(P_1\,\overline{\boxtimes}^s P_3) + 2C_{66}(P_1\,\overline{\boxtimes}^s P_2)\]

The Kelvin-Mandel matrix in the material frame (ordering $11,22,33,23,13,12$) is block-diagonal. Use KM_material(t) to retrieve it; KM(t) gives the matrix in the canonical frame.

The show method displays all 9 constants and the material frame, and isOrtho identifies the type:

julia> using TensND, Tensors
julia> ℬ = CanonicalBasis{3,Float64}() ;
julia> t = TensOrtho(10., 8., 9., 3., 2., 4., 2.5, 3., 1.5, ℬ) ;
julia> show(stdout, t)(10.0) P₁⊗P₁ + (8.0) P₂⊗P₂ + (9.0) P₃⊗P₃
  + (3.0)(P₁⊗P₂+P₂⊗P₁) + (2.0)(P₁⊗P₃+P₃⊗P₁) + (4.0)(P₂⊗P₃+P₃⊗P₂)
  + 2(2.5)(P₂⊠ˢP₃) + 2(3.0)(P₁⊠ˢP₃) + 2(1.5)(P₁⊠ˢP₂)
  frame: [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]
julia> KM_material(t)6×6 Matrix{Float64}:
 10.0  3.0  2.0  0.0  0.0  0.0
  3.0  8.0  4.0  0.0  0.0  0.0
  2.0  4.0  9.0  0.0  0.0  0.0
  0.0  0.0  0.0  5.0  0.0  0.0
  0.0  0.0  0.0  0.0  6.0  0.0
  0.0  0.0  0.0  0.0  0.0  3.0
julia> maximum(abs.(getarray(t) ⊡ getarray(inv(t)) - getarray(tensId4(Val(3), Val(Float64)))))1.249000902703301e-16
julia> isOrtho(t)true
julia> isTI(t)false
julia> isISO(t)false

TensTI (2nd-order transversely isotropic)

A 2nd-order transversely isotropic tensor with symmetry axis $\mathbf{n}$ is decomposed as:

\[\mathbf{A} = a\,\mathbf{n}_T + b\,\mathbf{n}_n\]

where $\mathbf{n}_n = \mathbf{n}\otimes\mathbf{n}$ and $\mathbf{n}_T = \mathbf{1} - \mathbf{n}_n$. The scalar $a$ is the transverse coefficient (in the plane $\perp\mathbf{n}$) and $b$ is the axial coefficient (along $\mathbf{n}$).

The type TensTI{order,T,N} mirrors the parametric design of TensISO{order,dim,T,N}. For order 2, N=2 and data = (a, b). When a = b, the tensor is isotropic.

julia> using TensND, LinearAlgebra
julia> n = [0., 0., 1.] ;
julia> A = TensTI{2}(5.0, 8.0, n)3×3 TensTI{2, Float64, 2}:
 5.0  0.0  0.0
 0.0  5.0  0.0
 0.0  0.0  8.0
julia> getarray(A)3×3 Matrix{Float64}:
 5.0  0.0  0.0
 0.0  5.0  0.0
 0.0  0.0  8.0
julia> tr(A)18.0
julia> isISO(A)false
julia> isTI(A)true
julia> inv(A).data(0.2, 0.125)

TensTI{2} supports all standard operations: addition, subtraction, scalar multiplication, inversion. Two TI tensors with the same axis can be combined:

julia> B = TensTI{2}(3.0, 2.0, n) ;
julia> (A + B).data(8.0, 10.0)
julia> (A - B).data(2.0, 6.0)
julia> (2.0 * A).data(10.0, 16.0)

When the axis is not aligned with $\mathbf{e}_3$, the full 3×3 matrix correctly reflects the off-diagonal structure:

julia> n45 = [1/√2, 0., 1/√2] ;
julia> C = TensTI{2}(3.0, 7.0, n45) ;
julia> getarray(C)3×3 Matrix{Float64}:
 5.0  0.0  2.0
 0.0  3.0  0.0
 2.0  0.0  5.0

TI convenience constructors and engineering parametrizations

Three parametrizations are available for constructing TI 4th-order tensors (stiffness or compliance). They all return a TensWalpole{T,5} and have corresponding extraction functions.

Direct component form: `tensTI` / `argTI`

Construct from the 5 independent components $C_{1111}, C_{1122}, C_{1133}, C_{3333}, C_{2323}$ (axis $\mathbf{n} = \mathbf{e}_3$). Works for both stiffness and compliance tensors:

julia> using TensND
julia> n = [0., 0., 1.] ;
julia> C = tensTI(10., 3., 2.5, 12., 2., n) ;
julia> argTI(C)(10.0, 3.0, 2.5, 12.0, 2.0)

Engineering form: `tensTI_eng` / `argTI_eng`

Construct the TI compliance tensor from 5 engineering constants commonly used in composite mechanics:

$E_1$: transverse Young's modulus (isotropic plane)
$E_3$: axial Young's modulus (symmetry axis)
$\nu_{12}$: in-plane Poisson's ratio
$\nu_{31}$: axial-transverse Poisson's ratio ($\nu_{31}/E_3 = \nu_{13}/E_1$)
$G_{31}$: axial shear modulus

To obtain the stiffness tensor, invert the result:

julia> 𝕊 = tensTI_eng(72., 50., 0.3, 0.25, 15., n) ;
julia> argTI_eng(𝕊)(72.0, 50.0, 0.29999999999999993, 0.25, 15.0)
julia> ℂ = inv(𝕊) ;
julia> maximum(abs.(getarray(ℂ ⊡ 𝕊) - getarray(tensId4(Val(3), Val(Float64)))))7.850462293418875e-17

Hoenig form: `tensTI_Hoenig` / `argTI_Hoenig`

An alternative parametrization (Hoenig, 1978) expressed as dimensionless ratios relative to the transverse Young's modulus $E$:

$E$: transverse Young's modulus ($= 1/S_{1111}$)
$\nu_1$: in-plane Poisson's ratio ($= -E\,S_{1122}$)
$\nu_2$: axial-transverse Poisson's ratio ($= -E\,S_{1133}$)
$H$: axial-to-transverse modulus ratio ($= 1/(E\,S_{3333})$)
$\Gamma$: shear anisotropy parameter ($= (1+\nu_1)/(2\,E\,S_{2323})$)

The Hoenig parametrization is useful when discussing anisotropy ratios independently of the overall stiffness scale:

julia> 𝕊h = tensTI_Hoenig(72., 0.3, 0.25, 0.7, 0.9, n) ;
julia> argTI_Hoenig(𝕊h)(72.0, 0.29999999999999993, 0.24999999999999994, 0.7000000000000001, 0.8999999999999998)

All three parametrizations are interconvertible through the TensWalpole representation; going from one to another simply requires calling the appropriate arg* extractor on a tensor built via the corresponding constructor.

Accessing Walpole coefficients

For any TensWalpole, the function get_ℓ returns the 6 Walpole coefficients as a tuple (for N=5, $\ell_3 = \ell_4$ is repeated):

julia> using TensND
julia> n = 𝐞(3) ;
julia> L = TensWalpole(2., 1., 0.5, 0.3, 0.8, n) ;
julia> get_ℓ(L)(2.00000000000000, 1.00000000000000, 0.500000000000000, 0.500000000000000, 0.300000000000000, 0.800000000000000)
julia> getaxis(L)(0, 0, 1)

For TensOrtho, the material frame is accessible via getframe:

julia> using TensND, Tensors
julia> ℬ = CanonicalBasis{3,Float64}() ;
julia> t = TensOrtho(10., 8., 9., 3., 2., 4., 2.5, 3., 1.5, ℬ) ;
julia> getframe(t)3×3 CanonicalBasis{3, Float64}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0
julia> getdata(t)(10.0, 8.0, 9.0, 3.0, 2.0, 4.0, 2.5, 3.0, 1.5)

Kelvin-Mandel representation

The Kelvin-Mandel (KM) representation maps a symmetric tensor to a vector (order 2) or a matrix (order 4) using the index ordering $11 \to 1,\; 22 \to 2,\; 33 \to 3,\; 23 \to 4,\; 13 \to 5,\; 12 \to 6$. Off-diagonal components are scaled by $\sqrt{2}$ so that the Frobenius norm is preserved: $\lVert\mathbf{A}\rVert^2 = A_{KM}^T A_{KM}$. This differs from Voigt notation which uses engineering shear (factor 2 instead of $\sqrt{2}$).

KM(t): KM matrix/vector in the canonical frame
KM_material(t::TensOrtho): KM matrix in the material frame (block-diagonal for orthotropic)
invKM(km): reconstruct a symmetric tensor from its KM representation

julia> KM_material(t)6×6 Matrix{Float64}:
 10.0  3.0  2.0  0.0  0.0  0.0
  3.0  8.0  4.0  0.0  0.0  0.0
  2.0  4.0  9.0  0.0  0.0  0.0
  0.0  0.0  0.0  5.0  0.0  0.0
  0.0  0.0  0.0  0.0  6.0  0.0
  0.0  0.0  0.0  0.0  0.0  3.0

Symmetry class predicates

The three predicates isISO, isTI, isOrtho form a consistent hierarchy across all specialized tensor types. Any value that is not a recognized tensor type returns false for all three:

julia> using TensND, Tensors
julia> 𝕀, 𝕁, 𝕂 = ISO(Val(3), Val(Float64)) ;
julia> n = 𝐞(3) ;
julia> L = TensWalpole(2., 1., 0.5, 3., 4., n) ;
julia> A2 = TensTI{2}(5.0, 8.0, n) ;
julia> ℬ = CanonicalBasis{3,Float64}() ;
julia> t = TensOrtho(10., 8., 9., 3., 2., 4., 2.5, 3., 1.5, ℬ) ;
julia> (isISO(𝕀), isTI(𝕀),  isOrtho(𝕀))(true, false, false)
julia> (isISO(L),  isTI(L),  isOrtho(L))(false, true, false)
julia> (isISO(A2), isTI(A2), isOrtho(A2))(false, true, false)
julia> (isISO(t),  isTI(t),  isOrtho(t))(false, false, true)

Projection onto symmetry subspaces

Given an arbitrary tensor (2nd or 4th order), one can project it onto the closest tensor with a prescribed symmetry class (ISO, TI, ORTHO). The projection minimises the Frobenius distance $\lVert B - A \rVert$ over all tensors $B$ of the target symmetry.

The function proj_tens provides this projection. It returns a 3-tuple (B, d, drel):

B: the projected tensor
d: absolute Frobenius distance $\lVert B - A \rVert$
drel: relative distance $d / \lVert A \rVert$

Fixed-axis TI projection (order 4)

Project a 4th-order tensor onto the TI subspace with a given axis $\mathbf{n}$. The result is a TensWalpole{T,5}:

julia> using TensND, LinearAlgebra
julia> n = [0., 0., 1.] ;
julia> C = tensTI(10., 3., 2.5, 12., 2., n) ;
julia> B, d, drel = proj_tens(:TI, getarray(C), n) ;
julia> drel < 1e-12true
julia> B isa TensWalpoletrue

Fixed-axis TI projection (order 2)

For a 2nd-order tensor, the projection returns a TensTI{2}:

julia> A = Float64[5 1 0; 1 5 0; 0 0 8] ;
julia> B2, d2, drel2 = proj_tens(:TI, A, n) ;
julia> B23×3 TensTI{2, Float64, 2}:
 5.0  0.0  0.0
 0.0  5.0  0.0
 0.0  0.0  8.0
julia> B2.data(5.0, 8.0)

Fixed-frame ORTHO projection (order 4)

Project onto the orthotropic subspace with a given material frame:

julia> frame = CanonicalBasis{3,Float64}() ;
julia> t = TensOrtho(10., 8., 12., 3., 2.5, 1.5, 2., 3., 3.5, frame) ;
julia> Bo, do_, drelo = proj_tens(:ORTHO, getarray(t), frame) ;
julia> drelo < 1e-12true
julia> Bo isa TensOrthotrue

Fixed-frame ORTHO projection (order 2)

julia> M = Float64[5 1 2; 1 8 3; 2 3 12] ;
julia> Bm, dm, drelm = proj_tens(:ORTHO, M, frame) ;
julia> Bm3×3 Matrix{Float64}:
 5.0  0.0   0.0
 0.0  8.0   0.0
 0.0  0.0  12.0
julia> dm > 0true

Best symmetry detection

The function best_sym_tens tries symmetries from the most restrictive to the least (ISO → TI → ORTHO) and returns the first whose relative projection error falls below a threshold $\varepsilon$ (default 1e-6). Pass an axis or frame for fixed-basis detection:

julia> C_ti = tensTI(10., 3., 2.5, 12., 2., n) ;
julia> _, _, _, sym = best_sym_tens(C_ti, n; proj = (:TI, :ORTHO)) ;
julia> sym:TI

Rotation-optimized projection (requires NLopt)

When no axis or frame is provided, proj_tens optimises the rotation angles to find the best approximation. This requires the NLopt package (declared as a weak dependency):

using NLopt   # triggers the TensNDNLoptExt extension

# Optimise the TI axis automatically
B_opt, d_opt, drel_opt = proj_tens(:TI, getarray(C))

# Optimise the ORTHO frame automatically
B_ort, d_ort, drel_ort = proj_tens(:ORTHO, getarray(C))

The optimiser uses a two-pass strategy (global + local refinement) inspired by the ECHOES library, with gradients computed via ForwardDiff.