Coordinate systems and differential operators

Symbolic coordinate systems (`CoorSystemSym`)

Symbolic coordinate systems support exact derivation of differential operators using SymPy. A CoorSystemSym stores the position vector OM, coordinate symbols, the natural basis $\mathbf{a}_i = \partial_i \mathbf{OM}$, the normalized (unit) basis $\mathbf{e}_i = \mathbf{a}_i / \|\mathbf{a}_i\|$, Lame coefficients $\chi_i = \|\mathbf{a}_i\|$, and Christoffel symbols $\Gamma_{ij}^k = \partial_i \mathbf{a}_j \cdot \mathbf{a}^k$.

Predefined symbolic systems

Constructor	Coordinates	Description
`coorsys_cartesian()`	$(x, y)$ or $(x, y, z)$	Cartesian
`coorsys_polar()`	$(r, \theta)$	Polar
`coorsys_cylindrical()`	$(r, \theta, z)$	Cylindrical
`coorsys_spherical()`	$(\theta, \varphi, r)$	Spherical
`coorsys_spheroidal(c)`	$(\varphi, p, q)$	Prolate spheroidal

Example: polar Laplacian

julia> using TensND, SymPy
julia> Polar = coorsys_polar() ; r, θ = getcoords(Polar) ; 𝐞ʳ, 𝐞ᶿ = unitvec(Polar) ;
julia> @set_coorsys Polar
julia> LAPLACE(SymFunction("f", real = true)(r, θ))                              2
                             ∂
               ∂             ───(f(r, θ))
 2             ──(f(r, θ))     2
∂              ∂r            ∂θ
───(f(r, θ)) + ─────────── + ────────────
  2                 r              2
∂r                                r
julia> n = symbols("n", integer = true)n
julia> simplify(HESS(r^n))2×2 TensND.TensRotated{2, 2, Sym{PyCall.PyObject}, Tensors.SymmetricTensor{2, 2, Sym{PyCall.PyObject}, 3}}:
 n^2*r^n/r^2 - n*r^n/r^2          0
                       0  n*r^n/r^2

Example: spherical divergence

julia> using TensND, SymPy
julia> Spherical = coorsys_spherical() ; θ, ϕ, r = getcoords(Spherical) ; 𝐞ᶿ, 𝐞ᵠ, 𝐞ʳ = unitvec(Spherical) ;
julia> @set_coorsys Spherical
julia> Christoffel(Spherical)3×3×3 Array{Sym{PyCall.PyObject}, 3}:
[:, :, 1] =
   0            0  1/r
   0  -sin(2*θ)/2    0
 1/r            0    0

[:, :, 2] =
        0  1/tan(θ)    0
 1/tan(θ)         0  1/r
        0       1/r    0

[:, :, 3] =
 -r            0  0
  0  -r*sin(θ)^2  0
  0            0  0
julia> ℬˢ = normalized_basis(Spherical)3×3 RotatedBasis{3, Sym{PyCall.PyObject}}:
 cos(θ)⋅cos(ϕ)  -sin(ϕ)  sin(θ)⋅cos(ϕ)
 sin(ϕ)⋅cos(θ)   cos(ϕ)  sin(θ)⋅sin(ϕ)
       -sin(θ)        0         cos(θ)
julia> σʳʳ = SymFunction("σʳʳ", real = true)(r) ;
julia> σᶿᶿ = SymFunction("σᶿᶿ", real = true)(r) ;
julia> σᵠᵠ = SymFunction("σᵠᵠ", real = true)(r) ;
julia> 𝛔 = σʳʳ * 𝐞ʳ ⊗ 𝐞ʳ + σᶿᶿ * 𝐞ᶿ ⊗ 𝐞ᶿ + σᵠᵠ * 𝐞ᵠ ⊗ 𝐞ᵠ3×3 TensND.TensRotated{2, 3, Sym{PyCall.PyObject}, Tensors.SymmetricTensor{2, 3, Sym{PyCall.PyObject}, 6}}:
 σᶿᶿ(r)       0       0
      0  σᵠᵠ(r)       0
      0       0  σʳʳ(r)
