design(grids): pick the non-box geometry seam — staircase mask vs phase-field vs Shortley–Weller

Start by reading src/Grids.jl, src/BlockForest.jl, src/boundaries.jl, src/operators/abstract.jl, src/multigrid.jl, and the cited DESIGN.md sections. Compare the four geometry seams against the listed invariants, then run the disk and sphere convergence studies and the materialized adjoint dot-product check. Done means a recommendation, one measured alternative, a passing adjoint identity, and a DESIGN.md §10 decision entry.

From Dennis Ogiermann, 2026-08-14: "There are papers from Pras Pathmanathan showing how to adopt FD to non-trivial geometries."

Where this sits relative to #61

#61 is the first-order answer; this is the second-order one. A binary active-cell mask makes a cell either in or out, which puts the domain boundary on a staircase — and a staircase boundary caps the geometric error at first order no matter how formally accurate the interior stencil is. For the hemisphere-in-a-box case that motivates #61 that may be perfectly adequate. For a curved tissue or ventricle boundary where the boundary flux is the physics, it is not, and the difference shows up as a convergence-rate ceiling rather than as a visible bug.

This issue is to decide, with numbers, which seam MatrixFreeOperators.jl adopts for non-box domains — not to build one speculatively.

What exists today

Nothing. There are exactly two AbstractGrid subtypes and both are axis-aligned boxes:

• CartesianGrid (src/Grids.jl:42) — uniform derived spacing (src/Grids.jl:67), no per-cell metric, no cell-volume or face-area accessor.
• BlockForest (src/BlockForest.jl:33) — a forest over the same box, whose constructor requires an exact root tiling (src/BlockForest.jl:49-54). AMR changes resolution, never participation.

Boundary conditions are per-dimension (low, high) face pairs attached to the grid (src/Grids.jl:56, src/boundaries.jl:12-71), implemented as ghost mirroring through _bc_sign (src/boundaries.jl:69-71) — six faces of a box, not per-cell and not per-face-fragment. The only mask in src/ is bounds-masking inside adjoint gathers (_maskedget, src/operators/abstract.jl:212), which is unrelated.

DESIGN.md is silent on this. It is not in §1 locked decisions, not in §9 forward-looking seams, not in §10 open decisions. The only relevant sentence is DESIGN.md:9-11 — "scoped initially to structured grids but built so any grid amenable to matrix-free application fits later" — plus the convention-based grid interface at DESIGN.md:113-114. #54 explicitly lists "unstructured, curvilinear, cut-cell, or embedded-boundary mesh framework" as out of scope, which is the right call for #54 and is exactly the boundary this issue is asking to revisit deliberately.

Candidate seams to evaluate

1. Staircase active-cell mask (#61). Cheapest, already scoped, keeps every stencil array-level. First-order geometric error.
2. Shortley–Weller / cut-cell one-sided stencils. Recovers second order for Dirichlet data on a curved boundary by shortening the arm that crosses it. Needs a per-cell distance-to-boundary and a per-cell stencil.
3. Generalized finite differences from a local neighbour set — weights solved per node against a Taylor system, typically via a pseudo-inverse. This is the family Dennis is pointing at. Most general, most expensive per cell, and the one that most clearly stops being "a stencil".
4. Level-set / embedded boundary with a signed-distance field carried as a Field and volume/aperture fractions derived from it. Composes best with AMR, and shares machinery with the geometry accessors #54 already proposes.

The tensions each option has to survive

These are the package invariants any answer must not quietly break, and they are the reason this is an investigation rather than a PR:

• Array-level leaf bodies. CLAUDE.md requires leaf bodies stay broadcast/slicing so they are device-agnostic and AD-friendly with no per-backend code. A per-cell variable stencil is a direct hit on that, and the KernelAbstractions @kernel escape hatch is supposed to be per-operator, not the default.
• Halo width. Every current stencil is compact — one ghost layer. Options 2–4 can widen the neighbour set, and halo width is baked into grid construction and into the exchange schedule.
• The adjoint is the gate. Whatever apply! does, apply_adjoint! must remain its exact discrete transpose. A shortened boundary arm has to be shortened identically on the way back, or ⟨Lx,y⟩ = ⟨x,Lᵀy⟩ fails. That test, not a plot, decides whether an implementation is correct.
• operator_diagonal becomes per-cell. Multigrid's smoothers and the coarsest-level LU (src/multigrid.jl:288) consume it, and the current implementations assume a diagonal that varies only by BC face.
• Conditioning. Small cut cells produce a small diagonal and a badly conditioned operator; that is the classic cut-cell failure mode and it lands on the Krylov path, not on the operator.
• Don't abstract ahead of evidence. DESIGN.md:862-870 — an abstraction is "deferred until a second topology backend … actually exists … Do not introduce it speculatively." This issue is allowed to produce a recommendation and one concrete implementation, not an AbstractGeometry layer.

Deliverable

• A short written comparison of the four options against the constraints above, naming the recommended one and why.
• A convergence study that makes the argument concrete rather than rhetorical: Poisson on a disk (2D) and a sphere (3D) with a manufactured analytic solution, measuring observed order for the staircase mask versus one chosen alternative on the same grids.
• The adjoint dot-product identity holding for whatever is implemented, on a materialized small problem.
• A DESIGN.md §10 entry recording the decision, so the next person does not re-open it from scratch.

Out of scope

• Building an unstructured or curvilinear mesh framework. The question is how far a structured grid can be pushed, not how to stop being one.
• Masked BlockForest or masked distributed slabs. Uniform CartesianGrid first, same boundary #61 and #54 drew.
• Moving or deforming geometry.

Open question — needs Dennis

The exact Pathmanathan references. A search turned up a generalized-finite-difference bidomain method on arbitrary irregular meshes that matches the description, but I could not confirm the author list, and there is more than one plausible line of work here. Worth getting the citations before spending a day on the wrong family.

_{🤖 Claude drew the staircase, measured it, and concluded stairs are not a curve.}