API Reference

abstract type AbstractSurface{M<:Manifold,C<:CRS} end

A surface of a PointSurface.

struct PointSurface{M,C,T,G} <: AbstractSurface{M,C}

This is a typical representation of a surface via points.

Type Parameters

• M<:Manifold - manifold type
• C<:CRS - coordinate reference system
• T<:AbstractTopology - topology type for surface-local connectivity
• G<:StructVector - storage type for surface elements

struct SurfaceElement{M,C,N,A}

Representation of a point on a <:PointSurface.

struct PointBoundary{M,C} <: Domain{M,C}

A boundary of points.

Fields

• surfaces: Named surfaces forming the boundary

Type Parameters

• M <: Manifold: The manifold type
• C <: CRS: The coordinate reference system

struct PointVolume{M,C,T,V} <: Domain{M,C}

Interior volume points with optional topology.

Type Parameters

• M<:Manifold - manifold type
• C<:CRS - coordinate reference system
• T<:AbstractTopology - topology type for volume-local connectivity
• V<:AbstractVector{Point{M,C}} - storage type (allows GPU arrays)

struct PointCloud{M,C,T} <: Domain{M,C}

A point cloud with optional topology (connectivity).

Type Parameters

• M<:Manifold - manifold type
• C<:CRS - coordinate reference system
• T<:AbstractTopology - topology type for cloud-level connectivity

points(surf::PointSurface)

Return vector of point coordinates for all surface elements.

points(vol::PointVolume)

Return vector of points from volume.

points(boundary::PointBoundary)

Return vector of all points from all surfaces in the boundary.

points(cloud::PointCloud)

Return vector of all points (boundary + volume).

normal(surf::PointSurface)

Return the vector of outward unit normal vectors for all surface elements.

area(surf::PointSurface)

Return the vector of surface areas for all surface elements.

topology(surf::PointSurface)

Return the topology of the surface.

topology(vol::PointVolume)

Return the topology of the volume.

topology(cloud::PointCloud)

Return the topology of the point cloud.

abstract type AbstractTopology{S}

Abstract base type for point cloud topology (connectivity). Type parameter S is the storage format for neighbor indices.

struct NoTopology <: AbstractTopology{Nothing}

Singleton type representing no topology. Default for PointCloud.

mutable struct KNNTopology{S} <: AbstractTopology{S}

k-nearest neighbors topology.

Fields

• neighbors::S - neighbor indices storage
• k::Int - number of neighbors per point

mutable struct RadiusTopology{S,R} <: AbstractTopology{S}

Radius-based topology where neighbors are all points within a given radius.

Fields

• neighbors::S - neighbor indices storage
• radius::R - search radius (scalar or function of position)

set_topology(surf::PointSurface, ::Type{KNNTopology}, k::Int)

Build and return new surface with k-nearest neighbor topology.

set_topology(surf::PointSurface, ::Type{RadiusTopology}, radius)

Build and return new surface with radius-based topology.

set_topology(vol::PointVolume, ::Type{KNNTopology}, k::Int)

Build and return new volume with k-nearest neighbor topology.

set_topology(vol::PointVolume, ::Type{RadiusTopology}, radius)

Build and return new volume with radius-based topology.

set_topology(cloud::PointCloud, ::Type{KNNTopology}, k::Int)

Build and return new cloud with k-nearest neighbor topology.

set_topology(cloud::PointCloud, ::Type{RadiusTopology}, radius)

Build and return new cloud with radius-based topology.

rebuild_topology!(topo::NoTopology, points)

No-op for NoTopology (nothing to rebuild).

rebuild_topology!(topo::KNNTopology, points)

Rebuild k-nearest neighbor topology in place.

rebuild_topology!(topo::RadiusTopology, points)

Rebuild radius-based topology in place.

rebuild_topology!(surf::PointSurface)

Rebuild topology in place using same parameters. No-op if NoTopology.

rebuild_topology!(vol::PointVolume)

Rebuild topology in place using same parameters. No-op if NoTopology.

rebuild_topology!(cloud::PointCloud)

Rebuild topology in place using same parameters. No-op if NoTopology.

neighbors(t::AbstractTopology)

Return the neighbor storage from a topology.

neighbors(t::AbstractTopology, i::Int)

Return neighbors of point i.

neighbors(surf::PointSurface)

Return all neighbor lists from the surface topology. Throws error if no topology.

neighbors(surf::PointSurface, i::Int)

Return neighbors of point i in surface-local indices. Throws error if no topology.

neighbors(vol::PointVolume)

Return all neighbor lists from the volume topology. Throws error if no topology.

neighbors(vol::PointVolume, i::Int)

Return neighbors of point i in volume-local indices. Throws error if no topology.

neighbors(cloud::PointCloud)

Return all neighbor lists from the topology. Throws error if no topology or invalid.

neighbors(cloud::PointCloud, i::Int)

Return neighbors of point i. Throws error if no topology or invalid.

hastopology(surf::PointSurface)

Check if surface has a topology (not NoTopology).

hastopology(vol::PointVolume)

Check if volume has a topology (not NoTopology).

hastopology(cloud::PointCloud)

Check if point cloud has a topology (not NoTopology).

discretize(bnd::PointBoundary, spacing; alg=auto, max_points=nothing)

Generate volume points for the given boundary and return a new PointCloud.

spacing can be either an AbstractSpacing object or a bare Unitful.Length value (which will be wrapped in ConstantSpacing).

Keyword Arguments

• alg: Discretization algorithm (default: SlakKosec() for 3D)
• max_points: Maximum number of volume points to generate. For the Octree algorithm, defaults to an automatic estimate from the spacing integral (∫ 1/h(x)³ dx) when nothing; other algorithms default to 10000000.

Example

mesh = import_mesh("model.stl", u"m")
boundary = PointBoundary(mesh)
octree = TriangleOctree(mesh; min_ratio=1e-6)
cloud = discretize(boundary, 3.0m; alg=SlakKosec(octree), max_points=100_000)

discretize(cloud::PointCloud, spacing; alg=auto, max_points=nothing)

Generate volume points for an existing cloud and return a new PointCloud with the volume populated.

For the Octree algorithm, max_points defaults to an automatic estimate from the spacing integral (∫ 1/h(x)³ dx) when nothing. Other algorithms default to 10000000.

AbstractSpacing

Interface for spacing functions that control node density during discretization.

Subtypes must be callable with a single Point or Vec argument and return a Unitful.Length representing the desired node spacing at that location.

(s::MySpacing)(p::Union{Point, Vec}) -> Unitful.Length

See ConstantSpacing, LogLike, and BoundaryLayerSpacing for concrete implementations.

SlakKosec <: AbstractNodeGenerationAlgorithm

Slak-Kosec algorithm for volume point generation with optional octree acceleration.

The algorithm generates candidate points on spheres around existing points and accepts them if they are inside the domain and sufficiently far from existing points.

Fields

• n::Int - Number of candidate points per sphere (default: 10)
• octree::Union{Nothing,TriangleOctree} - Optional octree for fast isinside queries

Constructors

SlakKosec()                          # Default: n=10, no octree
SlakKosec(20)                        # Custom n, no octree
SlakKosec(octree::TriangleOctree)    # Use octree acceleration with n=10
SlakKosec(20, octree)                # Custom n with octree acceleration

Performance

• Without octree: Uses Green's function for isinside (~50ms per query)
• With octree: Uses spatial indexing (~0.05ms per query, 1000× faster!)

Usage Examples

