JuliaMeshless · kylebeggs · Feb 18, 2026 · Feb 18, 2026 · Feb 18, 2026 · Feb 18, 2026
diff --git a/README.md b/README.md
@@ -19,18 +19,16 @@
 
 ---
 
-*Also known as: meshfree methods, meshless methods, RBF-FD (Radial Basis Function Finite Differences), scattered data interpolation, kernel methods, polyharmonic splines, collocation methods*
+## Why RadialBasisFunctions.jl?
 
----
+RBF methods approximate functions from scattered data without requiring a mesh. This package focuses on **local collocation** — building stencils from k-nearest neighbors instead of coupling all N points — so it scales to large problems without dense N×N solves.
 
-## Why RadialBasisFunctions.jl?
+Beyond interpolation, it provides differential operators (Laplacian, gradient, partials, custom) that act directly on scattered data, with support for Hermite interpolation to enforce boundary conditions in PDE applications.
 
-- **Truly Meshless** - Work directly with scattered points, no mesh generation required
-- **GPU Ready** - Seamless CPU/GPU execution via [KernelAbstractions.jl](https://github.com/JuliaGPU/KernelAbstractions.jl)
-- **Flexible Operators** - Built-in Laplacian, gradient, partial derivatives, or define your own custom operators
-- **Boundary Conditions** - Hermite interpolation for Dirichlet, Neumann, and Robin BCs in PDE applications
-- **Local Collocation** - k-nearest neighbor stencils for scalable large-domain problems
-- **Differentiable** - AD-compatible via Enzyme and Mooncake extensions
+Other things that might matter to you:
+- GPU-accelerated weight computation via [KernelAbstractions.jl](https://github.com/JuliaGPU/KernelAbstractions.jl)
+- Native autodiff rules for [Enzyme.jl](https://github.com/EnzymeAD/Enzyme.jl) and [Mooncake.jl](https://github.com/compintell/Mooncake.jl) (no generic fallbacks)
+- Operator algebra — combine operators with `+`, `-`, `*`
 
 <p align="center">
   <img src="docs/src/assets/interpolation_demo.png" alt="RBF Interpolation Demo" width="800">
@@ -43,42 +41,41 @@
 ```julia
 using RadialBasisFunctions, StaticArrays
 
-# Create scattered 2D points
+# Scattered 2D points — no mesh needed
 points = [SVector{2}(rand(2)) for _ in 1:500]
 f(x) = sin(4x[1]) * cos(3x[2])
 values = f.(points)
 
-# Interpolation - evaluate anywhere in the domain
+# Interpolation
 interp = Interpolator(points, values)
 interp(SVector(0.5, 0.5))
 
-# Differential operators - apply to scattered data
-lap = laplacian(points)       # Laplacian operator
-grad = gradient(points)       # Gradient operator
-∂x = partial(points, 1, 1)    # ∂/∂x₁
+# Differential operators on scattered data
+∇²  = laplacian(points)
+∇   = gradient(points)
+∂x  = partial(points, 1, 1)       # ∂/∂x₁
+∂²y = partial(points, 2, 2)       # ∂²/∂x₂²
 
-lap_values = lap(values)      # Apply Laplacian to data
-grad_values = grad(values)    # Returns Nx2 matrix of gradients
+∇²(values)                         # apply to data
+∇(values)                          # Nx2 matrix
+
+# Combine operators
+mixed = ∂x + ∂²y                   # operator algebra
+
+# Transfer data between point sets
+target = [SVector{2}(rand(2)) for _ in 1:1000]
+rg = regrid(points, target)
+rg(values)                         # interpolated onto target
 ```
 
 ## Supported Radial Basis Functions
 
 | Type | Formula | Best For |
 |------|---------|----------|
-| Polyharmonic Spline (PHS) | $r^n$ where $n = 1, 3, 5, 7$ | General purpose, no shape parameter tuning needed |
+| Polyharmonic Spline (PHS) | $r^n$ where $n = 1, 3, 5, 7$ | General purpose, no shape parameter tuning |
 | Inverse Multiquadric (IMQ) | $1 / \sqrt{(\varepsilon r)^2 + 1}$ | Smooth interpolation with tunable accuracy |
 | Gaussian | $e^{-(\varepsilon r)^2}$ | Infinitely smooth functions |
 
-All basis functions support polynomial augmentation (enabled by default) for enhanced accuracy.
-
-```julia
-# Basis function examples
-basis = PHS(3)                  # Cubic PHS (default: quadratic polynomial augmentation)
-basis = PHS(5; poly_deg=3)      # Quintic PHS with cubic polynomials
-basis = IMQ(1.0)                # IMQ with shape parameter ε=1.0
-basis = Gaussian(0.5)           # Gaussian with shape parameter ε=0.5
-```
-
 ## Installation
 
 ```julia
@@ -88,42 +85,15 @@ Pkg.add("RadialBasisFunctions")
 
 Requires Julia 1.10 or later.
 
-## Advanced Features
-
-### Hermite Interpolation for PDEs
-
-Handle boundary conditions properly with Hermite interpolation:
-
-```julia
-# Build operator with boundary awareness
-lap = laplacian(points, eval_points, PHS(3), k,
-                boundary_nodes, boundary_conditions, normals)
-```
-
-### GPU Acceleration
-
-Operators automatically leverage GPU when data is on device:
-
-```julia
-using CUDA
-points_gpu = cu(points)
-lap_gpu = laplacian(points_gpu)  # Weights computed on GPU
-```
-
-### Custom Operators
-
-Define any differential operator using automatic differentiation:
-
-```julia
-# Example: custom operator
-my_op = custom(points, basis -> (x, xc) -> your_derivative_function(basis, x, xc))
-```
+> [!TIP]
+> The examples above use sensible defaults (`PHS(3)` basis, quadratic polynomial augmentation, auto-selected stencil size). All of these are configurable — see the **[Getting Started guide](https://JuliaMeshless.github.io/RadialBasisFunctions.jl/stable/getting_started)** for details on choosing a basis function, stencil size (`k`), polynomial degree, Hermite boundary conditions, and more.
 
 ## Documentation
 
-- **[Getting Started](https://JuliaMeshless.github.io/RadialBasisFunctions.jl/stable/getting_started)** - Tutorials and examples
-- **[Theory](https://JuliaMeshless.github.io/RadialBasisFunctions.jl/stable/theory)** - Mathematical background on RBF methods
-- **[API Reference](https://JuliaMeshless.github.io/RadialBasisFunctions.jl/stable/api)** - Full function documentation
+- **[Getting Started](https://JuliaMeshless.github.io/RadialBasisFunctions.jl/stable/getting_started)** — tutorials covering interpolation, operators, boundary conditions, and GPU usage
+- **[Autodiff](https://JuliaMeshless.github.io/RadialBasisFunctions.jl/stable/autodiff)** — differentiating through operators with Enzyme and Mooncake
+- **[Theory](https://JuliaMeshless.github.io/RadialBasisFunctions.jl/stable/theory)** — mathematical background on RBF methods
+- **[API Reference](https://JuliaMeshless.github.io/RadialBasisFunctions.jl/stable/api)** — full function documentation
 
 ## Citation
 

diff --git a/docs/src/autodiff.md b/docs/src/autodiff.md
@@ -1,6 +1,6 @@
 # Automatic Differentiation
 
-RadialBasisFunctions.jl supports automatic differentiation (AD) through two package extensions:
+Both operators and interpolators can be differentiated with reverse-mode AD. Two backends are supported through package extensions:
 
 - **Mooncake.jl** - Reverse-mode AD with support for mutation
 - **Enzyme.jl** - Native EnzymeRules for high-performance reverse-mode AD
@@ -15,7 +15,7 @@ See: [Enzyme.jl#2699](https://github.com/EnzymeAD/Enzyme.jl/issues/2699)
 
 ## Implementation Status
 
-Enzyme.jl has native `EnzymeRules` for all supported operations. Mooncake.jl support currently uses the `@from_rrule` macro to convert ChainRulesCore rules. Native Mooncake rules are planned for a future release.
+Both backends have native AD rule implementations. Enzyme.jl uses `EnzymeRules` (`augmented_primal`/`reverse`) and Mooncake.jl uses native `rrule!!` with `@is_primitive`.
 
 ## Differentiating Through Operators
 
@@ -78,6 +78,10 @@ grad[1:5]
 
 When differentiating through interpolation, the `Interpolator` must be constructed inside the loss function since changing the input values changes the interpolation weights.
 
+!!! note
+    Differentiating through `Interpolator` construction currently requires Mooncake.
+    Enzyme does not yet support this code path via DifferentiationInterface.
+
 ```@example autodiff
 N_interp = 30
 points_interp = [SVector{2}(0.5 + 0.4 * cos(2π * i / N_interp), 0.5 + 0.4 * sin(2π * i / N_interp)) for i in 1:N_interp]
@@ -179,11 +183,14 @@ grad[1:6]
 | Component | Enzyme | Mooncake |
 |-----------|:------:|:--------:|
 | Operator evaluation (`op(values)`) | ✓ | ✓ |
+| Interpolator construction | ✗* | ✓ |
 | Interpolator evaluation | ✓ | ✓ |
 | Basis functions (PHS, IMQ, Gaussian) | ✓ | ✓ |
 | Weight construction (`_build_weights`) | ✓ | ✓ |
 | Shape parameter (ε) differentiation | ✓ | ✓ |
 
+*Enzyme does not support differentiating through `Interpolator` construction via DifferentiationInterface due to `factorize` limitations. Use Mooncake for this use case.
+
 ## Using Enzyme Backend
 
 Switch to Enzyme by changing the backend (requires Julia < 1.12):

diff --git a/docs/src/getting_started.md b/docs/src/getting_started.md
@@ -1,6 +1,6 @@
 # Getting Started
 
-First, let's load the package along with the [StaticArrays.jl](https://github.com/JuliaArrays/StaticArrays.jl) package which we use for each data point
+Data must be an `AbstractVector` of point vectors — each point needs an inferrable dimension (e.g., `SVector{2,Float64}` from [StaticArrays.jl](https://github.com/JuliaArrays/StaticArrays.jl)). This includes `CuVector` for GPU.
 
 ```@example overview
 using RadialBasisFunctions
@@ -26,7 +26,7 @@ interp = Interpolator(x, y)
 and evaluate it at a new point
 
 ```@example overview
-x_new = [rand(2) for _ in 1:5]
+x_new = [SVector{2}(rand(2)) for _ in 1:5]
 y_new = interp(x_new)
 y_true = f.(x_new)
 ```
@@ -37,31 +37,20 @@ and compare the error
 abs.(y_true .- y_new)
 ```
 
-Wow! The error is numerically zero! Well... we set ourselves up for success here. `Interpolator` (along with `RadialBasisOperator`) has an optional argument to provide the type of radial basis including the degree of polynomial augmentation. The default basis is a cubic polyharmonic spline with 2nd degree polynomial augmentation (which the constructor is `PHS(3, poly_deg=2)`) and given the underlying function we are interpolating is a 2nd order polynomial itself, we are able to represent it exactly (up to machine precision). Let's see what happens when we only use 1st order polynomial augmentation
+The error is numerically zero because the default basis — `PHS(3; poly_deg=2)` — includes quadratic polynomial augmentation, which can represent our 2nd-order polynomial `f` exactly. Reducing the polynomial degree shows the effect:
 
 ```@example overview
-interp = Interpolator(x, y, PHS(3, poly_deg=1))
+interp = Interpolator(x, y, PHS(3; poly_deg=1))
 y_new = interp(x_new)
 abs.(y_true .- y_new)
 ```
 
 ## Operators
 
-This package also provides an API for operators. There is support for several built-in operators along with support for user-defined operators:
-
-- partial derivative (1st and 2nd order)
-- laplacian
-- gradient / jacobian
-- directional derivative
-- custom (user-defined) operators
-- regridding (interpolation between point sets)
-
-Please make an issue or pull request for additional operators.
+Operators compute RBF-FD weights for differentiation on scattered data. They're lazy — weights are built on first evaluation and cached.
 
 ### Partial Derivative
 
-We can take the same data as above and build a partial derivative operator with a similar construction as the interpolator. For the partial we need to specify the order of differentiation we want along with the dimension for which to take the partial. We can also supply some optional arguments such as the basis and number of points in the stencil. The function inputs order is `partial(x, order, dim, basis; k=5)`
-
 ```@example overview
 df_x_rbf = partial(x, 1, 1)
 
@@ -74,8 +63,6 @@ all(abs.(df_x.(x) .- df_x_rbf(y)) .< 1e-10)
 
 ### Laplacian
 
-Building a laplacian operator is as easy as
-
 ```@example overview
 lap_rbf = laplacian(x)
 
@@ -219,6 +206,6 @@ typeof(result)
 
 ## Current Limitations
 
-1. **Data format**: The package requires `Vector{AbstractVector}` input (not matrices). Each point must have inferrable dimension, e.g., `SVector{2,Float64}` from StaticArrays.jl. Matrix input support is planned.
+1. **Data format**: The package requires `AbstractVector{<:AbstractVector}` input (not matrices). Each point must have inferrable dimension, e.g., `SVector{2,Float64}` from StaticArrays.jl. Matrix input support is planned.
 
 2. **Global interpolation**: `Interpolator` currently uses all points globally. Local collocation support (like the operators use) is planned for future releases.
diff --git a/docs/src/index.md b/docs/src/index.md
@@ -4,8 +4,8 @@ layout: home
 
 hero:
   name: RadialBasisFunctions.jl
-  text: Meshless Methods Made Easy
-  tagline: RBF interpolation, differential operators, and PDE tools for Julia
+  text: Meshless Computing in Julia
+  tagline: Radial basis functions for operators, machine learning, and beyond.
   actions:
     - theme: brand
       text: Get Started
@@ -15,49 +15,60 @@ hero:
       link: https://github.com/JuliaMeshless/RadialBasisFunctions.jl
 
 features:
-  - icon: 📐
-    title: Interpolation
-    details: Scattered data interpolation with polynomial augmentation
+  - icon: ☁️
+    title: No Mesh Required
+    details: Interpolate and differentiate directly on scattered point clouds. Local stencils from k-nearest neighbors scale to large problems.
   - icon: ∇
-    title: Differential Operators
-    details: Convenience constructors for common operators such as partial derivatives, laplacian, and gradient
+    title: Operators
+    details: "Laplacian, gradient, partials, directional derivatives, and custom operators — with operator algebra to combine them."
   - icon: 🚀
     title: GPU Ready
-    details: KernelAbstractions.jl enables seamless CPU and GPU execution
-  - icon: 🔧
-    title: Extensible
-    details: Custom operators and Hermite boundary conditions for PDE applications
+    details: Weight computation parallelizes over stencils via KernelAbstractions.jl. Same code runs on CPU and GPU.
+  - icon: ∂
+    title: Fully Differentiable
+    details: Native AD rules for Enzyme.jl and Mooncake.jl — differentiate through operators, interpolators, and weight construction.
 ---
 ```
 
-## Quick Example
+## Quick Start
 
 ```julia
 using RadialBasisFunctions, StaticArrays
 
-# Create scattered data points
-x = [SVector{2}(rand(2)) for _ in 1:1000]
-y = sin.(norm.(x))
+# Scattered data
+points = [SVector{2}(rand(2)) for _ in 1:500]
+f(x) = sin(4x[1]) * cos(3x[2])
+values = f.(points)
 
-# Build an interpolator
-interp = Interpolator(x, y)
-interp(SVector{2}(rand(2)))  # Evaluate at points
+# Interpolation
+interp = Interpolator(points, values)
+interp(SVector(0.5, 0.5))
 
-# Build differential operators
-lap = laplacian(x)
-grad = gradient(x)
-∂x = partial(x, 1, 1)  # ∂/∂x₁
+# Differential operators on scattered data
+∇²  = laplacian(points)
+∇   = gradient(points)
+∂x  = partial(points, 1, 1)       # ∂/∂x₁
+∂²y = partial(points, 2, 2)       # ∂²/∂x₂²
+
+∇²(values)                         # apply to data
+∇(values)                          # Nx2 matrix
+
+# Combine operators
+mixed = ∂x + ∂²y                   # operator algebra
+
+# Transfer data between point sets
+target = [SVector{2}(rand(2)) for _ in 1:1000]
+rg = regrid(points, target)
+rg(values)                         # interpolated onto target
 ```
 
 ## Supported Radial Basis Functions
 
-| Type                 | Function                              |
-| -------------------- | ------------------------------------- |
-| Polyharmonic Spline  | $r^n$ where $n = 1, 3, 5, 7$          |
-| Inverse Multiquadric | $1 / \sqrt{(\varepsilon r)^2 + 1}$    |
-| Gaussian             | $e^{-(\varepsilon r)^2}$              |
-
-All basis functions support optional polynomial augmentation for enhanced accuracy near boundaries.
+| Type | Formula | Best For |
+|------|---------|----------|
+| Polyharmonic Spline (PHS) | $r^n$ where $n = 1, 3, 5, 7$ | General purpose, no shape parameter tuning |
+| Inverse Multiquadric (IMQ) | $1 / \sqrt{(\varepsilon r)^2 + 1}$ | Smooth interpolation with tunable accuracy |
+| Gaussian | $e^{-(\varepsilon r)^2}$ | Infinitely smooth functions |
 
 ## Installation
 
@@ -66,8 +77,4 @@ using Pkg
 Pkg.add("RadialBasisFunctions")
 ```
 
-Or in the REPL:
-
-```julia
-] add RadialBasisFunctions
-```
+Requires Julia 1.10 or later.
diff --git a/docs/src/quickref.md b/docs/src/quickref.md
@@ -1,10 +1,8 @@
 # Quick Reference
 
-A compact cheat sheet for RadialBasisFunctions.jl.
-
 ## Data Format
 
-Data **must** be `Vector{AbstractVector}`, not a `Matrix`:
+Data **must** be an `AbstractVector` of point vectors, not a `Matrix`:
 
 ```julia
 using StaticArrays
@@ -25,8 +23,8 @@ points_3d = [SVector{3}(rand(3)) for _ in 1:100]
 | Type | Constructor | Formula | Shape Parameter |
 |------|-------------|---------|-----------------|
 | Polyharmonic Spline | `PHS(n)` | $r^n$ | None (scale-free) |
-| Inverse Multiquadric | `IMQ(ε)` | $\frac{1}{\sqrt{(\varepsilon r)^2 + 1}}$ | Required |
-| Gaussian | `Gaussian(ε)` | $e^{-(\varepsilon r)^2}$ | Required |
+| Inverse Multiquadric | `IMQ(ε)` | $\frac{1}{\sqrt{(\varepsilon r)^2 + 1}}$ | Optional (default: 1) |
+| Gaussian | `Gaussian(ε)` | $e^{-(\varepsilon r)^2}$ | Optional (default: 1) |
 
 ### PHS Orders
 
@@ -117,7 +115,7 @@ For PDE problems with boundary conditions:
 ```julia
 # Prepare boundary data
 is_boundary = [is_on_boundary(p) for p in points]
-bc = [Dirichlet() for _ in points]  # or Neumann(), Robin(α)
+bc = [Dirichlet() for _ in points]  # or Neumann(), Robin(α, β)
 normals = [compute_normal(p) for p in points]
 
 # Create operator with Hermite interpolation
@@ -133,7 +131,7 @@ lap = laplacian(
 |------|-------------|---------|
 | Dirichlet | `Dirichlet()` | Fixed value |
 | Neumann | `Neumann()` | Fixed normal derivative |
-| Robin | `Robin(α)` | $\alpha u + \frac{\partial u}{\partial n}$ |
+| Robin | `Robin(α, β)` | $\alpha u + \beta \frac{\partial u}{\partial n}$ |
 
 ## GPU Acceleration