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 SymPyPythonCall.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 (SymPyPythonCall.factor); identity on non-symbolic types. See tsimplify.

TensND.tsubs Function

tsubs(x, substitutions...)

Substitute into an expression elementwise (SymPyPythonCall.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 (SymPyPythonCall.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 project 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. SymPyPythonCall.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.carries_partials Function

carries_partials(x) -> Bool

Whether x carries a nonzero derivative, at any nesting level of ForwardDiff.Dual. A plain number, or a symbolic one, carries none.

The structural predicates isidentity and isdiagonal decide whether a basis collapses to a simpler one, and they compare values — ≈ and norm see only the value of a Dual. A rotation by a dual angle at zero is the identity in value, yet its derivative is the generator of the rotation: collapsing it to the canonical basis returned zero derivatives. The predicates therefore also require the derivatives to vanish.

TensND.primal Function

primal(x)

x stripped of every level of ForwardDiff.Dual: the value a choice that must not depend on derivatives is made on (a quadrature order, a scaling, a bracket, a pivot), so that the derivatives stay where the mathematics puts them. Any other number, symbolic ones included, is returned as it is.

See also carries_partials.

TensND.primal_type Function

primal_type(T) -> Type

The floating-point type under every level of ForwardDiff.Dual in T: the type of float(primal(x)) for x::T, to build a tolerance, a typemax or quadrature weights that carry no derivative. An integer or rational type gives its float type; any other (a complex or symbolic type) is returned as it is.

TensND.dual_depth Function

dual_depth(T) -> Int

The number of levels of ForwardDiff.Dual nested in T: 0 for a plain number, 2 for a second derivative taken through nested ForwardDiff calls.

TensND.dual_channels Function

dual_channels(x) -> Tuple

The value of x followed by every partial, at every nesting level of ForwardDiff.Dual, flattened into one tuple: (x,) for a plain number. A convergence test that must hold for the derivatives as well as for the value measures each channel separately.

TensND.float_like Function

float_like(T) -> Type

float(T) for an integer type, T itself otherwise: a float, a dual or a symbolic type stays as it is (float of a SymPy type answers Float64, which a symbolic quantity cannot be converted to).

TensND.isparallel Function

isparallel(a, b; rtol = 1.0e-10) -> Bool

Whether the directions a and b (3-vectors, tuples or 1st-order tensors, read in canonical components) are parallel, in the same or the opposite sense:    .

For comparable numbers (is_hard_numeric), within the relative tolerance rtol; a ForwardDiff.Dual is compared on its value. For symbolic ones the test is structural: the difference is simplified and must vanish exactly, so a symbolic answer is never a false positive, and is false when the backend cannot decide. No dot and no abs: both conjugate a SymPyPythonCall.Sym, and Symbolics does not fold abs of a literal.

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 project minimizes against.

Examples

julia> 𝕀, 𝕁, 𝕂 = iso_projectors(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.dotdot Function

dotdot(v1::AbstractTens{order1,dim}, S::AbstractTens{orderS,dim}, v2::AbstractTens{order2,dim})

Define a bilinear operator 𝐯₁⋅𝕊⋅𝐯₂

Examples

julia> n = Tens(Sym[0, 0, 1]) ;

julia> k, μ = symbols("k μ", real =true) ;

julia> ℂ = 3k * t𝕁() + 2μ * t𝕂() ;

julia> dotdot(n,ℂ,n) # Acoustic tensor
3×3 Tens{2, 3, Sym, Sym, Tensor{2, 3, Sym, 9}, CanonicalBasis{3, Sym}}:
 μ  0          0
 0  μ          0
 0  0  k + 4*μ/3

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. A matrix of ForwardDiff.Dual numbers is the identity only if its derivatives vanish too.

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