Coordinate systems (symbolic)

Symbolic differential calculus on a curvilinear chart. Theory: Curvilinear differential calculus; usage: Coordinate systems.

TensND.CoordinateSystem Type

CoordinateSystem{dim,T,G<:AbstractGeometry} <: AbstractCoorSystem{dim,T}

Curvilinear coordinate system of the dim-dimensional space, of one of two geometries:

Geometry	Fields	Operators
ExactGeometry	expressions of the coordinate symbols (Sym or Num)	take and return expressions, closed-form results valid everywhere
PointwiseGeometry	functions of the point	take a function of the point, return one, by automatic differentiation

The differential operators grad, symgrad, divergence, laplacian and hess have the same names and meaning on both.

Predefined systems: coorsys_cartesian, coorsys_polar, coorsys_cylindrical, coorsys_spherical, coorsys_spheroidal — exact by default, pointwise with PointwiseGeometry as first argument (e.g. coorsys_polar(PointwiseGeometry)). @set_coorsys makes an exact one the default so the operators take a single argument.

Construction

CoordinateSystem(OM, coords, tmp_coords = (), params = ();
                 rules = Dict(), tmp_var = Dict(), to_coords = Dict())
CoordinateSystem(OM, coords, bnorm, χᵢ, tmp_coords = (), params = (); ...)
CoordinateSystem(OM_func::Function, dim::Integer, T::Type = Float64)

The first two build an exact system from the position vector OM, an expression of the coordinate symbols coords. The first derives the natural basis from 𝐚ᵢ = ∂ᵢOM; the second takes the normalized basis and the Lamé coefficients directly, which is much faster when they are known in closed form. The third builds a pointwise system from OM_func, a function of the point returning its Cartesian components, by automatic differentiation.

The optional arguments of the exact forms drive symbolic simplification, without which a non-trivial chart produces correct but unusable nested radicals:

Argument	Role
tmp_coords	auxiliary symbols standing for compound expressions
params	constants appearing in OM
tmp_var	substitutions replacing expressions by those symbols
to_coords	how to eliminate them again before differentiating
rules	rewrite rules applied after each simplification

See Adding a coordinate system for the recipe and coorsys_spheroidal for a worked example.

One type since 0.8

The former CoorSystemSym and CoorSystemNum are the two geometries of this one type; their names were removed in 0.11.

TensND.AbstractGeometry Type

AbstractGeometry

What a CoordinateSystem knows of its chart, and so how its differential operators work: ExactGeometry holds closed-form expressions of the coordinate symbols, PointwiseGeometry functions of the point.

TensND.ExactGeometry Type

ExactGeometry{dim,T,VEC,BNORM,BNAT} <: AbstractGeometry

The geometry of a chart in closed form, as expressions of its coordinate symbols (SymPy's Sym or Symbolics' Num): the position vector, the coordinates, the two bases, the natural and dual frames, the Lamé coefficients, the connection array, and the symbolic-simplification settings. Differential operators on it take and return expressions.

Held by both a CoordinateSystem and a Submanifold — which differ only in what they add to it (a submanifold adds the two fundamental forms) and in how many directions they differentiate along (nderiv). Declaring the block once keeps the two in step: a fix to one is a fix to the other.

Reached through core (internal); user code should go through the accessors (get_coords, lame_coefficients, christoffel, …) instead.

TensND.@set_coorsys Macro

@set_coorsys CS [vec] [coords]

Make CS the default coordinate system for the differential operators, so they can be called with one argument.

Equivalent to set_coorsys!(CS; vec, coords); kept as a macro for familiarity. Unlike earlier versions it does not define any method — it only stores the system — and it leaves the meaning of ∂(t, x) untouched.

Examples

julia> Spherical = coorsys_spherical() ; θ, ϕ, r = get_coords(Spherical) ;

julia> @set_coorsys Spherical

julia> laplacian(1 / r)      # one argument
0

