Boundary Conditions

For PDE applications, operators support Hermite interpolation with boundary conditions. Use this when you need to enforce Dirichlet, Neumann, or Robin conditions at boundary nodes while solving the PDE only at interior nodes.

For the mathematics — why standard collocation breaks down near boundaries and how the Hermite formulation restores symmetry — see Hermite Approach for Boundary Stencils in the theory reference.

Boundary Condition Types

Type	Constructor	Meaning
Dirichlet	Dirichlet()	Value specified:
Neumann	Neumann()	Normal derivative specified:
Robin	Robin(α, β)	Mixed:
Internal	Internal()	Interior point (no boundary condition)

Dirichlet, Neumann, and Robin are public but not exported — downstream physics packages (e.g. Macchiato.jl) export their own boundary-condition types with these names. Import them explicitly:

using RadialBasisFunctions: Dirichlet, Neumann, Robin

Setting Up a Domain

Hermite stencils need a real boundary with meaningful outward normals. Below is a unit disk: interior nodes scattered inside, boundary nodes placed on the circle itself.

using RadialBasisFunctions
using RadialBasisFunctions: Dirichlet, Neumann, Robin
using StaticArrays
using LinearAlgebra: norm, normalize
using Random

Random.seed!(42)

# Boundary nodes on the unit circle
n_boundary = 120
θ = range(0, 2π; length=n_boundary + 1)[1:(end - 1)]
boundary_pts = map(t -> SVector(cos(t), sin(t)), θ)

# Interior nodes scattered inside the disk
interior_pts = SVector{2,Float64}[]
while length(interior_pts) < 700
    p = 2 * rand(SVector{2,Float64}) .- 1
    norm(p) < 0.95 && push!(interior_pts, p)
end

points = vcat(interior_pts, boundary_pts)
is_boundary = [i > length(interior_pts) for i in eachindex(points)]
boundary_indices = findall(is_boundary)

# Outward unit normals — radial on a disk
normals = normalize.(points[boundary_indices])
length(points), count(is_boundary)

(820, 120)

The hermite Keyword

Pass a named tuple with three fields to any operator constructor. hermite is keyword-only (since v0.6):

• is_boundary — a Bool per data point, marking which points lie on the boundary

• bc — a boundary condition per boundary point, ordered as findall(is_boundary)

• normals — an outward unit normal per boundary point, same order

bcs = fill(Dirichlet(), length(boundary_indices))

lap = laplacian(points; hermite=(
    is_boundary=is_boundary,
    bc=bcs,
    normals=normals
))
typeof(lap)

RadialBasisOperator{Laplacian, SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{StaticArraysCore.SVector{2, Float64}}, Vector{StaticArraysCore.SVector{2, Float64}}, Vector{Vector{Int64}}, PHS3{Int64}, KernelAbstractions.CPU}

We can check it against a harmonic function,   , for which  :

u(p) = p[1]^2 - p[2]^2

result = lap(u.(points))
maximum(abs, result[.!is_boundary])   # ≈ 0 at interior nodes

2.9103830456733704e-11

What Goes Into the Input Vector

This is the part that trips people up. A Hermite operator assembles as

so entries of the vector you pass mean different things depending on the node:

Node	What that entry must hold
Interior	the field value
Dirichlet boundary	— same thing, so passing u everywhere just works
Neumann boundary	— not
Robin boundary

Dirichlet is forgiving because its boundary data is the field value. Neumann is not. On the unit circle,     :

∂u∂n(p) = 2 * p[1]^2 - 2 * p[2]^2

lap_neumann = laplacian(points; hermite=(
    is_boundary=is_boundary,
    bc=fill(Neumann(), length(boundary_indices)),
    normals=normals
))

v = map((b, p) -> b ? ∂u∂n(p) : u(p), is_boundary, points)
maximum(abs, lap_neumann(v)[.!is_boundary])

2.9103830456733704e-11

Passing u at the Neumann nodes instead gives an error many orders of magnitude larger — the operator is doing exactly what it was built to do, just with the wrong data.

Condition types can be mixed freely; build the bc vector however the geometry requires, and populate the input vector per the table above.

Where It Applies

The Hermite treatment is applied only to stencils that include boundary nodes. Interior stencils far from the boundary use the standard RBF-FD formulation unchanged, so there is no cost for problems where boundary effects don't reach.

The same hermite keyword works on operators built with @operator — see Building PDE Operators — and on Custom Operators.

Alternative: Boundary Nodes as Unknowns

For multi-region or coupled problems it can be preferable to solve the governing equation at boundary nodes too, treating all nodes as unknowns. That approach keeps the standard RBF basis everywhere and modifies only the right-hand side at boundary evaluation points — it is simpler than the Hermite method. See Constructing an Operator Treating Boundary Nodes as Unknowns for the formulation.