Standard Usage (Green's function)

using WhatsThePoint

# Load boundary
boundary = PointBoundary("model.stl", u"m")
cloud = PointCloud(boundary)

# Discretize without octree (slow for large domains)
spacing = ConstantSpacing(1.0u"m")
result = discretize(cloud, spacing; alg=SlakKosec(), max_points=10_000)

Octree-Accelerated Usage (Recommended for large domains)

using WhatsThePoint

# Load the mesh once; boundary and octree share it
mesh = import_mesh("model.stl", u"m")
boundary = PointBoundary(mesh)
cloud = PointCloud(boundary)
octree = TriangleOctree(mesh; min_ratio=1e-6, classify_leaves=true)

# Use octree-accelerated discretization (100-1000× faster!)
spacing = ConstantSpacing(1.0u"m")
alg = SlakKosec(octree)  # Pass octree to algorithm
result = discretize(cloud, spacing; alg=alg, max_points=100_000)

References

Šlak J, Kosec G. "On generation of node distributions for meshless PDE discretizations" (2019)

VanDerSandeFornberg <: AbstractNodeGenerationAlgorithm

3D volume discretization algorithm that projects a 2D grid onto the shadow plane and fills the volume layer by layer using sphere packing heights. Requires ConstantSpacing.

See: Van der Sande, K. & Fornberg, B. (2021). SIAM J. Sci. Comput., 43(1).

FornbergFlyer <: AbstractNodeGenerationAlgorithm

2D volume discretization algorithm using a height-field approach projected onto the x-axis. This is the default and only algorithm for 2D boundaries. Requires ConstantSpacing.

See: Fornberg, B. & Flyer, N. (2015). Comput. Math. Appl., 69(7).

Octree <: AbstractNodeGenerationAlgorithm

Spacing-driven volume discretization algorithm.

Note: This is not solution-adaptive (AMR). Refinement is determined a priori by a prescribed spacing function, not by computed solution features.

Algorithm

Uses dual octrees internally:

• Triangle octree: Captures geometry (surfaces, curvature)
• Node octree: Spacing-driven subdivision where h_box ≤ alpha * h_spacing(x)

Points are generated by weighted allocation across octree leaves based on local spacing requirements.

Key Parameters

• min_ratio: Triangle octree resolution (default: auto from mesh complexity)
• node_min_ratio: Node octree resolution (default: auto from spacing if provided)
• alpha: Subdivision aggressiveness, h_box ≤ alpha * h_spacing (default: 2.0, use 1.0 for fine boundary layers)
• placement: :random, :jittered, :lattice, or :bridson (default: :bridson — global graded Poisson-disk, the recommended production sampler)
• bridson_factor: Poisson-disk radius relative to h(x) for :bridson (default: 0.75)
• boundary_oversampling: Oversampling near boundaries (default: 2.0)
• max_growth: Lipschitz cap on the spacing gradient |∇h| (default: 0.0 = off). When > 0, the prescribed spacing is gradient-limited so neighbouring points differ in spacing by no more than this rate — steep variations stay sharp where the geometry forces them but transition smoothly, which RBF-FD stencils need. 0.1–0.2 matches CFD boundary-layer growth ratios of 1.1–1.2. See "Gradient-limited spacing" below.

Gradient-limited spacing (max_growth)

A raw spacing function may vary faster than a meshless stencil can tolerate: two adjacent points with very different target spacings give an asymmetric, ill-conditioned neighbourhood. With max_growth = g > 0 the algorithm replaces the prescribed field h₀(x) with its g-Lipschitz envelope h(x) = minᵧ (h₀(y) + g·‖x − y‖) — the steepest field that is everywhere ≤ h₀ and grows no faster than g. The limiter runs on the node-octree leaves (a multi-source min-plus relaxation over a k-NN graph of leaf centres) and, because limiting can make the field finer than h₀ in a transition band, it then refines any leaf the new field out-resolves and re-limits, to a fixpoint. The sampler, the grid resolution, and the point-count estimate all then read the limited field, so the delivered point distribution grades smoothly. The limiter is a no-op (no refinement) when h₀ is already g-smooth.

Placement modes

:random, :jittered, and :lattice sample each octree leaf independently (point counts allocated by leaf volume / local spacing). :bridson instead runs a single global advancing-front Poisson-disk pass (Bridson 2007, graded to h(x)) seeded from the boundary points: every generated point keeps a distance of at least min(rᵢ, rⱼ) with r = bridson_factor · h(x) from every other point — including the boundary — by construction. The front saturates at the disk-packing density, so max_points acts as a cap, not a target; a warning is emitted if the cap truncates the front before saturation (which would leave unfilled regions).

bridson_factor (default 0.75) sets the disk radius r = bridson_factor·h relative to the local spacing. A saturated graded front packs ≈ 0.46–0.52 points per r³ (measured on this implementation, k = 30 attempts; geometry-dependent — ≈ 1.09/h³ on a convex box, ≈ 1.22/h³ on the non-convex Stanford bunny). At the 0.75 default that is ≈ 1.1–1.2× the prescribed 1/h³ density, i.e. the front slightly over-fills the nominal budget; the automatic max_points estimate carries matching headroom (see _BRIDSON_CAP_HEADROOM) so the inward-advancing front saturates rather than truncating — a truncated front leaves the deep interior (filled last) empty. Use 1.0 for strict d_NN ≥ h Poisson-disk sampling — ≈ 50% fewer points than 1/h³, which starts repel's spacing-equilibrium force in its attractive branch and degrades rather than polishes the seeding (measured on the cavity gate).

Examples

# Automatic (recommended)
alg = Octree(mesh; spacing, alpha=1.0)
cloud = discretize(boundary, spacing; alg, max_points=100_000)

# Manual geometry resolution
alg = Octree(mesh; min_ratio=1e-3, spacing, alpha=1.0)

ConstantSpacing{L<:Unitful.Length} <: AbstractSpacing

Constant node spacing.

BoundaryLayerSpacing <: VariableSpacing

Smooth spacing transition from fine spacing at the boundary to coarse spacing in the bulk.

Uses physical boundary layer intuition with clear parameters:

• at_wall: Spacing at the boundary surface (fine)
• bulk: Spacing far from boundaries (coarse)
• layer_thickness: Distance over which transition occurs

Example

# Fine 0.5m spacing at walls, coarse 10m in bulk, 8m boundary layer
spacing = BoundaryLayerSpacing(boundary, at_wall=0.5m, bulk=10m, layer_thickness=8m)

Internally uses sigmoid: h(d) = at_wall + (bulk - at_wall) * σ(d) where σ(d) = 1 / (1 + exp(-(d - δ/2) / (δ/6))) and δ = layer_thickness.

LogLike <: VariableSpacing

Node spacing based on a log-like function of the distance to nearest boundary $x/(x+a)$ where $x$ is the distance to the nearest boundary and $a$ is a parameter to control the growth rate as $a = 1 - (g - 1)$ where $g$ is the conventional growth rate parameter.

compute_normals(surf::PointSurface{𝔼{N},C}; k::Int=5) where {N,C<:CRS}

Estimate the normals of a set of points that form a surface. Uses the PCA approach from "Surface Reconstruction from Unorganized Points" - Hoppe (1992).

Requires Euclidean manifold (𝔼{2} or 𝔼{3}). This function assumes flat space geometry.

compute_normals(search_method::KNearestSearch, surf::PointSurface{𝔼{N},C}) where {N,C<:CRS}

Estimate the normals of a set of points that form a surface. Uses the PCA approach from "Surface Reconstruction from Unorganized Points" - Hoppe (1992).

Requires Euclidean manifold (𝔼{2} or 𝔼{3}). This function assumes flat space geometry.

orient_normals!(search_method::KNearestSearch, normals::AbstractVector{<:AbstractVector}, points)

Correct the orientation of normals on a surface as the compute_normals function does not guarantee if the normal is inward or outward facing. Uses the approach from "Surface Reconstruction from Unorganized Points" - Hoppe (1992).

orient_normals!(normals::AbstractVector{<:AbstractVector}, points::AbstractVector{<:Point{𝔼{N}}}; k::Int=5) where {N}

Correct the orientation of normals on a surface as the compute_normals function does not guarantee if the normal is inward or outward facing. Uses the approach from "Surface Reconstruction from Unorganized Points" - Hoppe (1992).

Requires Euclidean manifold (𝔼{2} or 𝔼{3}). This function uses Euclidean dot products for orientation consistency.

update_normals!(surf::PointSurface{𝔼{N},C}; k::Int=5) where {N,C<:CRS}

Update the normals of the boundary of a surf. This is necessary whenever the points change for any reason.

Requires Euclidean manifold (𝔼{2} or 𝔼{3}). This function assumes flat space geometry.

split_surface!(cloud, angle; k=10)
split_surface!(cloud, target_surf, angle; k=10)

Split a surface into sub-surfaces based on normal angle discontinuities. Builds a k-nearest neighbor graph, removes edges where adjacent normals differ by more than angle, and labels each connected component as a separate named surface.

When called on a cloud/boundary with a single surface, that surface is split automatically. When multiple surfaces exist, specify target_surf by name.

combine_surfaces!(boundary::PointBoundary, surfs...)

Merge multiple named surfaces into one. The first name is kept and subsequent surfaces are merged into it. All original surfaces are removed and replaced by the combined surface.

ShadowPoints(Δ, order=1)
ShadowPoints(Δ::Number, order)

Shadow point configuration for generating virtual points offset inward from the boundary. Δ is the offset distance (constant or a function of position). order is the derivative order for Hermite-type boundary condition enforcement.

generate_shadows(points, normals, shadow::ShadowPoints)
generate_shadows(surf::PointSurface, shadow::ShadowPoints)
generate_shadows(cloud::PointCloud, shadow::ShadowPoints)

Generate shadow points offset inward from the boundary along the normal direction by the distance specified in shadow. Returns a vector of Point objects.

Fast interior/exterior test using octree spatial index.

isinside(testpoint::Point{𝔼{2}}, pts::AbstractVector{<:Point{𝔼{2}}}) -> Bool
isinside(testpoint::Point{𝔼{N}}, cloud::Union{PointCloud, PointBoundary}) -> Bool

Test whether testpoint lies inside the closed domain defined by the boundary points.

For 2D, uses the winding number algorithm — pts must be ordered sequentially around the polygon boundary (clockwise or counter-clockwise). An ArgumentError is thrown if the points do not form a valid ordered polygon.

For 3D, uses a Green's function approach over the boundary surfaces.

Octree spatial index for triangle mesh queries. Accelerates isinside(), signed distance, etc.

Check if triangle faces are consistently oriented (manifold orientation test).

emptyspace(testpoint, points)

Check if a point occupies empty space within a certain tolerance.

boundingbox(pts::AbstractVector{<:Point})

Compute the axis-aligned bounding box of a collection of points.

centroid(pts::AbstractVector{<:Point})

Compute the centroid (geometric center) of a collection of points.

repel(cloud::PointCloud, spacing; kwargs...) -> PointCloud

Optimize the point distribution via node repulsion. Only volume points move; boundary points form a fixed wall. Points pushed outside the domain are discarded (isinside filter). Returns a new cloud with NoTopology.

Each iteration rebuilds the k-NN graph from the current positions, computes the repel force on every point, and moves it by the adaptive step α_i = clamp(1/|F_i|, α_min, α_max) scaled by the local spacing and capped at one spacing unit. Convergence is the force norm max_i(|F_i|·s_i), which vanishes at equilibrium.

Keywords

• force_model = ClippedSpacingForce(β): force law, any RepelForceModel; the default is repulsive below r = s and zero beyond, so a cloud that already satisfies the Poisson-disk criterion is preserved rather than re-packed. β = 0.2 feeds the default and is ignored when force_model is passed explicitly.
• α, α_min: step-size bounds. Defaults: α = 0.05·min(spacing), α_min = α/100.
• k = 21: neighborhood size.
• max_iters = 1000, tol = 1e-6: iteration and convergence limits.
• cv_target = 0.0: quality-based stop — end the relaxation once the movable points' d_NN/s coefficient of variation drops to this value (read off the sweep's nearest-neighbor data, no extra cost). The natural setting is the raw quality of the direct generation pipeline (≈ 0.07 on the cavity): relaxing past the quality a re-seed would give is wasted budget. The stop returns the configuration the measurement describes (the pre-sweep snapshot), so a cloud already at target comes back unchanged. Off when 0.
• stall_after = 50: stop when that same CV has not improved by ≥0.1 % for this many consecutive iterations. The force residual of a saturated repulsion-only packing plateaus at a nonzero value instead of reaching tol, so cv_target/stall_after are the practical stops for the default force; CV keeps creeping down for hundreds of iterations, making stall_after the backstop (on by default so default runs terminate instead of burning max_iters) and cv_target the primary. Pass 0 to disable and rely on tol/max_iters alone.
• rebuild_every = 1: iterations between k-NN graph rebuilds (larger = cheaper, staler).
• kick_after = 0: if the closest pair freezes at the same r/s for this many iterations (a balanced standoff), kick one point by 0.1·s in a random direction to break the symmetry. Off when 0; 10–20 is reasonable.
• cull_ratio = 0.0: after relaxation, drop near-duplicates closer than cull_ratio·spacing to a kept point. A safety net — a healthy relaxation leaves nothing to cull, so a @warn is emitted whenever it fires.
• convergence: pass an empty float vector (e.g. Float64[]) to collect the per-iteration force norm; entries are computed in the cloud's machine type and convert on insertion.
• trace: pass a NamedTuple[] to record the closest pair each iteration (global boundary-then-volume indices, measured on that iteration's snapshot).

repel(cloud::PointCloud, spacing, octree::TriangleOctree; kwargs...) -> PointCloud

Node repulsion with boundary projection: all points move, boundary points are re-projected onto the mesh surface every iteration, and volume points that escape the domain bounce back — or stick to the surface (see deposit_ratio). Step size, convergence, and the shared keywords are as in the method without octree.

The returned boundary is a single surface named :boundary (use split_surface! to re-establish surface distinctions); topology is NoTopology.

Additional keywords

• deposit_ratio = 0.0: when > 0, an escaped volume point is deposited — projected onto the nearest triangle and converted into a boundary point, accepted only if no boundary point already lies within deposit_ratio·spacing of the landing site. Surface sampling then emerges from volume containment instead of the mesh tessellation; the acceptance test keeps the deposited density self-limiting (conversion is one-way). Deposited points carry the landing triangle's normal and spacing² as area. 0.5–0.7 is reasonable.
• cull_ratio = 0.0: as in the volume-only method; here it also targets boundary pairs that the deterministic projection parks on a shared edge/vertex.

RepelForceModel

Abstract type for node-repulsion force laws used by repel.

A concrete subtype M <: RepelForceModel must implement compute_force(m::M, u::Real).

ClippedSpacingForce(β=0.2, u0=1.0)

Repulsion-only force law: F(u) = (u0² − u²) / (u² + β)² for u < u0, zero beyond. The compact support makes any configuration whose pairwise distances all exceed u0·s an exact equilibrium — the Poisson-disk property — so an already-blue-noise cloud is preserved rather than pulled toward a different packing.

This is SpacingEquilibriumForce with the attractive branch removed (identical for u < u0 when u0 = 1). The attractive branch acts as a cohesion force whose preferred bond length s is unreachable at the prescribed density 1/s³: the cloud condenses into locally denser clusters plus voids (measured on the cavity gate: spacing CV and coordination rise, separation falls, starting from an already-good cloud, at a rate proportional to the step size — an instability, not noise). Clipping the force at its root removes the mechanism; the same 300 iterations then improve a constructed blue-noise cloud (CV 0.072 → 0.044) instead of degrading it.

u0 sets the support radius (and root) in units of the local spacing s.

The force residual of a saturated repulsion-only packing plateaus at a small nonzero value instead of vanishing (a frustrated glass), so combine with the stall_after stopping criterion of repel rather than relying on tol alone.

InverseDistanceForce(β=0.2)

Original Miotti (2023) force law F(u) = 1 / (u² + β)². Purely repulsive and monotonically decreasing. β > 0 softens the force near u = 0 and keeps it finite. Has no root, so equilibrium is reached only through damping (α).

SpacingEquilibriumForce(β=0.2)

Force law F(u) = (1 − u²) / (u² + β)² with a zero at u = 1. Repulsive for u < 1 (points closer than the target spacing push apart), attractive for u > 1 (points farther than the target pull together), and zero at the target spacing itself. β > 0 softens the amplitude near u = 0.

Shares the β-softened denominator with InverseDistanceForce; the numerator (1 − u²) introduces the equilibrium at u = 1 without changing the small-u behavior.

StrongSpacingForce(β=0.2, γ=3)

Force law F(u) = (1 − u²) / (u² + β)^γ with a zero at u = 1 and a configurable singularity strength γ. Like SpacingEquilibriumForce but with a stronger repulsive core: at small u the force scales as u^(-2γ) instead of u^(-4). This breaks balanced standoffs where neighbor forces cancel the weaker default core.

γ = 2 recovers SpacingEquilibriumForce. γ = 3 (default) is strong enough to break typical standoffs in a few iterations. Higher values increase the repulsive kick at close range but may require more iterations to settle; the displacement cap in repel prevents runaway.

compute_force(model::RepelForceModel, u::Real) -> Real

Evaluate the force magnitude at normalized separation u = r / s, where r is the distance between two points and s is the local target spacing.

The returned scalar multiplies the unit vector (xᵢ − xⱼ) / r in the repulsion step, so positive values push xᵢ away from xⱼ and negative values pull it toward xⱼ. Concrete subtypes of RepelForceModel must implement a method for this function.

metrics(cloud::PointCloud; k=20)

Compute distance statistics (mean, std, max, min) to the k nearest neighbors for all points in the cloud. Useful for assessing point distribution quality before and after repulsion.

Also reports the global separation and fill distances and their ratio, the quasi-uniformity quality measure most relevant to meshless stencil conditioning:

• separation — the smallest nearest-neighbor distance anywhere in the cloud. Small values signal near-coincident points (the source of singular RBF-FD stencils).
• fill — the largest nearest-neighbor distance (a proxy for the worst covering gap).
• mesh_ratio — fill / separation (≥ 1). Closer to 1 means a more uniform, blue-noise-like cloud; large values indicate clustering and voids coexisting.

Returns a NamedTuple with fields avg, std, max, min, separation, fill, mesh_ratio, and k.

spacing_metrics(cloud::PointCloud, spacing::AbstractSpacing; k=20)

Measure how closely the point distribution matches the target spacing function.

For each point xᵢ, the local actual spacing is estimated as the mean distance to its k nearest neighbors (self excluded). The per-point relative error is

errorᵢ = |r̄ᵢ − s(xᵢ)| / s(xᵢ)

Returns a NamedTuple (max_error, mean_error, std_error, k). Use before and after repel (or any placement step) to quantify spacing preservation.

spacing_fidelity_metrics(cloud::PointCloud, spacing::AbstractSpacing; k=30, coord_radius=1.4)

Per-point spacing fidelity: how well each point's nearest-neighbor distance matches the prescribed spacing h(x).

Computes d_NN(i) / h(x_i) for every point and returns:

• mean_dnn_h — mean of the distribution (ideal ≈ 0.74 for 3D blue-noise)
• cv — coefficient of variation std / mean (lower = more uniform)
• p05, p50, p95 — percentiles of d_NN/h (tight spread = good)
• coordination — mean count of neighbors within coord_radius · h (ideal ≈ 12–14 for 3D blue-noise packing)
• k, coord_radius

suggest_spacing(mesh; n_points=nothing, bridson_factor=0.75, verbose=true)
suggest_spacing(boundary; ...)
suggest_spacing("model.stl", u"mm"; ...)

Quick geometry probe that recommends a baseline node spacing — the "step 0" before discretize. Reports the domain extent, enclosed volume, and three spacing landmarks, and (with verbose=true) prints a short summary.

The recommendation is driven by the shortest bounding-box axis L_min, because that axis sets how coarse a Poisson-disk fill can be before its interior collapses. With the bridson disk radius r = bridson_factor·h, the interior along L_min has width L_min − 2r, so it is empty once h ≥ L_min/(2·bridson_factor) — the reported h_ceiling. Returned landmarks:

• h_ceiling — coarsest spacing that still hosts any interior; spacings at or above it yield an empty bridson cloud. Stay well below this.
• h_baseline — recommended starting point: ≈10 points across the shortest axis (or, when n_points is given, cbrt(volume / n_points) capped to stay fillable). Good enough to run a first simulation, then refine where needed.
• h_fine — h_baseline/2, a second rung for resolving features.

n_baseline/n_fine are rough volume-point counts (volume / h³).

Returns a NamedTuple with extent, min_extent, max_extent, diagonal, volume, n_triangles, bridson_factor, h_ceiling, h_baseline, h_fine, n_baseline, and n_fine (all spacings/lengths carry units).

Example

mesh = import_mesh("bunny.stl", u"m")
g = suggest_spacing(mesh)
cloud = discretize(PointBoundary(mesh), g.h_baseline; alg=Octree(mesh))

geometry_info(filepath, filepaths...; verbose=true) -> Vector{NamedTuple}

Inspect the raw bounding box of one or more mesh files before constructing anything — the "what do these numbers mean?" probe for files without unit metadata (STL coordinates are just numbers; GeoIO.jl assigns meters by default). Returns one (file, min, max, extent) named tuple per file, with values as raw (unitless) coordinate tuples. With verbose=true prints each bounding box and, for multiple files, their union — useful when parts must fit together.

Once you know the unit, pass it to import_mesh or PointBoundary — the raw numbers are reinterpreted in that unit.

Example

geometry_info("intake.stl", "exhaust.stl")
# ─── intake.stl ───
#   min:    (0.0, 0.0, 0.0)
#   max:    (120.5, 87.3, 42.0)
#   extent: (120.5, 87.3, 42.0)
# ─── exhaust.stl ───
#   ...
# ─── Union ───
#   ...

sample_surface(mesh::SimpleMesh, spacing; factor=0.75, max_points=10_000_000,
               stall_limit=2000) -> PointSurface

Sample the surface of mesh with a graded Poisson-disk distribution: blue-noise points with pairwise separation at least min(rᵢ, rⱼ) where r = factor·spacing(x), by construction. An alternative to the face-center sampling of PointBoundary(mesh) whose point density follows the prescribed spacing instead of the mesh tessellation.

Each sample carries its parent triangle's normal; point areas preserve the total mesh surface area, distributed proportionally to r² (equal shares for constant spacing).

Dart throwing runs until saturation — stall_limit consecutive rejections — or until max_points is reached (with a warning, since a truncated pass leaves the surface under-sampled).

import_mesh(filepath, unit::Unitful.Units) -> SimpleMesh

Load a surface mesh from a file (STL, OBJ, or any format supported by GeoIO.jl) and reinterpret its raw coordinates in unit: a stored 46 becomes 46 mm with unit = u"mm" — no conversion happens, because mesh files carry no unit metadata (GeoIO's default of meters is discarded). Topology and coordinate machine type (Float32 stays Float32) are preserved.

This is the single gateway for file geometry: the returned mesh feeds PointBoundary, TriangleOctree, and Octree. Use geometry_info to probe the raw extents when unsure of the unit.

import_surface(filepath, unit::Unitful.Units)

Load a surface mesh from a file (STL, OBJ, or any format supported by GeoIO.jl). Returns a tuple of (points, normals, areas, mesh) where points are face centers. Coordinates are reinterpreted in unit (see import_mesh); areas therefore carry unit^2.

save(filename::String, cloud::PointCloud; format=:jld2)

Save a point cloud to a file.

• format=:jld2 (default): Serialize via FileIO.jl.
• format=:vtk: Export to VTK format with boundary and volume points, normals, and a point type indicator (1 = boundary, 2 = volume).

save(filename::String, boundary::PointBoundary; format=:jld2)

Save a boundary to a file.

• format=:jld2 (default): Serialize via FileIO.jl.
• format=:vtk: Export to VTK format with boundary points and normals.

save(filename::String, surf::PointSurface; format=:jld2)

Save a surface to a file.

• format=:jld2 (default): Serialize via FileIO.jl.
• format=:vtk: Export to VTK format with surface points, normals, and areas.

export_vtk(filename, cloud::PointCloud; fields=(), verbose=false)

Write cloud to a ParaView-ready .vtu (one VTK_VERTEX cell per point). Open it in ParaView and set Representation to Point Gaussian (or Points).

Always-attached point data:

• point_type — 1 = boundary, 2 = volume (colour to separate wall from bulk).
• surface_id — 1..N in names(cloud) order, 0 for volume (colour by named surface). Pass verbose=true to print the integer→name legend.
• normals — boundary normals (zero on volume points).

fields attaches solution data, so a .vtu can be re-exported after solving and viewed like any CAE result. It is an iterable of name => values pairs (e.g. a Dict or a tuple of pairs); values may be scalars or per-point vectors and must be ordered like points(cloud) — boundary points first, then volume (the natural global DOF order). Units are stripped automatically.

Examples

export_vtk("cloud", cloud)                                 # geometry only
export_vtk("sol", cloud; fields = ("T" => temp, "U" => velocity))

AbstractOctree{E,T} = AbstractSpatialTree{3,E,T}

Convenience alias for 3D spatial trees (octrees).

AbstractSpatialTree{N,E,T}

Abstract type for N-dimensional spatial trees.

Type Parameters

• N::Int: Spatial dimensionality (2 for quadtree, 3 for octree)
• E: Element type stored in tree (e.g., Int for indices)
• T<:Real: Coordinate numeric type (Float64, Float32, etc.)

Interface Requirements

Trees must implement:

• find_leaf(tree, point) - Locate leaf containing point
• bounding_box(tree) - Get overall bounds
• num_elements(tree) - Total elements stored

Optional:

• find_neighbors(tree, box_idx, direction) - Neighbor queries
• balance!(tree) - Enforce refinement constraints

AndCriterion{T<:Tuple} <: SubdivisionCriterion

Combine multiple criteria - all must be satisfied for subdivision.

Fields

• criteria::T: Tuple of criteria to combine (parametrized for type stability)

MaxElementsCriterion <: SubdivisionCriterion

Subdivide box if number of elements exceeds threshold.

Fields

• max_elements::Int: Maximum elements before subdivision

MeshPseudonormals{T}

Angle-weighted pseudonormals (Bærentzen & Aanæs 2005) for exact sign determination in signed-distance queries. For a watertight, consistently outward-oriented triangle mesh, sign(dot(p - cp, n_feature)) is provably correct for every query point p, where cp is the closest point on the mesh and n_feature is the normal of the closest feature: the face normal on a face interior, the sum of the two incident face normals on an edge, and the incidence-angle-weighted sum of face normals at a vertex. (Only the sign of the dot product is used, so the edge/vertex sums are stored unnormalized.)

This replaces the previous distance-weighted "local sign vote" over the triangles of one leaf, which was heuristic — it could come back ambiguous (mapped to EXTERIOR by isinside) and could mix votes from the two sides of a thin sheet — and slower (a second find_leaf plus a closest-point computation per vote triangle on every query).

Edges and vertices are keyed by exact coordinates, so triangle-soup input (e.g. binary STL with duplicated vertices) needs no topology cleanup.

SizeCriterion{T<:Real} <: SubdivisionCriterion

Subdivide box if size exceeds threshold.

Fields

• h_min::T: Minimum box size (stop subdividing when reached)

SpacingCriterion{T<:Real, S} <: SubdivisionCriterion

Octree subdivision criterion based on local spacing requirements.

Subdivides boxes where h_box > alpha * h_spacing(center), ensuring the octree resolution is fine enough to properly represent the spacing function.

Fields

• spacing::S: Spacing function object
• alpha::T: Subdivision aggressiveness factor
• absolute_min::T: Absolute minimum box size (prevents infinite subdivision)

Algorithm

For each octree box:

1. Query h_local = spacing(box_center)
2. If h_box > alpha * h_local, subdivide
3. Stop if h_box ≤ absolute_min

Smaller alpha values create finer octrees (more aggressive subdivision).

SpatialOctree{E,T<:Real} <: AbstractOctree{E,T}

Concrete octree implementation using integer coordinate system for efficient neighbor finding.

Uses (i,j,k,N) coordinate system where:

• (i,j,k) are integer coordinates at refinement level N
• Box center = origin + (2*[i,j,k] + 1) * (root_size / N) / 2
• Enables O(1) neighbor calculation from coordinates

Type Parameters

• E: Type of elements stored (e.g., Int for triangle IDs)
• T: Coordinate numeric type (e.g., Float64, Float32)

Fields

• parent::Vector{Int}: Parent box index for each node (0 = root has no parent)
• children::Vector{SVector{8,Int}}: 8 child indices per box (0 = no child)
• coords::Vector{SVector{4,Int}}: (i,j,k,N) coordinates per box where N is refinement level
• origin::SVector{3,T}: Spatial origin of root box
• root_size::T: Size of root box (assumed cubic)
• element_lists::Vector{Vector{E}}: Elements in each box
• num_boxes::Ref{Int}: Current number of boxes (mutable counter)

Interface Implementation

Implements AbstractSpatialTree interface:

• find_leaf(tree, point) - O(log n) point location
• bounding_box(tree) - Root box bounds
• num_elements(tree) - Total stored elements
• find_neighbor(tree, box, dir) - 6-directional neighbor finding
• balance_octree!(tree) - 2:1 refinement constraint

Example

using StaticArrays

origin = SVector(0.0, 0.0, 0.0)
octree = SpatialOctree{Int,Float64}(origin, 10.0)

# Subdivide root
subdivide!(octree, 1)

# Find leaf containing point
point = SVector(2.0, 2.0, 2.0)
leaf = find_leaf(octree, point)

SpatialOctree{E,T}(origin::SVector{3,T}, size::T; initial_capacity=1000)

Create empty octree with root node.

Arguments

• origin: Minimum corner of root box
• size: Edge length of root box (assumed cubic)
• initial_capacity: Initial array capacity (will grow as needed)

Example

origin = SVector(0.0, 0.0, 0.0)
octree = SpatialOctree{Int,Float64}(origin, 10.0)

SubdivisionCriterion

Abstract type for subdivision decision logic.

Allows pluggable subdivision criteria via dispatch.

Built-in Criteria

• MaxElementsCriterion(max_elements) - Subdivide if too many elements
• SizeCriterion(h_min) - Subdivide until box small enough
• AndCriterion(criteria...) - All criteria must be satisfied

Example

criterion = AndCriterion((
    MaxElementsCriterion(50),
    SizeCriterion(0.1)
))

Subdivide based on vertex density within box bounds.

_ClampedSpacing(inner, hmax) <: AbstractSpacing

Spacing that caps inner at hmax: min(inner(p), hmax). Used by _guard_coarse_spacing to make an everywhere-too-coarse spacing fillable without discarding the user's variation where it is already fine.

_LeafSpacing(node_tree, field, fallback, len_unit) <: AbstractSpacing

Spacing that reads a precomputed per-leaf field (the gradient-limited envelope) by find_leaf lookup, falling back to the original spacing outside the tree. Lets the sampler, grid sizing, and point-count estimate all see the limited field through the usual _spacing_value / call interface.

filter(f::Function, vol::PointVolume)

Return new PointVolume with only points satisfying predicate f. Topology is stripped since point indices change.

isvalid(t::AbstractTopology)

Check if topology is valid. With immutable design, topology is always valid if it exists.

_allocate_counts_by_volume(volumes, total_count; ensure_one=false)

Proportionally allocate counts to volumes using probabilistic rounding. Fractional parts become selection probabilities, reducing clustering.

_apply_gradient_limit(node_tree, classification, spacing, alg, tri_octree)
    -> (node_tree, classification, spacing_used)

Gradient-limit the spacing on the node-octree leaves with alg.max_growth, and refine any leaf the limited field now out-resolves (box_size > alpha·h), re-limiting to a fixpoint. Returns the (possibly further-subdivided) tree, the refreshed classification, and a _LeafSpacing that serves the limited field. A no-op tree-wise when h₀ is already g-smooth.

_auto_min_ratio(::Type{T}, mesh) where {T}

Default triangle octree resolution: 1 / (4 * cbrt(n_triangles)).

Factor of 4 (vs. 2) ensures accurate geometry in high-curvature regions. Override with explicit min_ratio parameter if needed.

_box_corners(bbox_min, bbox_max) -> Tuple

Generate 8 corner points of a 3D bounding box.

_bridson_h_min(node_tree, classification, spacing) -> T

Minimum spacing over non-exterior leaf centers — sets the background grid resolution for graded spacing. T is the node tree's coordinate type.

_bridson_inside(c, node_tree, classification, tri_octree) -> Bool

Domain test for a Bridson candidate, mirroring the trust rules of the per-leaf path: interior node-octree leaves are accepted outright, boundary leaves get an exact isinside check, exterior leaves (and anything outside the mesh bbox) are rejected.

_bridson_separated(grid, pts, rs, c, r_c) -> Bool

true when candidate c (local disk radius r_c) is at least min(r_c, r_q) away from every accepted point q. Any rejecting point lies within r_c of c, so scanning buckets within r_c suffices.

_build_knn_neighbors(points, k::Int) -> Vector{Vector{Int}}

Build k-nearest neighbor adjacency list from points.

_build_radius_neighbors(points, radius) -> Vector{Vector{Int}}

Build radius-based adjacency list from points.

_classify_leaf_conservative(tree, leaf_idx, geometry_query, tol_rel, tol_abs) -> Int8

Conservative 9-point classification (center + 8 corners).

_classify_point_octree(point, octree; tol=0) -> Int8

Shared classification of a point against a TriangleOctree. Returns LEAF_INTERIOR, LEAF_EXTERIOR, or LEAF_BOUNDARY. Uses the cached leaf classification for fast INTERIOR/EXTERIOR dispatch; only BOUNDARY probes fall through to the full signed-distance computation.

tol expands the mesh bounding box for the exterior fast-path (set to 0 for exact bbox checks, or a positive tolerance for conservative classification).

_closest_pair(nn_dist, nn_id, spacings, n_fixed) -> NamedTuple

Closest pair read off the sweep's per-point nearest-neighbor data (no extra search). Returns (; r, s, r_over_s, idx_a, idx_b) in snapshot-global indices. A frozen r_over_s across iterations indicates a balanced standoff; an oscillating one, an overshoot limit-cycle.

Signed distance to the mesh: distance to the closest point, signed by the angle-weighted pseudonormal of the closest feature (exact for watertight, consistently outward-oriented meshes — Bærentzen & Aanæs 2005). Returns 0 only for points exactly on the surface (or at a degenerate fold where the feature pseudonormal vanishes).

_constrain_octree(id, xi, x_proposed, is_bnd, escaped, tri_indices,
                  octree, offset_dist, len_unit) -> Point

Wall rule for one point: boundary points are re-projected onto the mesh (falling back to projecting their previous position), volume points keep the proposed position while it stays inside, and escapees revert — flagged in escaped as deposition candidates.

_cull(pts, spacing, ratio) -> BitVector

Near-duplicate keep-mask plus the defect warning: the cull is a safety net, so any non-zero removal is surfaced.

_deposit_escaped!(p, tree, kq, escaped, is_bnd, tri_indices, octree, spacing,
                  deposit_ratio, offset_dist, len_unit) -> n_deposited

One deposition pass: each escaped volume point is projected onto its nearest triangle and converted to a boundary point, unless another boundary point already sits within deposit_ratio·spacing of the landing site. tree is the sweep's snapshot kd-tree and kq the neighbor count to inspect. Serial on purpose — earlier deposits must be visible to later candidates, because conversion is one-way and a parallel pass would over-deposit when a whole layer escapes in one iteration.

_dnn_cv(nn_dist, spacings, n_fixed) -> Real

Coefficient of variation of d_NN/s over the movable points, read off the sweep's per-point nearest-neighbor data. The quality monitor behind stall_after: it tracks the gate's binding spacing-CV metric for free. Computed in the promoted distance/spacing type.

_estimate_volume_points(node_tree, classification, spacing) -> Int

Estimate the auto max_points cap for the Octree algorithm (used when the caller leaves max_points unset) as ⌈1.5 · ∑ box_volume/h(x)³⌉ over non-exterior leaves — the discrete spacing integral ∫ 1/h³ dx with headroom for the super-1/h³ saturated Poisson-disk packing density (see _estimate_volume_points). It is a non-truncating ceiling, not a target.

_extract_min_spacing(spacing)

Extract minimum spacing value from spacing object (field access, no sampling).

Pseudonormal of the feature (face / edge / vertex) the closest point lies on. Falls back to the face normal if the feature key is missing — cannot happen for features built from this mesh's own triangles, but keeps the query total.

_generate_bridson(node_tree, classification, tri_octree, spacing, seeds,
                  max_points; factor=0.75, k_attempts=30)
    -> Vector{SVector{3,T}}

Graded Bridson Poisson-disk sampling of the domain volume with disk radius factor · h(x) (see the Octree docstring for the choice of factor). seeds (boundary points) initialize the advancing front and occupy the background grid so volume points keep their distance from the wall, but are not returned. The front runs until saturation or until max_points volume points exist; truncation warns because it leaves the far side of the front unfilled.

_generate_points_in_box(bbox_min, bbox_max, n, placement)

Generate n points: :random, :jittered (stratified), or :lattice (grid).

_gradient_limit_field(node_tree, leaves, h0_field, g; k, tol, max_sweeps)
    -> Vector{T}

Box-indexed g-Lipschitz envelope of h0_field, restricted to leaves. Multi-source min-plus relaxation over a k-NN graph of leaf centres: every leaf is a source carrying its own h₀, and h[i] ← min(h[i], min_j h[j] + g·dᵢⱼ) is swept to a fixpoint. The k-NN edges approximate Euclidean distance through multi-hop paths, so the result converges to the Euclidean envelope as the leaf sampling refines. Propagation crosses only the leaf graph, so thin exterior gaps between interior regions are bridged only if their centres are k-NN close (acceptable; a geodesic limiter would be exact but is not needed here).

_grid_insert!(grid, p)

Insert the point with the next sequential index into its bucket. Points must be inserted in index order (1, 2, …) — the link vector grows by one per call.

_guard_coarse_spacing(spacing, tri_octree, bridson_factor) -> spacing

Bridson safety net. When the finest prescribed spacing over the domain is at or above the Poisson-disk ceiling L_min/(2·bridson_factor) — i.e. a saturated front would leave the interior empty — emit a loud @warn and return a _ClampedSpacing capped at a usable baseline (L_min/10) so generation still yields a cloud. Otherwise returns spacing unchanged (the request is viable and is respected).

_leaf_class_from_signed_distance(sd, tol) -> Int8

Convert signed distance to leaf classification.

Returns:

• LEAF_BOUNDARY if |sd| ≤ tol (on or near surface)
• LEAF_INTERIOR if sd < -tol (inside)
• LEAF_EXTERIOR if sd > tol (outside)

_leaf_spacing_field(::Type{T}, node_tree, leaves, spacing) -> Vector{T}

Box-indexed vector of the prescribed spacing evaluated at each leaf centre (zero for non-leaf / exterior boxes). Parallel over leaves.

_maybe_kick!(p, pair, state, kick_after, spacings, n_fixed, n_protected, len_unit)
    -> (state, kicked)

Kick one point of the closest pair by 0.1·s in a random direction once the pair has stayed frozen (same indices and r/s) for kick_after consecutive iterations. Prefers an index past n_protected (a volume point) and always picks a movable one (past n_fixed). state is the (pair, rs, count) tuple the caller threads through.

_near_duplicate_keep_mask(pts, spacings, ratio) -> BitVector

Greedy, order-preserving keep-mask: drop any point closer than ratio·spacing to a kept, lower-indexed point (so boundary points, indexed first, survive over volume points). The ball search at the largest cull radius sees every point inside the threshold, so the guarantee holds for clusters of any size.

_non_exterior_leaves(node_tree, classification) -> Vector{Int}

Indices of every leaf that is not classified LEAF_EXTERIOR (interior + boundary). The shared leaf set behind the spacing-integral helpers below.

_project_to_boundary(sv, octree, offset_dist) -> (SVector, tri_idx)

Nearest point on the mesh surface, nudged offset_dist inward along the triangle normal. tri_idx == 0 means no triangle was found.

_rand_point_in_box(bbox_min, bbox_max)

Generate a random point uniformly distributed in a bounding box.

_raw_point(pt) -> SVector

Unitless coordinates of a point in its native machine type (dimension-generic).

_reconstruct_cloud(cloud, p, tri_indices, is_bnd, n_boundary, octree, spacing, keep)

Rebuild a PointCloud after octree repel: kept points are partitioned by is_bnd into a single :boundary surface and the volume. Projected boundary points take the landing triangle's normal; imported ones (id ≤ n_boundary) keep their original area, deposited ones get spacing².

_relax!(p, p_old, snap, spacing, force_model, constrain; kwargs...) -> Vector{T}

Shared relaxation loop behind both repel methods. p holds the movable points (updated in place); snap is the search snapshot whose first n_fixed entries are static and whose tail mirrors p, refreshed every rebuild_every iterations. constrain(id, xi, x_proposed) maps a proposed position to the final one (identity, or the octree wall rule). deposit!(p, method, i), when given, runs serially after each sweep. Kick targets prefer indices past n_protected. cv_target > 0 and stall_after > 0 add the quality-based stops: end when the movable points' dNN/s CV (read off the sweep's nearest-neighbor data, no extra search) reaches the target, or has not improved for that many iterations. Returns the force-norm convergence history `maxi(|Fi|·si)in the points' machine typeT`.

_safe_direction(xi, xj, r) -> Vec

Unit direction from xj to xi; a random unit vector when r == 0, avoiding the 0/0 NaN that traps coincident points.

Signed volume of a closed triangle mesh (divergence theorem, Σ dot(v1, v2 × v3) / 6). Positive iff the winding orients normals outward. Catches globally inside-out meshes, which has_consistent_normals cannot (a perfectly consistent but inverted mesh classifies its complement as interior — exactly the 2026-06-11 cavity corruption #2). Only meaningful for closed surfaces.

_triangle_axis_test(
    axis::SVector{3,T},
    v0::SVector{3,T},
    v1::SVector{3,T},
    v2::SVector{3,T},
    half::SVector{3,T}
) where {T<:Real} -> Bool

Internal helper for triangle-box intersection separating axis test.

Tests if the projection intervals of the triangle vertices and box overlap along the given axis.

Arguments

• axis: Separating axis direction
• v0, v1, v2: Triangle vertices in box-centered coordinates
• half: Box half-extents

Returns

true if intervals overlap (potential intersection), false if separated

add_box!(octree::SpatialOctree, i::Int, j::Int, k::Int, N::Int, parent_idx::Int) -> Int

Add new box to octree. Returns box index.

Automatically grows arrays if capacity exceeded.

Arguments

• i, j, k: Integer coordinates at level N
• N: Refinement level (N=1 is root, N=2 is first subdivision, etc.)
• parent_idx: Index of parent box

Returns

Index of newly created box

all_boxes(octree::SpatialOctree) -> Vector{Int}

Return indices of all boxes (leaves and internal nodes).

all_leaves(octree::SpatialOctree) -> Vector{Int}

Return indices of all leaf boxes.

any_leaf_overlapping(tree::SpatialOctree, bbox_min, bbox_max, predicate) -> Bool

Return true if any leaf whose bounding box overlaps [bbox_min, bbox_max] satisfies predicate(leaf_idx). Descends from root and prunes subtrees that cannot overlap, giving O(log L) expected cost instead of O(L) scan.

balance_octree!(octree::SpatialOctree, criterion::SubdivisionCriterion)

Enforce 2:1 balance constraint across entire octree.

Iteratively subdivides boxes that violate the 2:1 constraint until all adjacent boxes differ by at most one refinement level.

Arguments

• criterion: Subdivision criterion (only size constraints are enforced)

Algorithm

1. Collect all leaves
2. Check each leaf for balance violations
3. Subdivide violating neighbors (only respecting physical minimum-size limits)
4. Repeat until no violations

Note

• Uses can_subdivide (not should_subdivide) to ignore element count criteria
• Balancing is a geometric constraint, not an optimization decision
• Maximum iterations limit prevents infinite loops. If hit, tree may not be fully balanced.

bounding_box(octree::SpatialOctree) -> (SVector{3}, SVector{3})

Get overall bounding box of tree (root box).

box_bounds(octree::SpatialOctree, box_idx::Int) -> (SVector{3}, SVector{3})

Get (mincorner, maxcorner) of box.

Returns

Tuple of (mincorner, maxcorner) as SVector{3,C}

box_center(octree::SpatialOctree, box_idx::Int) -> SVector{3}

Compute spatial center of box using (i,j,k,N) coordinates.

Formula

center = origin + (2*[i,j,k] + 1) * box_size / 2

where box_size = root_size / N

box_size(tree::AbstractSpatialTree, box_idx::Int) -> Real

Get edge length of box.

Required

Trees must implement this function.

box_size(octree::SpatialOctree, box_idx::Int) -> Real

Get edge length of box.

Box size at refinement level N is: root_size / N

build_node_octree(triangle_octree, spacing, alpha, node_min_ratio)

Build a spacing-driven node octree from an existing triangle octree.

Creates a new SpatialOctree that subdivides based on a spacing function, enabling spacing-aware point distribution. The node octree is:

1. Recursively subdivided using SpacingCriterion
2. Balanced to maintain 2:1 refinement ratio
3. Independent of the triangle octree resolution

Arguments

• triangle_octree: Base TriangleOctree for geometry
• spacing: Spacing function (e.g., ConstantSpacing, BoundaryLayerSpacing)
• alpha: Subdivision aggressiveness (h_box ≤ alpha * h_spacing)
• node_min_ratio: Minimum box size ratio relative to domain

Returns

SpatialOctree{Int, T} with spacing-driven subdivision, where T is the triangle octree's coordinate type (the mesh CRS machine type)

Example

tri_octree = TriangleOctree(mesh; classify_leaves=true)
spacing = BoundaryLayerSpacing(points; at_wall=0.5m, bulk=5m, layer_thickness=2m)
node_tree = build_node_octree(tri_octree, spacing, 1.0, 1e-6)

can_subdivide(criterion::SubdivisionCriterion, tree, box_idx) -> Bool

Check if box CAN be subdivided based on physical constraints only.

Unlike should_subdivide, this ignores content-based criteria (like element count) and only checks physical limits (like minimum size). Used during balancing where subdivision is required for geometric correctness, not optimization.

Arguments

• criterion: Subdivision criterion (only size constraints are checked)
• tree: Spatial tree
• box_idx: Index of box to check

Returns

true if box can physically be subdivided, false if at minimum size limit.

Example

# For balancing, we only respect size limits
if needs_balancing(leaf)
    if can_subdivide(criterion, tree, leaf)
        subdivide!(tree, leaf)
    end
end

classify_leaves!(
    tree::SpatialOctree,
    geometry_query::Function;
    tolerance_relative=1e-6,
    tolerance_absolute=1e-8
) -> Vector{Int8}

Generic conservative classification of octree leaves using a geometry query function.

Arguments

• tree: The octree to classify
• geometry_query: Function (point::SVector{3,T}, tol::T) -> Int8 returning LEAF_INTERIOR, LEAF_BOUNDARY, or LEAF_EXTERIOR
• tolerance_relative: Relative tolerance scaled by box size
• tolerance_absolute: Absolute tolerance floor

Classification Strategy

Uses 9-point conservative probing (center + 8 actual box corners):

• LEAF_INTERIOR: All 9 points are interior
• LEAF_EXTERIOR: All 9 points are exterior
• LEAF_BOUNDARY: Mixed results or any boundary point

Note: Uses actual box corners (not inset) for accurate classification in high-curvature regions.

Returns

Vector of Int8 classifications for each box:

• -1 (LEAF_UNKNOWN): Non-leaf nodes
• 0 (LEAF_EXTERIOR): Exterior leaf
• 1 (LEAF_BOUNDARY): Boundary leaf
• 2 (LEAF_INTERIOR): Interior leaf

classify_node_octree(node_tree, triangle_octree)

Classify node octree leaves as interior, boundary, or exterior using the triangle octree's geometry query. Correctness of downstream sampling (skipping isinside on LEAF_INTERIOR points) relies on the mesh-bbox early return inside _mesh_geometry_query, which prevents sign-vote flips from promoting far-exterior leaves into LEAF_INTERIOR.

Returns

Vector of Int8 classifications indexed by node-octree box index.

closest_point_on_triangle(
    P::SVector{3,T},
    v1::SVector{3,T},
    v2::SVector{3,T},
    v3::SVector{3,T}
) where {T<:Real} -> SVector{3,T}

Compute the closest point on triangle (v1, v2, v3) to point P.

Uses barycentric coordinate method from Ericson's "Real-Time Collision Detection". The closest point is computed by:

1. Projecting P onto the triangle plane
2. Computing barycentric coordinates
3. Clamping to triangle if outside

Algorithm

The triangle can be parameterized as: T(u,v) = v1 + u(v2-v1) + v(v3-v1) for u,v ≥ 0, u+v ≤ 1

We find the closest point by solving a constrained minimization problem.

Returns

Point on triangle surface closest to P (may be on edge or vertex).

References

Ericson, "Real-Time Collision Detection", Chapter 5.1.5

closest_point_on_triangle_feature(P, v1, v2, v3) -> (SVector{3,T}, Int8)

Same as closest_point_on_triangle but also reports which triangle feature the closest point lies on (FEATURE_FACE, FEATURE_VERTEX_k, FEATURE_EDGE_jk). The feature is a free by-product of the Ericson region classification and selects the correct angle-weighted pseudonormal for exact signed-distance computation (Bærentzen & Aanæs 2005).

distance_point_triangle(
    P::SVector{3,T},
    v1::SVector{3,T},
    v2::SVector{3,T},
    v3::SVector{3,T}
) where {T<:Real} -> T

Compute unsigned distance from point P to triangle (v1, v2, v3).

Returns

Unsigned distance (always positive or zero)

distance_point_triangle(
    P::SVector{3,T},
    v1::SVector{3,T},
    v2::SVector{3,T},
    v3::SVector{3,T},
    normal::SVector{3,T}
) where {T<:Real} -> T

Compute signed distance from point P to triangle (v1, v2, v3).

The distance is:

• Positive if P is on the side of the triangle that the normal points to
• Negative if P is on the opposite side
• Zero if P is on the triangle plane

Algorithm

1. Find closest point Q on triangle to P
2. Compute distance ||P - Q||
3. Determine sign based on which side of triangle P is on

Arguments

• P: Query point
• v1, v2, v3: Triangle vertices in counterclockwise order
• normal: Outward-pointing unit normal vector

Returns

Signed distance (positive = outside, negative = inside for closed surface)

Example

using StaticArrays

# Triangle in xy-plane
v1 = SVector(0.0, 0.0, 0.0)
v2 = SVector(1.0, 0.0, 0.0)
v3 = SVector(0.0, 1.0, 0.0)
normal = SVector(0.0, 0.0, 1.0)

# Point above triangle
P = SVector(0.25, 0.25, 1.0)
d = distance_point_triangle(P, v1, v2, v3, normal)
# d ≈ 1.0 (positive, on normal side)

emptyspace(testpoint, points)

Check if a point occupies empty space within a certain tolerance.

find_boxes_at_coords(octree::SpatialOctree, i_target::Int, j_target::Int, k_target::Int, N_target::Int) -> Vector{Int}

Find box(es) at given (i,j,k,N) coordinates.

• If exact match found at level Ntarget, returns [boxidx]
• If location is covered by coarser box, returns [coarserboxidx]
• If location is subdivided finer, returns all descendants at that location

Returns

Vector of box indices covering the target location

find_leaf(octree::SpatialOctree, point::SVector{3}) -> Int

Find leaf box containing point. Returns box index.

Traverses tree from root to leaf in O(log n) time.

Arguments

• point: Query point in same coordinate system as octree

Returns

Index of leaf box containing point

Throws

AssertionError if point is outside octree bounds

find_neighbor(octree::SpatialOctree, box_idx::Int, direction::Int) -> Vector{Int}

Find neighbor(s) in given direction. Returns vector of neighbor indices.

Handles 2:1 refinement level difference:

• If neighbor exists at same level: returns [neighbor_idx]
• If neighbor is subdivided (finer): returns children on shared face (up to 4)
• If neighbor doesn't exist (boundary): returns empty vector

Arguments

• box_idx: Box to find neighbor of
• direction: Direction code (1-6, see neighbor_direction)

Returns

Vector of neighbor box indices (may be empty if at boundary)

global_to_local(boundary::PointBoundary, global_idx::Int) -> (Symbol, Int)

Convert a boundary-global index to a (surface_name, local_index) tuple.

global_to_local(cloud::PointCloud, global_idx::Int) -> (Symbol, Int)

Convert a cloud-global index to a (component, local_index) tuple. Returns (:volume, local_idx) for volume indices, or (surface_name, local_idx) for boundary indices.

has_children(tree::AbstractSpatialTree, box_idx::Int) -> Bool

Check if box has been subdivided.

Required

Trees must implement this function.

has_children(octree::SpatialOctree, box_idx::Int) -> Bool

Check if box has been subdivided.

is_leaf(tree::AbstractSpatialTree, box_idx::Int) -> Bool

Check if box is a leaf (has no children).

Required

Trees must implement this function.

is_leaf(octree::SpatialOctree, box_idx::Int) -> Bool

Check if box has no children.

local_to_global(boundary::PointBoundary, name::Symbol, local_idx::Int) -> Int

Convert a surface-local index to a boundary-global index.

local_to_global(cloud::PointCloud, name::Symbol, local_idx::Int) -> Int

Convert a surface-local index to a cloud-global index.

needs_balancing(octree::SpatialOctree, box_idx::Int) -> Bool

Check if subdividing this box would violate 2:1 balance with any neighbor.

Returns true if any neighbor has grandchildren (2-level refinement difference).

neighbor_direction(direction::Int) -> (Int, Int, Int)

Convert direction code to (di, dj, dk) offset.

Direction Codes

• 1: -x (left)
• 2: +x (right)
• 3: -y (bottom)
• 4: +y (top)
• 5: -z (front)
• 6: +z (back)

num_elements(octree::SpatialOctree) -> Int

Total number of elements stored in tree.

should_subdivide(criterion::SubdivisionCriterion, tree, box_idx) -> Bool

Determine if box should be subdivided based on criterion.

Arguments

• criterion: Subdivision criterion to evaluate
• tree: Spatial tree
• box_idx: Index of box to check

Returns

true if box should be subdivided, false otherwise.

subdivide!(octree::SpatialOctree, box_idx::Int) -> SVector{8,Int}

Subdivide box into 8 children. Returns child indices [1:8].

Child Ordering (Standard Octree Convention)

1: (0,0,0) - bottom-left-front (x-, y-, z-) 2: (1,0,0) - bottom-right-front (x+, y-, z-) 3: (0,1,0) - top-left-front (x-, y+, z-) 4: (1,1,0) - top-right-front (x+, y+, z-) 5: (0,0,1) - bottom-left-back (x-, y-, z+) 6: (1,0,1) - bottom-right-back (x+, y-, z+) 7: (0,1,1) - top-left-back (x-, y+, z+) 8: (1,1,1) - top-right-back (x+, y+, z+)

Arguments

• box_idx: Index of box to subdivide (must be leaf)

Returns

SVector{8,Int} of child indices in standard order

surface_offset(boundary::PointBoundary, name::Symbol) -> Int

Return the global offset for the start of surface name within the boundary. The first point of surface name has boundary-global index surface_offset(...) + 1.

triangle_box_intersection(
    v1::SVector{3,T},
    v2::SVector{3,T},
    v3::SVector{3,T},
    box_min::SVector{3,T},
    box_max::SVector{3,T}
) where {T<:Real} -> Bool

Test if triangle (v1, v2, v3) intersects axis-aligned box.

Uses the Separating Axis Theorem (SAT) with 13 potential separating axes:

• 3 box face normals (x, y, z axes)
• 1 triangle normal
• 9 edge-edge cross products

If any axis separates the triangle and box, they don't intersect.

Algorithm

1. Translate triangle and box so box is centered at origin
2. Test each potential separating axis
3. Return false if any axis separates, true otherwise

References

• Akenine-Möller, "Fast 3D Triangle-Box Overlap Testing" (2001)
• Ericson, "Real-Time Collision Detection", Chapter 5.2.9

Performance

Optimized with early-out tests. Average case is much faster than worst case.

volume_to_global(cloud::PointCloud, local_idx::Int) -> Int

Convert a volume-local index to a cloud-global index.