Geometry

AbstractCartesianGeometry{T} <: AbstractGeometry{T}

Cartesian metrics. Default: CartesianGeometry.

AbstractEllipsoidalGeometry{T} <: AbstractGeometry{T}

Oblate-spheroid metrics. Default: SpheroidGeometry. Coordinate names match the spherical convention, (λ, φ[, h]), with φ the GEODETIC latitude.

A subtype implements semimajor_axis and flattening; the rest — semiminor_axis, eccentricity², the curvature radii, distance, area_element, scale_factors — follow from those two.

AbstractGeometry{T<:AbstractFloat}

Supertype for coordinate metrics (T is the float type).

AbstractSphericalGeometry{T} <: AbstractGeometry{T}

Spherical metrics. Default: SphericalGeometry.

A subtype implements radius; every other method here is written in terms of it.

CartesianGeometry()
CartesianGeometry(T)
CartesianGeometry{T}()

Flat metric in T, of any dimension.

It carries no grid spacing. A cell's extent belongs to the grid's axes, which already state it per cell and per direction; a nominal dx/dy/dz on the geometry would be a second, unchecked copy that could contradict them. The dimension likewise lives on the grid.

PoleRotation(λp, φp)

The frame whose north pole sits at (λp, φp) of the original one — a rotated-pole grid's coordinate change. Apply it with rotate and undo it with unrotate.

SphericalGeometry{T} <: AbstractSphericalGeometry{T}

Default spherical geometry with sphere radius R (meters; default Earth 6.371e6).

SpheroidGeometry(a, f)
SpheroidGeometry()

Oblate spheroid of equatorial radius a and flattening f = (a-b)/a. The no-argument form is WGS 84, a = 6378137.0, f = 1/298.257223563.

distance, area_element, volume_element and scale_factors differ from a sphere. Grid directions are (λ, φ, h), with h the height above the ellipsoid — not an absolute radius, which is what the spherical third direction is.

