Symbolic helpers and array algebra

Element-wise symbolic operations — no-ops on numeric types — and the raw-array products the whole library is built on. Usage: Symbolic and numeric; theory: Tensor algebra.

TensND.tsimplify Function

tsimplify(x, args...; kwargs...)

Simplify a scalar, an array or a tensor elementwise.

Dispatches to SymPy.simplify for Sym, to Symbolics.simplify for Num, and is the identity on every other type — so a routine written to be simplification-aware runs unchanged, and at full speed, on numeric input.

See also tfactor, tsubs, tdiff, ttrigsimp, texpand_trig.

Examples

julia> tsimplify(3.0)          # no-op on numbers
3.0

julia> tsimplify([1.0, 2.0])
2-element Vector{Float64}:
 1.0
 2.0

TensND.tfactor Function

tfactor(x, args...; kwargs...)

Factorize an expression elementwise (SymPy.factor); identity on non-symbolic types. See tsimplify.

TensND.tsubs Function

tsubs(x, substitutions...)

Substitute into an expression elementwise (SymPy.subs / Symbolics.substitute); identity on non-symbolic types.

Substitutions are given as old => new pairs. See tsimplify.

TensND.tdiff Function

tdiff(y, x...; kwargs...)

Differentiate an expression elementwise with respect to one or more symbols (SymPy.diff / Symbolics.Differential); identity on non-symbolic types.

Repeated symbols give higher derivatives. See tsimplify.

TensND.ttrigsimp Function

ttrigsimp(x, args...; kwargs...)

Trigonometric simplification elementwise (sympy.trigsimp); identity on non-symbolic types. Often succeeds where tsimplify does not, on expressions built from a rotated basis.

TensND.texpand_trig Function

texpand_trig(x, args...; kwargs...)

Trigonometric expansion elementwise (sympy.expand_trig); identity on non-symbolic types. The counterpart of ttrigsimp.

TensND.tlimit Function

tlimit(x, s, v)
tlimit(x, s, v, dir)

Limit of an expression elementwise as the symbol s tends to v (sympy.limit), optionally one-sided with dir = "+" or "-"; identity on non-symbolic types, a numeric value having no free symbol to send anywhere.

This is the pass a degenerate limit of a closed form needs — a vanishing regularization, an incompressible phase (k => oo), a flat inclusion — and it keeps the tensor's structure: applied to a TensISO, a TensTI or a TensOrtho it takes the limit of the few canonical coefficients and rebuilds the same type, never going through the dim^order components.

Unlike its siblings there is no Symbolics.Num method: Symbolics has no limit, and returning the input unchanged would be a silent wrong answer, so a Num raises instead.

See also tsimplify, tsubs.

Examples

julia> tlimit(3.0, :whatever, 0)      # no-op on numbers
3.0

TensND.ApproxType Type

ApproxType

Element types whose symmetry predicates use a tolerance rather than exact equality: floating-point numbers, and ForwardDiff.Dual numbers built on them.

ForwardDiff.Dual <: Real but not <: AbstractFloat, so before this union existed a Dual-valued tensor fell through to the exact fallback, where a round-off of a few ulp was enough to report it as non-minor-symmetric. The consequence was not local: _KM_of_array then built a 9×9 matrix instead of a 6×6 one, and every proj_tens call on a Dual input died with a DimensionMismatch — automatic differentiation through a projection was impossible. Symbolic types keep their own exact methods and are unaffected.

Companion of is_hard_numeric, which answers a different question: ApproxType says use a tolerance rather than exact equality, is_hard_numeric says a comparison is allowed at all.

TensND.is_hard_numeric Function

is_hard_numeric(::Type) -> Bool
is_hard_numeric(x)      -> Bool

Whether comparisons on that scalar type yield an honest Bool, so a value of it may drive if, && or a tolerance test.

T <: Real is not this predicate. Symbolics.Num <: Real, yet <, ≈ and iszero on a Num return symbolic expressions, and using one in a boolean context throws TypeError: non-boolean (Num) used in boolean context. ForwardDiff.Dual is also <: Real and does compare to a Bool, so the two cannot be separated by <: Real either way. SymPy.Sym is not <: Real at all, which is why the bug it guards against only ever showed up on the Symbolics side.

The list is explicit and closed, and the default is the safe answer: an unknown scalar type is assumed not comparable, which pushes callers onto their structural (isequal) branch. Add a type here only after checking that its < really returns a Bool.

This is the companion of ApproxType: ApproxType says use a tolerance rather than exact equality, is_hard_numeric says a comparison is allowed at all. Int and Rational are hard numeric but not ApproxType.

julia> TensND.is_hard_numeric(Float64), TensND.is_hard_numeric(Sym)
(true, false)

TensND.contract Function

contract(t::AbstractArray, i::Integer, j::Integer)

Contract (trace over) indices i and j of a single tensor, lowering its order by two. For an order-2 array this is the ordinary trace.

Examples

julia> contract([1.0 2.0; 3.0 4.0], 1, 2)
5.0

contract(t::AbstractTens{order,dim}, i::Integer, j::Integer)

Calculate the tensor obtained after contraction with respect to the indices i and j

TensND.qcontract Function

qcontract(t1::AbstractArray, t2::AbstractArray)
t1 ⊙ t2

Quadruple contraction: the last four indices of t1 against the first four of t2, lowering the total order by eight.

