Quick Reference

Data Format

Data must be an AbstractVector of point vectors, not a Matrix:

using StaticArrays

# CORRECT: Vector of static vectors
points = rand(SVector{2,Float64}, 100)

# Converting from matrix
matrix = rand(100, 2)
points = map(SVector{2}, eachrow(matrix))

# 3D points
points_3d = rand(SVector{3,Float64}, 100)

Basis Functions

Type	Constructor	Formula	Shape Parameter
Polyharmonic Spline	PHS(n)		None (scale-free)
Inverse Multiquadric	IMQ(ε)		Optional (default: 1)
Gaussian	Gaussian(ε)		Optional (default: 1)

PHS Orders

Order	Constructor	Smoothness	Use Case
1	PHS(1)	C⁰	Rough data
3	PHS(3)	C²	General purpose (default)
5	PHS(5)	C⁴	Smooth functions
7	PHS(7)	C⁶	Very smooth functions

Polynomial Augmentation

# Default: quadratic (poly_deg=2)
basis = PHS(3)

# Custom polynomial degree
basis = PHS(3; poly_deg=0)   # Constant only
basis = PHS(3; poly_deg=1)   # Linear
basis = PHS(3; poly_deg=3)   # Cubic
basis = PHS(3; poly_deg=-1)  # No polynomial (not recommended)

Operators

Creating Operators

using RadialBasisFunctions, StaticArrays

points = rand(SVector{2,Float64}, 100)

# Scalar field operators
lap = laplacian(points)           # ∇²f (scalar output)
grad = gradient(points)           # ∇f (vector output)
∂x = partial(points, 1, 1)        # ∂f/∂x (1st order, dim 1)
∂²y = partial(points, 2, 2)       # ∂²f/∂y² (2nd order, dim 2)
∂xy = mixed_partial(points, 1, 2) # ∂²f/∂x∂y
H = hessian(points)               # Full Hessian matrix
∂v = directional(points, v)       # ∇f·v
∂n = normal_derivative(points, normals)  # ∇f·n̂
diff = (@operator ∇ ⋅ (κ * ∇))(points)          # ∇⋅(κ∇f)

# Vector field operators (input: N×D matrix)
div_op = divergence(points)       # ∇⋅u (scalar output)
curl_op = curl(points)            # ∇×u (2D: scalar, 3D: vector)
S = strain_rate(points)           # ½(∇u + (∇u)ᵀ)
R = rotation_rate(points)         # ½(∇u − (∇u)ᵀ)

# Interpolation operators
rg = regrid(source, target)       # Local interpolation

Applying Operators

values = sin.(getindex.(points, 1))

# Apply to data
lap_values = lap(values)          # Vector of scalars
grad_values = grad(values)        # Matrix (N × dim)

Common Options

# Custom basis
lap = laplacian(points; basis=PHS(5; poly_deg=3))

# Custom stencil size
lap = laplacian(points; k=30)

# Different evaluation points
lap = laplacian(points; eval_points=other_points)

# Precomputed neighbors
adjl = find_neighbors(points, k)
lap = laplacian(points; adjl=adjl)

Global Interpolation

# Create interpolator (uses all points)
interp = Interpolator(points, values)

# Evaluate at single point
result = interp(SVector(0.5, 0.5))

# Evaluate at multiple points
results = interp(new_points)

Hermite Boundary Conditions

Type	Constructor	Meaning
Dirichlet	Dirichlet()	Fixed value
Neumann	Neumann()	Fixed normal derivative
Robin	Robin(α, β)
Internal	Internal()	Interior point, no condition

using RadialBasisFunctions: Dirichlet, Neumann, Robin  # not exported
lap = laplacian(points; hermite=(is_boundary=…, bc=…, normals=…))

See Boundary Conditions for the full walkthrough.

GPU Evaluation

Weights are always built on CPU. Once built, cu(op) (via Adapt) moves the operator's StencilWeights — both the dense weight-value matrix and the Int32 neighbor-index matrix — to the device, and applying the operator runs a KernelAbstractions gather kernel directly on the GPU:

using CUDA, Adapt

# Build operator on CPU
lap = laplacian(points)

# Move to GPU — both the weight values and the stencil indices transfer
lap_gpu = cu(lap)
values_gpu = cu(values)
result_gpu = lap_gpu(values_gpu)  # KernelAbstractions kernel on device

The adjoint apply (weights(op)' * x) runs on device through the same mechanism — a transpose-map gather kernel. Known limitations: update_weights! on a device-adapted operator throws (weights are always built on CPU — rebuild there and re-adapt), and differentiating through a device-adapted operator is untested; keep AD workloads on CPU.

GPU weight building is planned — see #88.

Common Errors

Error messages, causes, and fixes are covered in Troubleshooting.

Performance Tips

2. Reuse adjacency lists for multiple operators on same points

3. Use StaticArrays (SVector) for best performance

4. Batch operations - create operator once, apply many times

5. GPU for large problems — build operators on CPU, then cu(operator); evaluation runs a device-native KernelAbstractions kernel