Rectilinear grids and the index-space connectivity built on them work as they do for a sphere; the samplings are purely angular and so carry over unchanged. The spherical area routines (unstructured_grid's Voronoi areas, the cell areas behind spherical_grid) do not: they are built on spherical excess and 4πR²/n, which are sphere identities, so they stay restricted to AbstractSphericalGeometry rather than silently returning sphere areas for an ellipsoid.

Component names of an ambient-frame vector (as returned by local_tangent_basis).

Convert positional point values to the geometry's own float type T.

area_element(geo::AbstractCartesianGeometry, dx, dy)
area_element(geo::AbstractSphericalGeometry, φ, dλ, dφ)

Local cell area from the cell's own extents: dx·dy on the plane, R²·cosφ·dλ·dφ on the sphere.

area_element(geo::AbstractEllipsoidalGeometry, φ, dλ, dφ)

M(φ)·N(φ)cosφ·dλ·dφ, the exact ellipsoidal surface element — equal to a²(1-e²)cosφ / (1-e²sin²φ)²·dλ·dφ.

as_ntuple(p) -> Tuple

Normalize a point (Tuple, NamedTuple, or AbstractVector of length 1–3) to a plain Tuple.

as_tensor6(τ) -> NTuple{6}

Normalize a symmetric rank-2 tensor to a plain 6-tuple of its independent components, from a Tuple, NamedTuple or AbstractVector — the rank-2 counterpart of as_ntuple.

build_point(S, names::NTuple{N,Symbol}, vals::NTuple{N}) -> S

Assemble a point/vector in representation S from positional vals and their names. Bare NamedTuple takes names; a parameterized NamedTuple type takes its own. Everything else is built from vals alone, covering Tuple, NTuple{N,T}, Vector{T}, SVector{N,T}, MVector{N,T}, and any user type with a positional constructor.

cartesian_to_spherical(geo, xyz) -> NamedTuple{(:λ,:φ,:r)}
cartesian_to_spherical(S, geo, xyz) -> S

Inverse of spherical_to_cartesian: Cartesian (x, y, z) to (λ, φ, r) with λ ∈ (-π, π], φ ∈ [-π/2, π/2], and r the absolute radius from the origin. At the origin λ = φ = 0.

distance(geo::AbstractEllipsoidalGeometry, p1, p2)

2-D (λ, φ): geodesic distance by Vincenty's inverse method, iterated to 1e-12 in the auxiliary longitude, which holds the result to well under a millimetre at Earth scale.

Vincenty's iteration converges slowly for very nearly antipodal point pairs. It is capped, and on reaching the cap the great-circle distance on a sphere of the mean radius is returned instead of a half-converged number.

3-D (λ, φ, h): the chord through Cartesian (ECEF), mirroring the spherical 3-D convention. A lone (λ,) is the shorter arc of the equator.

distance(geo::AbstractGeometry, pt1, pt2)

Distance between two points. Spherical 2D uses great-circle (Haversine); spherical 3D uses the chord through Cartesian, with r the absolute radius from the origin; a lone (λ,) is the shorter arc of its circle.

First eccentricity squared, e² = f(2-f).

flattening(geo) -> T

Flattening f = (a-b)/a. Defaults to the f field.

float_type(geo) -> Type{<:AbstractFloat}

The element type geo computes in. Anything that pairs a geometry with an optional element type should default to this rather than to Float64, or supplying a geometry of one width and no element type silently rebuilds it at the other.

geodetic_to_cartesian(geo, coords) -> (; x, y, z)

(λ, φ, h) — geodetic latitude and height above the ellipsoid — to Earth-centred Cartesian: x = (N(φ)+h)cosφ·cosλ, y = (N(φ)+h)cosφ·sinλ, z = (N(φ)(1-e²)+h)sinφ. The 2-D form takes h = 0, the surface.

jacobian(geo, point) -> T

∏ of scale_factors: the volume element per unit coordinate volume at point.

local_tangent_basis(geo, coords) -> NamedTuple
local_tangent_basis(S, geo, coords) -> NamedTuple of `S`

Coordinate-aligned unit vectors at coords, in the ambient Cartesian frame:

• Cartesian: (; x = ê_x, y = ê_y), 2-component
• Spherical: (; λ = ê_λ, φ = ê_φ), 3-component

The basis vectors themselves are Tuples by default; pass a leading S to get them as SVectors (or any other representation) so they can be used with vector arithmetic directly.

meridional_radius(geo, φ) -> T

M(φ) = a(1-e²) / (1 - e²sin²φ)^{3/2}: the radius of curvature along the meridian. A unit step in latitude covers M(φ).

named_point(geo, vals::NTuple{N}) -> NamedTuple

Build the geometry-appropriate named point from positional values (axis order).

nonuniform_first_derivative(f_m, f_0, f_p, h_m, h_p)

Standard 3-point, 2nd-order-accurate centered finite-difference approximation of the first derivative at the middle node on a possibly nonuniform stencil. f_m, f_0, f_p are the function values at the minus/center/plus nodes, and h_m = x_0 - x_{-}, h_p = x_{+} - x_0 > 0 are the (physical) left/right spacings.

Reduces exactly to the uniform central difference (f_p - f_m)/(2h) when h_m == h_p == h, and is exact for linear and quadratic f for any h_m, h_p.

point_names(geo, Val(N)) -> NTuple{N,Symbol}

Field names for an N-coordinate point in this geometry. Cartesian: (:x,), (:x,:y), (:x,:y,:z), then (:x1, …, :xN). Spherical: (:λ,), (:λ,:φ), (:λ,:φ,:r), then (:λ,:φ,:r,:q4,…).

prime_vertical_radius(geo, φ) -> T

N(φ) = a / √(1 - e²sin²φ): the radius of curvature perpendicular to the meridian. A unit step in longitude covers N(φ)cosφ.

project_to_tangent_plane(geo, center, neighbor) -> NamedTuple
project_to_tangent_plane(S, geo, center, neighbor) -> S

Tangent-plane displacement of neighbor relative to center along the local coordinate basis from local_tangent_basis.

radius(geo) -> T

The sphere radius. The one method a spherical geometry supplies; defaults to the R field.

rotate!(λ, φ, rot) -> (λ, φ)
unrotate!(λ, φ, rot) -> (λ, φ)

Rotate a whole point set in place — the form a sampling's spherical_points output takes. λ and φ are any arrays of matching shape, so this covers a scattered node set and a grid's 2-D coordinate fields alike. Allocates nothing.

rotate(rot, λ, φ) -> (λ′, φ′)

(λ, φ) expressed in the rotated frame. The rotation's own pole maps to φ′ = π/2.

scale_factors(geo, point) -> NTuple

Physical length of a unit coordinate step in each direction at point — see FlowGeometries.Discretization.scale_factors.

Cartesian: 1 in every direction, the metric being the identity. Spherical: (R·cosφ, R) on the surface of radius R, and (r·cosφ, r, 1) where a radius direction is present, r being that point's own radius.

semimajor_axis(geo) -> T

Equatorial radius a. Defaults to the a field.

Semi-minor axis b = a(1-f).

similar_geometry(T, geo) -> AbstractGeometry{T}

geo with its element type changed to T, keeping its shape parameters.

A geometry fixes the width of everything computed against it — coordinates, distances, metric factors — so building a grid at one element type around a geometry at another silently promotes the whole grid back. Anything that takes an element type and a geometry together should pass the geometry through here first.

A geometry already at T is returned unchanged, whatever its type, so a geometry defined outside this package inherits the whole stack without supplying anything. Only an actual change of width needs a method, and only a geometry that wants to support one has to define it.

spherical_excess(a, b, c) -> T

Spherical excess of the triangle spanned by three UNIT vectors, via Van Oosterom & Strackee (1983):

tan(E/2) = |a · (b × c)| / (1 + a·b + b·c + c·a)

One atan and no other transcendental, versus L'Huilier's three great-circle distances (each its own trig) plus four tangents. It takes directions rather than (λ, φ) so a mesh can convert each vertex once — see unit_vector — instead of re-deriving them per triangle. Multiply by R² for an area, which is what triangle_area does.

spherical_to_cartesian(geo, coords) -> NamedTuple{(:x,:y,:z)}
spherical_to_cartesian(S, geo, coords) -> S

(λ, φ) on the reference sphere of radius R, or (λ, φ, r) with r the absolute radius from the origin, mapped to Cartesian coordinates of the same space.

tensor_from_local(geo, τλλ, τφφ, τrr, τλφ, τλr, τφr, λ, φ) -> NamedTuple
tensor_from_local(geo, τ, λ, φ) -> NamedTuple
tensor_from_local(S, geo, args...) -> S

The inverse of tensor_to_local: τ = Rᵀ τ' R, returning the ambient Cartesian components (; xx, yy, zz, xy, xz, yz). Composing the two in either order is the identity.

tensor_to_local(geo, τxx, τyy, τzz, τxy, τxz, τyz, λ, φ) -> NamedTuple
tensor_to_local(geo, τ, λ, φ) -> NamedTuple
tensor_to_local(S, geo, args...) -> S

Rotate a symmetric rank-2 tensor from the ambient Cartesian frame into the local (ê_λ, ê_φ, ê_r) frame at (λ, φ): τ' = R τ Rᵀ, returning (; λλ, φφ, rr, λφ, λr, φr).

The rank-1 counterpart is vector_from_cartesian; between them they cover essentially every field anyone rotates — a stress, a strain rate, a covariance, a subfilter flux. τ may be given as the six independent components or as any 6-component representation, in the order above.

R is the same basis local_tangent_basis defines, so the two cannot drift apart.

triangle_area(geo, p1, p2, p3) -> T

Exact area of the spherical triangle through three (λ, φ) points, each in any accepted representation. R² times spherical_excess.

Convenient rather than fast: it converts all three points every call. A mesh that already holds vertex directions should use spherical_excess on those and scale once.

triangle_area_from_unit_vectors(geo, u1, u2, u3) -> T
triangle_area_from_unit_vectors(R², u1, u2, u3) -> T

triangle_area for vertices already held as unit vectors — R² times spherical_excess, with no coordinate conversion.

A separate name rather than a triangle_area method because the two cannot be told apart by dispatch: a unit vector and a 3-D spherical point (λ, φ, r) are both 3-tuples.

Pass R² itself to walk many triangles: a mesh shares vertices between cells, so the squaring belongs outside the loop alongside the one-per-vertex unit_vector call. That is the form the cell-area and Voronoi paths here use.

unit_vector(T, p) -> NTuple{3,T}

The direction of the (λ, φ) point p on the unit sphere, in any accepted point representation. This is spherical_to_cartesian with the radius divided out, as a bare tuple: the form the spherical-triangle kernels want, computed once per vertex and reused across every triangle sharing it.

unrotate(rot, λ′, φ′) -> (λ, φ)

Inverse of rotate.

vector_from_cartesian(geo, ux, uy, uz, λ, φ) -> NamedTuple{(:λ,:φ,:r)}
vector_from_cartesian(geo, v, λ, φ) -> NamedTuple{(:λ,:φ,:r)}
vector_from_cartesian(S, geo, args...) -> S

Map Cartesian velocity (ux, uy, uz) at (λ, φ) into spherical components (λ = u_λ, φ = u_φ, r = u_r). v may be any 3-component point representation.

vector_to_cartesian(geo, u_λ, u_φ, u_r, λ, φ) -> NamedTuple{(:x,:y,:z)}
vector_to_cartesian(geo, u_λ, u_φ, λ, φ) -> NamedTuple{(:x,:y,:z)}
vector_to_cartesian(S, geo, args...) -> S

Map the components of a vector given in the local spherical basis at (λ, φ) into the Cartesian basis. This transforms VECTOR COMPONENTS at a point, not a position — see spherical_to_cartesian for positions. The 2-component form assumes u_r = 0.

volume_element(geo::AbstractCartesianGeometry, dx, dy, dz)
volume_element(geo::AbstractSphericalGeometry, r, φ, dλ, dφ, dr)

Local cell volume from the cell's own extents. The spherical form uses the LOCAL radius r at this level, not the reference radius, so it is the genuine shell element r²·cosφ·dλ·dφ·dr.

volume_element(geo::AbstractEllipsoidalGeometry, φ, h, dλ, dφ, dh)

(N(φ)+h)cosφ · (M(φ)+h) · dλ·dφ·dh, the geodetic volume element at ellipsoidal height h.

Unlike the spherical form this does not factor into a function of φ times a function of h: both curvature radii are offset by h, so the two directions are coupled.