julia> div𝛔 = simplify(DIV(𝛔))3-element TensND.TensRotated{1, 3, Sym{PyCall.PyObject}, Tensors.Vec{3, Sym{PyCall.PyObject}}}:
 r*(-σᵠᵠ(r)*sin(2*θ)/(2*r^2*sin(θ)^2) + σᶿᶿ(r)/(r^2*tan(θ)))
                                                           0
    Derivative(σʳʳ(r), r) + 2*σʳʳ(r)/r - σᵠᵠ(r)/r - σᶿᶿ(r)/r

Differential operators

The following operators are available after calling @set_coorsys CS:

Operator	Signature	Description
`GRAD(f)`	scalar or tensor field	Gradient
`SYMGRAD(v)`	vector field	Symmetric gradient $(\nabla \mathbf{v} + \nabla \mathbf{v}^T)/2$
`DIV(t)`	tensor field (order >= 1)	Divergence
`LAPLACE(f)`	scalar or tensor field	Laplacian $\nabla^2 f$
`HESS(f)`	scalar field	Hessian $\nabla \nabla f$

Accessors

getcoords(CS) / getcoords(CS, i): coordinate symbols
getOM(CS): position vector
normalized_basis(CS): orthonormal basis $(\mathbf{e}_i)$
natural_basis(CS): natural basis (Basis type including $\mathbf{a}_i$ and $\mathbf{a}^i$)
unitvec(CS) / unitvec(CS, i): unit basis vectors as tensors
natvec(CS, i, :cov) / natvec(CS, i, :cont): natural basis vectors
Lame(CS) / Lame(CS, i): Lame coefficients $\chi_i$
Christoffel(CS): Christoffel symbols as a 3D array

Numerical coordinate systems (`CoorSystemNum`)

Numerical coordinate systems evaluate differential operators pointwise using automatic differentiation via ForwardDiff. They do not require a symbolic setup and work with any numeric type (including ForwardDiff.Dual for nested differentiation).

Predefined numerical systems

Constructor	Coordinates	Description
`coorsys_cartesian_num(dim)`	$(x_1, \ldots, x_d)$	Cartesian
`coorsys_polar_num()`	$(r, \theta)$	Polar
`coorsys_cylindrical_num()`	$(r, \theta, z)$	Cylindrical
`coorsys_spherical_num()`	$(\theta, \varphi, r)$	Spherical

Custom numerical systems

A CoorSystemNum can be built from any position vector function OM(x):

CS = CoorSystemNum(x -> [x[1]*cos(x[2]), x[1]*sin(x[2])], 2)  # polar coordinates

Pointwise evaluation

All operators are functions that return functions (lazy evaluation). They must be called with the coordinate point as argument:

CS = coorsys_polar_num()
f(x) = x[1]^2              # r^2
lap_f = LAPLACE(CS, f)      # returns a function
lap_f([2.0, π/4])           # evaluate at (r=2, theta=pi/4) => 4.0

Operators (same as symbolic, pointwise)

Operator	Signature	Returns
`GRAD(CS, f)`	scalar/tensor field function	function `x -> gradient`
`SYMGRAD(CS, v)`	vector field function	function `x -> sym. gradient`
`DIV(CS, t)`	tensor field function	function `x -> divergence`
`LAPLACE(CS, f)`	scalar/tensor field function	function `x -> Laplacian`
`HESS(CS, f)`	scalar field function	function `x -> Hessian`

Pointwise accessors

normalized_basis(CS, x0): orthonormal basis at point x0
unitvec(CS, x0, i): unit vector $\mathbf{e}_i$ at point x0
natvec(CS, x0, i, :cov): natural covariant vector at point x0
Lame(CS, x0) / Lame(CS, x0, i): Lame coefficients at point x0

For a detailed tutorial with validation examples (polar, spherical, cylindrical, Lame problem on a hollow sphere), see Numerical differential operators.