See also set_coorsys!, default_coorsys, unset_coorsys!.

TensND.set_coorsys! Function

set_coorsys!(CS; vec = '𝐞', coords = nothing) → CS

Make CS the default coordinate system, so that grad, symgrad, divergence, laplacian, hess and pprint may be called with a single argument.

vec and coords control how pprint names the basis vectors. The macro @set_coorsys is a thin wrapper over this function.

∂ is deliberately not affected

∂(t, x) always means the plain derivative of t with respect to the symbol x, whether or not a default chart is set. For the covariant derivative, pass the system explicitly: ∂(t, x, CS).

See also default_coorsys, unset_coorsys!.

TensND.default_coorsys Function

default_coorsys() → AbstractCoorSystem

The coordinate system installed by set_coorsys! or @set_coorsys. Throws if none has been set.

TensND.unset_coorsys! Function

unset_coorsys!()

Forget the default coordinate system installed by set_coorsys!, with the letter and the coordinate names it gave pprint.

TensND.grad Function

grad(t::Union{T,AbstractTens{order,dim,T}}, CS::CoordinateSystem{dim,T,<:ExactGeometry}) where {order,dim,T<:SymType}

Gradient of the scalar or tensor field t with respect to the coordinate system CS, raising the order by one:

grad(t) = Σᵢ ∂ᵢt ⊗ 𝐚ⁱ

The derivative index comes last

The dual natural vector 𝐚ⁱ is appended on the right, so for a vector field 𝐯 the components of grad(𝐯) are (∇𝐯)ᵢⱼ = ∂ⱼvᵢ. This is the convention that makes divergence the contraction of the last index and laplacian = divergence ∘ grad come out right; a library using the opposite convention differs by a transpose.

If CS has been made the default with @set_coorsys, the second argument may be omitted.

Examples

julia> Spherical = coorsys_spherical() ; θ, ϕ, r = get_coords(Spherical) ;

julia> @set_coorsys Spherical

julia> grad(r)          # = 𝐞ʳ

julia> 𝐞ᶿ, 𝐞ᵠ, 𝐞ʳ = unitvec(Spherical) ; grad(𝐞ʳ)   # = (𝐞ᶿ⊗𝐞ᶿ + 𝐞ᵠ⊗𝐞ᵠ)/r

See also symgrad, divergence, laplacian, hess.

grad(f::Function, CS::_PointwiseSystem) -> Function

Return a function x₀ -> gradient as an AbstractTens in the normalized basis at x₀.

• If f returns a scalar, the gradient is a rank-1 AbstractTens.

• If f returns a rank-n AbstractTens or array, the gradient is a rank-n+1 AbstractTens.

TensND.symgrad Function

symgrad(t::Union{T,AbstractTens{order,dim,T}}, CS::CoordinateSystem{dim,T,<:ExactGeometry}) where {order,dim,T<:SymType}

Symmetrized gradient of t with respect to the coordinate system CS:

symgrad(t) = Σᵢ ∂ᵢt ⊗ˢ 𝐚ⁱ

Applied to a displacement field this is the linearized strain tensor, 𝛆 = (∇𝛏 + ᵗ∇𝛏)/2, which is why it is a primitive rather than a composition of grad and a transpose.

If CS has been made the default with @set_coorsys, the second argument may be omitted.

Examples

julia> Cylindrical = coorsys_cylindrical() ; r, θ, z = get_coords(Cylindrical) ;

julia> 𝐞ʳ, 𝐞ᶿ, 𝐞ᶻ = unitvec(Cylindrical) ; @set_coorsys Cylindrical

julia> 𝛏 = SymFunction("ξʳ", real = true)(r, z) * 𝐞ʳ + SymFunction("ξᶻ", real = true)(r, z) * 𝐞ᶻ ;

julia> symgrad(𝛏)       # axisymmetric strain tensor, with εᶿᶿ = ξʳ/r