Between two order-4 tensors it is the Frobenius scalar product T ⊙ T' = Tᵢⱼₖₗ T'ᵢⱼₖₗ — the inner product every projection in proj_tens minimizes against.

Examples

julia> 𝕀, 𝕁, 𝕂 = ISO(Val(3), Val(Float64));

julia> (𝕁 ⊙ 𝕁, 𝕂 ⊙ 𝕂, 𝕁 ⊙ 𝕂)
(1.0, 5.0, 0.0)

𝕁 ⊙ 𝕁 = 1 and 𝕂 ⊙ 𝕂 = 5 are the dimensions of the spherical and deviatoric subspaces in 3-D.

See also dcontract, contract.

qcontract(t1::AbstractTens{order1,dim}, t2::AbstractTens{order2,dim})

Define a quadruple contracted product between two tensors

𝔸 ⊙ 𝔹 = AᵢⱼₖₗBⁱʲᵏˡ

Examples

julia> 𝕀 = t𝕀(Sym) ; 𝕁 = t𝕁(Sym) ; 𝕂 = t𝕂(Sym) ;

julia> 𝕀 ⊙ 𝕀
6

julia> 𝕁 ⊙ 𝕀
1

julia> 𝕂 ⊙ 𝕀
5

julia> 𝕂 ⊙ 𝕁
0

Tensors.otimesu Function

otimesu(t1::AbstractTens{order1,dim}, t2::AbstractTens{order2,dim})

Define a special tensor product between two tensors of at least second order

(𝐚 ⊠ 𝐛) ⊡ 𝐩 = 𝐚⋅𝐩⋅𝐛 = aⁱᵏbʲˡpₖₗ eᵢ⊗eⱼ

TensND.otimesul Function

otimesul(t1::AbstractArray, t2::AbstractArray)
sboxtimes(t1, t2)
t1 ⊠ˢ t2

Symmetrized modified tensor product, (a ⊠ˢ b)ᵢⱼₖₗ = (aᵢₖ bⱼₗ + aᵢₗ bⱼₖ) / 2, i.e. the half-sum of otimesu and otimesl.

𝕀 = 𝟏 ⊠ˢ 𝟏 is the identity of the symmetric order-2 tensors, which is what tens_Id4 returns.

It does not invert termwise

(a ⊠ˢ b)⁻¹ ≠ a⁻¹ ⊠ˢ b⁻¹ unless a and b are proportional. Commuting is not sufficient, and taking b = 𝟏 does not help. This is why inversion is implemented per symmetry class rather than by a single generic formula.

otimesul(t1::AbstractTens{order1,dim}, t2::AbstractTens{order2,dim})

Define a special tensor product between two tensors of at least second order

(𝐚 ⊠ˢ 𝐛) ⊡ 𝐩 = (𝐚 ⊠ 𝐛) ⊡ (𝐩 + ᵗ𝐩)/2 = 1/2(aⁱᵏbʲˡ+aⁱˡbʲᵏ) pₖₗ eᵢ⊗eⱼ

TensND.sotimes Function

sotimes(t1::AbstractArray, t2::AbstractArray)
t1 ⊗ˢ t2

Tensor product symmetrized over the last index of t1 and the first of t2. For two vectors, u ⊗ˢ v = (u ⊗ v + v ⊗ u) / 2.

It is the product used by SYMGRAD, which is why that operator returns a linearized strain tensor directly.

See also otimes, otimesul.

sotimes(t1::AbstractTens{order1,dim}, t2::AbstractTens{order2,dim})

Define a symmetric tensor product between two tensors

(aⁱeᵢ) ⊗ˢ (bʲeⱼ) = 1/2(aⁱbʲ + aʲbⁱ) eᵢ⊗eⱼ

TensND.sboxtimes Function

otimesul(t1::AbstractArray, t2::AbstractArray)
sboxtimes(t1, t2)
t1 ⊠ˢ t2

Symmetrized modified tensor product, (a ⊠ˢ b)ᵢⱼₖₗ = (aᵢₖ bⱼₗ + aᵢₗ bⱼₖ) / 2, i.e. the half-sum of otimesu and otimesl.

𝕀 = 𝟏 ⊠ˢ 𝟏 is the identity of the symmetric order-2 tensors, which is what tens_Id4 returns.

It does not invert termwise

(a ⊠ˢ b)⁻¹ ≠ a⁻¹ ⊠ˢ b⁻¹ unless a and b are proportional. Commuting is not sufficient, and taking b = 𝟏 does not help. This is why inversion is implemented per symmetry class rather than by a single generic formula.

otimesul(t1::AbstractTens{order1,dim}, t2::AbstractTens{order2,dim})

Define a special tensor product between two tensors of at least second order

(𝐚 ⊠ˢ 𝐛) ⊡ 𝐩 = (𝐚 ⊠ 𝐛) ⊡ (𝐩 + ᵗ𝐩)/2 = 1/2(aⁱᵏbʲˡ+aⁱˡbʲᵏ) pₖₗ eᵢ⊗eⱼ

TensND.isidentity Function

isidentity(a::AbstractMatrix) → Bool

Test whether a is the identity matrix, within the tolerance of ≈ for floating-point element types and exactly for symbolic ones.

Examples

julia> isidentity([1.0 0.0; 0.0 1.0])
true

julia> isidentity([1.0 1e-20; 0.0 1.0])
true

julia> isidentity([1.0 0.5; 0.0 1.0])
false