Theory
Bases and variance, the tensor products, Kelvin-Mandel storage, rotations, and the symmetry classes with their Walpole algebra.
Bases and variance, structured tensor types, symmetry projection, and differential operators in curvilinear coordinates — symbolic or numeric, one implementation.

TensND handles tensors of any order in any coordinate system. A basis carries its metric, so covariant and contravariant components are related rather than assumed; the differential operators follow from the chart, symbolically or by automatic differentiation.
Three structured types — TensISO, TensTI and TensCubic, TensOrtho — store the 2, 5, 3 and 9 scalars a symmetry class really has, and compute products and inverses in closed form. The same code runs on Float64, ForwardDiff.Dual, SymPy.Sym and Symbolics.Num.
The design is inspired by the Maple library Tens3d of Jean Garrigues.
using TensND, SymPy
Spherical = coorsys_spherical()
θ, ϕ, r = getcoords(Spherical)
𝐞ᶿ, 𝐞ᵠ, 𝐞ʳ = unitvec(Spherical)
@set_coorsys Spherical
σʳʳ = SymFunction("σʳʳ", real = true)(r)
σᶿᶿ = SymFunction("σᶿᶿ", real = true)(r)
𝛔 = σʳʳ * 𝐞ʳ ⊗ 𝐞ʳ + σᶿᶿ * (𝐞ᶿ ⊗ 𝐞ᶿ + 𝐞ᵠ ⊗ 𝐞ᵠ)
pprint(DIV(𝛔))(Derivative(σʳʳ(r), r) + 2*σʳʳ(r)/r - 2*σᶿᶿ(r)/r)𝐞ʳThe equilibrium equation of a spherically symmetric stress state, derived rather than transcribed.