See also grad, divergence.

symgrad(f::Function, CS::_PointwiseSystem) -> Function

Return a function x₀ -> symmetric gradient as an AbstractTens{2} in the normalized basis at x₀. Applies to vector-valued functions f.

Tensors.divergence Function

divergence(t::AbstractTens{order,dim,T}, CS::CoordinateSystem{dim,T,<:ExactGeometry}) where {order,dim,T<:SymType}

Divergence of the tensor field t (of order ≥ 1) with respect to the coordinate system CS, lowering the order by one:

divergence(t) = Σᵢ ∂ᵢt ⋅ 𝐚ⁱ

The contracted index is the last one

Consistently with grad, the contraction acts on the last index of t: for an order-2 field, divergence(𝛔)ᵢ = ∂ⱼσᵢⱼ. For a symmetric field such as a stress tensor the distinction is immaterial, but it is not for a general order-2 field.

If CS has been made the default with @set_coorsys, the second argument may be omitted.

Examples

julia> Spherical = coorsys_spherical() ; θ, ϕ, r = get_coords(Spherical) ;

julia> 𝐞ᶿ, 𝐞ᵠ, 𝐞ʳ = unitvec(Spherical) ; @set_coorsys Spherical

julia> divergence(𝐞ʳ)          # = 2/r

julia> 𝛔 = SymFunction("σʳʳ", real = true)(r) * 𝐞ʳ ⊗ 𝐞ʳ ;

julia> simplify(divergence(𝛔)) # radial equilibrium operator

See also grad, laplacian.

divergence(f::Function, CS::_PointwiseSystem) -> Function

Return a function x₀ -> divergence as a scalar (if f is a vector field) or AbstractTens (if f is a tensor field of order ≥ 2).

f may return a plain Array or an AbstractTens (components are extracted via Array).

TensND.laplacian Function

laplacian(t::Union{T,AbstractTens{order,dim,T}}, CS::CoordinateSystem{dim,T,<:ExactGeometry}) where {order,dim,T<:SymType}

Laplacian of the scalar or tensor field t, defined as the composition

laplacian(t) = divergence(grad(t))

and therefore preserving the order of t.

If CS has been made the default with @set_coorsys, the second argument may be omitted.

Examples

julia> Polar = coorsys_polar() ; r, θ = get_coords(Polar) ; @set_coorsys Polar

julia> laplacian(SymFunction("f", real = true)(r, θ))   # the polar Laplacian

julia> n = symbols("n", integer = true) ; simplify(laplacian(r^n * cos(n*θ)))   # = 0

See also grad, divergence, hess.

laplacian(f::Function, CS::_PointwiseSystem) -> Function

Return a function x₀ -> Laplacian (scalar). Uses the raw (plain-array) gradient and divergence internally, so it is fully compatible with ForwardDiff differentiation.

TensND.hess Function

hess(t::Union{T,AbstractTens{order,dim,T}}, CS::CoordinateSystem{dim,T,<:ExactGeometry}) where {order,dim,T<:SymType}

Hessian of t, defined as the second gradient

hess(t) = grad(grad(t))

and therefore raising the order of t by two. For a scalar field the trace of hess(t) is laplacian(t).

If CS has been made the default with @set_coorsys, the second argument may be omitted.

Examples

julia> Spherical = coorsys_spherical() ; θ, ϕ, r = get_coords(Spherical) ;

julia> @set_coorsys Spherical

julia> simplify(hess(1/r))    # the kernel of the 3-D Laplace equation

See also grad, laplacian.

hess(f::Function, CS::_PointwiseSystem) -> Function

Return a function x₀ -> Hessian as an AbstractTens{2} in the normalized basis. Uses the raw (plain-array) double gradient internally.

TensND.∂ Function

∂(t::AbstractTens{order,dim,T,A},xᵢ::T)

Return the derivative of the tensor t with respect to the variable x_i

Examples

julia> (θ, ϕ, r), (𝐞ᶿ, 𝐞ᵠ, 𝐞ʳ), ℬˢ = init_spherical() ;

julia> ∂(𝐞ʳ, ϕ) == sin(θ) * 𝐞ᵠ
true

julia> ∂(𝐞ʳ ⊗ 𝐞ʳ,θ)
Tens.TensRotated{2, 3, Sym, SymmetricTensor{2, 3, Sym, 6}}
# data: 3×3 SymmetricTensor{2, 3, Sym, 6}:
 0  0  1
 0  0  0
 1  0  0
# basis: 3×3 Tensor{2, 3, Sym, 9}:
 cos(θ)⋅cos(ϕ)  -sin(ϕ)  sin(θ)⋅cos(ϕ)
 sin(ϕ)⋅cos(θ)   cos(ϕ)  sin(θ)⋅sin(ϕ)
       -sin(θ)        0         cos(θ)
# var: (:cont, :cont)

∂(f::Function, i::Integer, CS::_PointwiseSystem, x₀)

Covariant derivative of f with respect to the i-th coordinate, evaluated at x₀.

f must accept an AbstractVector of length dim and return either a scalar or a plain Array of physical (normalized-frame) components of the tensor field. Returns a plain scalar or plain Array (suitable for further ForwardDiff differentiation).

See also grad, divergence, laplacian, hess, symgrad for high-level operators that wrap results in AbstractTens.

TensND.get_coords Function

get_coords(CS::AbstractCoorSystem) → NTuple
get_coords(CS::AbstractCoorSystem, i::Integer)

Coordinate symbols of CS, or the i-th one.

Note the ordering of the spherical system, (θ, ϕ, r) and not (r, θ, ϕ), so that θ = ϕ = 0 reproduces the canonical basis in the canonical order.

TensND.position_vector Function

position_vector(CS::AbstractCoorSystem) → AbstractTens{1}

Position vector of CS as a function of its coordinates — the map the natural basis, the Lamé coefficients and the Christoffel symbols are all derived from.

TensND.unitvec Function

unitvec(CS::_PointwiseSystem, x₀, i) → AbstractTens

Return the i-th unit vector of the normalized basis at point x₀.

TensND.natvec Function

natvec(CS::_PointwiseSystem, x₀, i, var=:cov) → AbstractTens

Return the i-th natural basis vector at point x₀, either covariant (var=:cov, proportional to χᵢ) or contravariant (var=:cont, proportional to 1/χᵢ).

TensND.normalized_basis Function

normalized_basis(CS::_PointwiseSystem, x₀) → AbstractBasis

Return the normalized (orthonormal) basis at point x₀ as a RotatedBasis or CanonicalBasis (when the system is Cartesian at that point).

TensND.natural_basis Function

natural_basis(CS::_PointwiseSystem, x₀) → AbstractBasis

Return the natural (non-normalized) basis at point x₀ as an OrthogonalBasis (i.e. a RotatedBasis scaled by the Lamé coefficients).

TensND.lame_coefficients Function

lame_coefficients(CS::AbstractCoorSystem) → NTuple
lame_coefficients(CS::CoordinateSystem{dim,T,<:PointwiseGeometry}, x₀::AbstractVector) → Vector

Lamé coefficients χᵢ = ‖𝐚ᵢ‖, the norms of the natural basis vectors.

They relate the natural and normalized bases, 𝐞ᵢ = 𝐚ᵢ/χᵢ, and for an orthogonal system give the line element ds² = Σᵢ χᵢ² (dqⁱ)². Symbolic systems return expressions; a PointwiseGeometry evaluates them at a point.

lame_coefficients(coorsys_spherical()) is (r, r sin(θ), 1).

See also Curvilinear differential calculus.

TensND.christoffel Function

christoffel(CS::AbstractCoorSystem) → Array{T,3}
christoffel(CS::CoordinateSystem{dim,T,<:PointwiseGeometry}, x₀::AbstractVector) → Array{T,3}

Christoffel symbols of CS, Γᵏᵢⱼ = ∂ᵢ𝐚ⱼ ⋅ 𝐚ᵏ, symmetric in (i,j).

Storage convention

The array is indexed Γ[i,j,k] = Γᵏᵢⱼ — the contravariant index last. The same convention is used by the Γ_func closure of PointwiseGeometry.

They are what distinguishes a derivative on a curvilinear chart from a plain partial derivative; every one of them vanishes exactly for a Cartesian chart.

See also lame_coefficients, Curvilinear differential calculus.

TensND.nderiv Function

nderiv(CS::AbstractCoorSystem) → Int

Number of independent differentiation directions of CS.

dim for a chart of the whole space, dim-1 for a Submanifold, which differentiates along its surface coordinates only. Every differential operator loops to nderiv(CS), which is what lets a single implementation serve both.

TensND.coorsys_cartesian Function

coorsys_cartesian(coords = symbols("x y z", real = true))

Return the cartesian coordinate system, in which the natural and normalized bases both coincide with the canonical one and all Christoffel symbols vanish.

Examples

The divergence of a general symmetric order-2 field, which in cartesian coordinates reduces to the plain sum of partial derivatives:

julia> Cartesian = coorsys_cartesian() ; 𝐗 = get_coords(Cartesian) ;

julia> ℬ = normalized_basis(Cartesian) ;

julia> 𝛔 = Tens(SymmetricTensor{2,3}((i, j) -> SymFunction("σ$i$j", real = true)(𝐗...))) ;

julia> get_array(divergence(𝛔, Cartesian))
 Derivative(σ11(x, y, z), x) + Derivative(σ21(x, y, z), y) + Derivative(σ31(x, y, z), z)
 Derivative(σ21(x, y, z), x) + Derivative(σ22(x, y, z), y) + Derivative(σ32(x, y, z), z)
 Derivative(σ31(x, y, z), x) + Derivative(σ32(x, y, z), y) + Derivative(σ33(x, y, z), z)

which is divergence(𝛔)ᵢ = ∂ⱼσᵢⱼ.

See also coorsys_polar, coorsys_cylindrical, coorsys_spherical, @set_coorsys.

coorsys_cartesian(PointwiseGeometry, T = Float64)

Numerical Cartesian coordinate system in 3D. Coordinates: (x, y, z).

TensND.coorsys_polar Function

coorsys_polar(coords = (symbols("r", positive = true), symbols("θ", real = true)); canonical = false)

Return the polar coordinate system

Examples

julia> Polar = coorsys_polar() ; r, θ = get_coords(Polar) ; 𝐞ʳ, 𝐞ᶿ = unitvec(Polar) ; ℬᵖ = get_basis(Polar)

julia> f = SymFunction("f", real = true)(r, θ) ;

julia> laplacian(f, Polar)
                               2
                              ∂
                             ───(f(r, θ))
                               2
               ∂             ∂θ
  2            ──(f(r, θ)) + ────────────
 ∂             ∂r                 r
───(f(r, θ)) + ──────────────────────────
  2                        r
∂r

coorsys_polar(PointwiseGeometry, T = Float64)

Numerical polar coordinate system. Coordinates: (r, θ). Lamé coefficients: (1, r).

TensND.coorsys_cylindrical Function

coorsys_cylindrical(coords = (symbols("r", positive = true), symbols("θ", real = true), symbols("z", real = true)); canonical = false)

Return the cylindrical coordinate system

Examples

julia> Cylindrical = coorsys_cylindrical() ; rθz = get_coords(Cylindrical) ; 𝐞ʳ, 𝐞ᶿ, 𝐞ᶻ = unitvec(Cylindrical) ; ℬᶜ = get_basis(Cylindrical)

julia> 𝐯 = Tens(Vec{3}(i -> SymFunction("v$(rθz[i])", real = true)(rθz...)), ℬᶜ) ;

julia> divergence(𝐯, Cylindrical)
                                                  ∂
                                    vr(r, θ, z) + ──(vθ(r, θ, z))
∂                 ∂                               ∂θ
──(vr(r, θ, z)) + ──(vz(r, θ, z)) + ─────────────────────────────
∂r                ∂z                              r

coorsys_cylindrical(PointwiseGeometry, T = Float64)

Numerical cylindrical coordinate system. Coordinates: (r, θ, z). Lamé coefficients: (1, r, 1).

TensND.coorsys_spherical Function

coorsys_spherical(coords = (symbols("θ", real = true), symbols("ϕ", real = true), symbols("r", positive = true)); canonical = false)

Return the spherical coordinate system

Examples

julia> Spherical = coorsys_spherical() ; θ, ϕ, r = get_coords(Spherical) ; 𝐞ᶿ, 𝐞ᵠ, 𝐞ʳ = unitvec(Spherical) ; ℬˢ = get_basis(Spherical)

julia> for σⁱʲ ∈ ("σʳʳ", "σᶿᶿ", "σᵠᵠ") @eval $(Symbol(σⁱʲ)) = SymFunction($σⁱʲ, real = true)($r) end ;

julia> 𝛔 = σʳʳ * 𝐞ʳ ⊗ 𝐞ʳ + σᶿᶿ * 𝐞ᶿ ⊗ 𝐞ᶿ + σᵠᵠ * 𝐞ᵠ ⊗ 𝐞ᵠ ;

julia> div𝛔 = divergence(𝛔, Spherical)
Tens.TensRotated{1, 3, Sym, Vec{3, Sym}}
# data: 3-element Vec{3, Sym}:
                              (-σᵠᵠ(r) + σᶿᶿ(r))*cos(θ)/(r*sin(θ))
                                                                 0
 Derivative(σʳʳ(r), r) + (σʳʳ(r) - σᵠᵠ(r))/r + (σʳʳ(r) - σᶿᶿ(r))/r
# basis: 3×3 Tensor{2, 3, Sym, 9}:
 cos(θ)⋅cos(ϕ)  -sin(ϕ)  sin(θ)⋅cos(ϕ)
 sin(ϕ)⋅cos(θ)   cos(ϕ)  sin(θ)⋅sin(ϕ)
       -sin(θ)        0         cos(θ)
# var: (:cont,)

julia> div𝛔 ⋅ 𝐞ʳ
d            σʳʳ(r) - σᵠᵠ(r)   σʳʳ(r) - σᶿᶿ(r)
──(σʳʳ(r)) + ─────────────── + ───────────────
dr                  r                 r

coorsys_spherical(PointwiseGeometry, T = Float64)

Numerical spherical coordinate system. Coordinates: (θ, ϕ, r) — same convention as coorsys_spherical(). Lamé coefficients: (r, r·sin(θ), 1).

TensND.coorsys_spheroidal Function

coorsys_spheroidal(coords = (symbols("ϕ", real = true),symbols("p", real = true),symbols("q", positive = true),),
                        c = symbols("c", positive = true),tmp_coords = (symbols("p̄ q̄", positive = true)...,),)

Return the spheroidal coordinate system

Examples

julia> Spheroidal = coorsys_spheroidal() ; OM = position_vector(Spheroidal)
Tens.TensCanonical{1, 3, Sym, Vec{3, Sym}}
# data: 3-element Vec{3, Sym}:
 c⋅p̄⋅q̄⋅cos(ϕ)
 c⋅p̄⋅q̄⋅sin(ϕ)
          c⋅p⋅q
# basis: 3×3 Tens.LazyIdentity{3, Sym}:
 1  0  0
 0  1  0
 0  0  1
# var: (:cont,)

julia> laplacian(OM[1]^2, Spheroidal)
2