@@ -8,3 +8,4 @@ artifacts/
 dist/
 build/
 *.egg-info/
+data/ibamr_rbc_3d/
@@ -8,29 +8,28 @@
 **Meshfree geometry, RBF-FD, partition-of-unity methods, and PDE solvers for
 Python.**
 
-`kernelpack-python` provides Python implementations of meshfree geometry,
-scattered-node discretizations, and PDE solvers on embedded fixed domains. It
-combines reusable geometry and node-generation tools with standard and
-overlapped PHS+poly RBF-FD, weighted least squares, localized
-partition-of-unity approximations, and sparse elliptic and diffusion solvers.
-
-The package is the Python member of the KernelPack family. The companion
-[`kernelpack-matlab`](https://github.com/VarShankar/kernelpack-matlab)
-repository includes the moving-surface implementation accompanying the
-preprint [*A high-order, meshless, Lagrangian--Eulerian RBF-FD method for
-advection--diffusion--reaction on moving manifolds*](https://arxiv.org/abs/2608.19384)
-by Matthew Lowery, Grady B. Wright, and Varun Shankar. This Python release is
-focused on the fixed-domain numerical core and does not claim that
-moving-surface solver.
+`kernelpack-python` provides NumPy/SciPy implementations of meshfree geometry,
+scattered-node discretizations, and PDE solvers on fixed domains, moving
+domains, and evolving surfaces. It combines reusable geometry and node
+generation with standard and overlapped PHS+poly RBF-FD, weighted least
+squares, localized partition-of-unity approximations, and sparse iterative
+solvers.
 
-![Poisson solution and nodal error on an embedded domain](docs/readme_assets/poisson_solution.png)
+The moving-domain and moving-surface implementations follow the same published
+methods and numerical conventions as
+[`kernelpack-matlab`](https://github.com/VarShankar/kernelpack-matlab) and
+[`kernelpack-jax`](https://github.com/VarShankar/kernelpack-jax), while keeping
+the execution path in NumPy and SciPy.
+
+![Geometry and a boundary-refined meshfree node cloud](docs/readme_assets/geometry_domain.png)
 
-The figure shows an end-to-end Poisson solve on a geometry-defined scattered
-node set: the geometric model supplies the boundary and normals, the node
-generator fills the domain, and PHS+poly RBF-FD supplies the differential and
-boundary operators.
+The figure shows the geometry pipeline used by the PDE solvers: an embedded
+boundary and its normals define a level set, and the node generator fills the
+interior while retaining boundary and ghost-node structure.
 
 [Install](#installation) | [First solve](#first-solve) |
+[Workflows](#fixed-domain-workflows) | [Moving domains](#moving-domain-pdes) |
+[Moving surfaces](#moving-surface-pdes) |
 [Examples](#examples) | [Tests](#verification) |
 [Papers](#research-foundations) | [Citation](#citation)
 
@@ -43,10 +42,15 @@ want to:
 - generate scattered nodes and differential operators on embedded domains;
 - compare standard and overlapped local assembly;
 - solve elliptic and diffusion problems without constructing a volume mesh;
+- solve semi-Lagrangian ADR problems with moving embedded boundaries;
+- solve conservative ADR problems on stationary or evolving closed surfaces;
+- evolve point-cloud surfaces by mean curvature;
 - build localized PU or divergence-free RBF approximations; or
 - extend a tested, inspectable numerical research codebase.
 
-It is research software, not a general-purpose finite-element package.
+It is research software, not a general-purpose finite-element package. The
+implementation favors transparent SciPy sparse algebra and targeted Numba
+kernels; the JAX repository is the accelerator-oriented member of the family.
 
 ## At a glance
 
@@ -56,19 +60,28 @@ It is research software, not a general-purpose finite-element package.
 | Node generation | Seeded fixed- and variable-radius Poisson sampling in boxes, clipping by embedded geometry, boundary and ghost nodes, boundary-zone outer refinement, and dual node sets |
 | Local approximation | Centered and scaled Legendre polynomial bases, standard and overlapped PHS+poly RBF-FD, weighted-least-squares stencils, and local divergence-free PHS interpolation |
 | Fixed-domain solvers | Poisson, variable-coefficient and nonlinear variable-coefficient Poisson, BDF1--BDF3 diffusion, localized PU diffusion, and homogeneous or heterogeneous multispecies diffusion |
-| Numerical infrastructure | SciPy KD trees and sparse matrices, cached local polynomial templates, and targeted Numba kernels for repeated geometry, polynomial, and stencil calculations |
+| Moving-domain ADR | RK3 boundary motion, geometric reconstruction, carve/refill node updates, local PHS+Legendre departure-point interpolation, BDF1--BDF3 transport, and implicit diffusion--reaction solves in two and three dimensions |
+| Moving-surface ADR | Tangent-plane PHS+Legendre RBF-FD, cached-factor defect updates, adaptive surface hyperviscosity, spherical or toroidal geometry models, Lagrangian BDF1--BDF3 stepping, mass projection, and semi-Lagrangian history backfill |
+| Geometric surface evolution | Explicit or matrix-free semi-implicit tangent-plane RBF-FD mean-curvature-flow steps |
+| Numerical infrastructure | SciPy KD trees, batched LAPACK and sparse iterative solvers, plus cached Numba-parallel kernels for tangent-plane systems, stencil matvecs, SBF geometry, polynomial evaluation, and control-point selection |
 
 The main namespaces are `kernelpack.geometry`, `kernelpack.nodes`,
-`kernelpack.domain`, `kernelpack.poly`, `kernelpack.rbffd`,
+`kernelpack.domain`, `kernelpack.manifold`, `kernelpack.poly`, `kernelpack.rbffd`,
 `kernelpack.divfree`, and `kernelpack.solvers`.
 
+Repeated numerical kernels run outside the Python interpreter: Numba builds
+tangent-plane local systems in parallel and applies stencil operators inside
+GMRES, while SciPy supplies compiled KD-tree searches, batched LU solves, ILU,
+and sparse Krylov iterations. Python remains responsible for user callbacks,
+time-level orchestration, and topology-changing node-set updates.
+
 ## Requirements
 
 - Python 3.11 or newer
 - NumPy 2.0 or newer
 - SciPy 1.14 or newer
 - Numba 0.61 or newer
-- Matplotlib 3.9 or newer for examples and figures
+- Matplotlib 3.9 or newer is optional and used by examples and figure scripts
 
 ## Installation
 
@@ -153,9 +166,9 @@ plt.title("Poisson solution")
 plt.show()
 ```
 
-![Geometry-clipped interior, boundary, and ghost nodes](docs/readme_assets/geometry_domain.png)
+![Poisson solution and nodal error on an embedded domain](docs/readme_assets/poisson_solution.png)
 
-## Solver workflows
+## Fixed-domain workflows
 
 All solvers use a `DomainDescriptor`, so geometry, node generation, and local
 operator construction remain separate from the PDE definition. A target
@@ -172,16 +185,94 @@ reuses time-independent operators and preconditioners where possible. The PU
 variants localize approximation and support single- or multispecies diffusion,
 including distinct diffusivities by species.
 
+The same geometry and node descriptors feed every solver. This makes it
+possible to compare standard RBF-FD, overlapped RBF-FD, weighted least squares,
+and localized PU assembly without rewriting the domain construction or
+boundary callbacks.
+
 ![Diffusion solution, error, and time history](docs/readme_assets/diffusion_solution.png)
 
+## Moving-domain PDEs
+
+`MovingDomainADRSolver` advances
+
+\[
+\frac{D c}{D t}=\nu\Delta c+\lambda c+f
+\]
+
+on domains with moving embedded boundaries. RK3 advances boundary markers;
+the geometric model reconstructs the boundary; the background cloud is carved
+and refilled; and local PHS+Legendre interpolation supplies semi-Lagrangian
+BDF1--BDF3 history values. Diffusion and reaction are solved implicitly with
+overlapped RBF-FD, GMRES, and ILU. The same implementation supports the 2D
+moving-hole and 3D moving-cavity examples.
+
+```bash
+python examples/moving_domain_adr_example.py
+python examples/moving_domain_adr_3d_example.py
+```
+
+## Moving-surface PDEs
+
+The surface solver discretizes
+
+\[
+\partial^\bullet c+c\,\nabla_\Gamma\!\cdot\mathbf{w}
+=\nu\Delta_\Gamma c+R+f
+\]
+
+with target-centered tangent-plane PHS+Legendre RBF-FD. The implementation
+provides direct assembly, cached-factor defect updates, adaptive
+hyperviscosity, spherical and toroidal geometry models with analytic normals
+and quadrature, matrix-free implicit diffusion, BDF1--BDF3 material histories,
+mass projection, RK3 marker motion, local or SBF transfer, and semi-Lagrangian
+history backfill. The same operators also advance
+\(\mathbf{X}_t=\Delta_\Gamma\mathbf{X}\) by explicit or semi-implicit
+mean-curvature-flow steps.
+
+```bash
+python examples/stationary_surface_adr_example.py
+python examples/moving_surface_adr_example.py
+python examples/mean_curvature_flow_sphere_example.py
+python examples/mean_curvature_flow_ellipsoid_example.py
+python examples/surface_rearrangement_example.py
+```
+
+The three-dimensional red-blood-cell example replays the public IBAMR
+trajectory used by the moving-surface paper. The downloader verifies the
+release archive before installing the trajectory locally; the solver then
+reconstructs SBF geometry and normals, advances a source-free diffusing tracer
+with the same tangent-plane ADR machinery, enforces the quadrature mass law,
+and writes the final point cloud and concentration to `artifacts/`.
+
+```bash
+python scripts/download_rbc_capstone_data.py
+python examples/moving_surface_rbc_capstone.py
+```
+
+The default command processes the complete 401-frame trajectory. Use
+`--steps 2` only for a short installation check.
+
 ## Examples
 
-Complete convergence drivers live in [`examples`](examples):
+Complete workflows live in [`examples`](examples):
 
 | Goal | Example |
 | --- | --- |
+| Solve Poisson's equation | [`poisson_solver_example.py`](examples/poisson_solver_example.py) |
 | Verify a two-dimensional pure-Neumann Poisson solve | [`poisson_convergence_2d_neumann.py`](examples/poisson_convergence_2d_neumann.py) |
 | Verify a three-dimensional pure-Neumann Poisson solve | [`poisson_convergence_3d_neumann.py`](examples/poisson_convergence_3d_neumann.py) |
+| Solve variable and nonlinear variable-coefficient Poisson problems | [`variable_poisson_solver_example.py`](examples/variable_poisson_solver_example.py) |
+| Compare FD and localized-PU diffusion | [`diffusion_solver_example.py`](examples/diffusion_solver_example.py) |
+| Compare standard, heterogeneous, and PU multispecies diffusion | [`multispecies_diffusion_example.py`](examples/multispecies_diffusion_example.py) |
+| Solve moving-domain ADR in two dimensions | [`moving_domain_adr_example.py`](examples/moving_domain_adr_example.py) |
+| Exercise moving-domain ADR in three dimensions | [`moving_domain_adr_3d_example.py`](examples/moving_domain_adr_3d_example.py) |
+| Verify stationary-surface ADR | [`stationary_surface_adr_example.py`](examples/stationary_surface_adr_example.py) |
+| Verify moving-surface ADR on a breathing sphere | [`moving_surface_adr_example.py`](examples/moving_surface_adr_example.py) |
+| Evolve a sphere by mean curvature | [`mean_curvature_flow_sphere_example.py`](examples/mean_curvature_flow_sphere_example.py) |
+| Smooth an ellipsoid by mean curvature | [`mean_curvature_flow_ellipsoid_example.py`](examples/mean_curvature_flow_ellipsoid_example.py) |
+| Remap markers and semi-Lagrangian-backfill BDF history | [`surface_rearrangement_example.py`](examples/surface_rearrangement_example.py) |
+| Transport and diffuse a tracer on an IBAMR red-blood-cell trajectory | [`moving_surface_rbc_capstone.py`](examples/moving_surface_rbc_capstone.py) |
 
 Run a study from the repository root, for example:
 
@@ -220,11 +311,8 @@ your work.
 | Overlapped RBF-FD assembly (`kernelpack.rbffd.FDODiffOp`) | V. Shankar, [*The overlapped radial basis function-finite difference (RBF-FD) method: A generalization of RBF-FD*](https://doi.org/10.1016/j.jcp.2017.04.037), Journal of Computational Physics 342 (2017), 211--228 |
 | PHS geometric models and Poisson node generation (`kernelpack.geometry`, `kernelpack.nodes`) | V. Shankar, R. M. Kirby, and A. L. Fogelson, [*Robust node generation for mesh-free discretizations on irregular domains and surfaces*](https://doi.org/10.1137/17M114090X), SIAM Journal on Scientific Computing 40 (2018), A2584--A2608 |
 | Lower odd-degree PHS selection used by the solver stencil defaults | V. Shankar and A. L. Fogelson, [*Hyperviscosity-based stabilization for radial basis function-finite difference (RBF-FD) discretizations of advection-diffusion equations*](https://doi.org/10.1016/j.jcp.2018.06.036), Journal of Computational Physics 372 (2018), 616--639 |
-
-For moving-surface ADR, surface hyperviscosity, and the associated research
-drivers, see the public
-[`kernelpack-matlab`](https://github.com/VarShankar/kernelpack-matlab)
-release and its research-foundations table.
+| Moving-domain ADR | V. Shankar, G. B. Wright, and A. L. Fogelson, [*An efficient high-order meshless method for advection-diffusion equations on time-varying irregular domains*](https://doi.org/10.1016/j.jcp.2021.110633), Journal of Computational Physics 445 (2021), 110633 |
+| Moving-surface ADR | M. Lowery, G. B. Wright, and V. Shankar, [*A high-order, meshless, Lagrangian--Eulerian RBF-FD method for advection--diffusion--reaction on moving manifolds*](https://doi.org/10.48550/arXiv.2608.19384), arXiv:2608.19384 (2026) |
 
 ## Citation
 
 
@@ -0,0 +1,36 @@
+"""Shared domain and manufactured fields for the public PDE examples."""
+
+import numpy as np
+
+from kernelpack.domain import DomainDescriptor
+
+
+def build_square_domain(grid_size: int = 11) -> DomainDescriptor:
+    coordinates = np.linspace(-1.0, 1.0, grid_size)
+    xx, yy = np.meshgrid(coordinates, coordinates, indexing="ij")
+    points = np.column_stack((xx.ravel(), yy.ravel()))
+    boundary_mask = np.any(np.isclose(np.abs(points), 1.0), axis=1)
+    boundary = points[boundary_mask]
+    interior = points[~boundary_mask]
+    normals = np.where(np.isclose(np.abs(boundary), 1.0), np.sign(boundary), 0.0)
+    normals /= np.linalg.norm(normals, axis=1, keepdims=True)
+    ghosts = boundary + 0.25 * normals
+
+    domain = DomainDescriptor()
+    domain.set_nodes(interior, boundary, ghosts)
+    domain.set_normals(normals)
+    domain.set_sep_rad(float(coordinates[1] - coordinates[0]))
+    domain.build_structs()
+    return domain
+
+
+def exact_solution(time: float, points: np.ndarray) -> np.ndarray:
+    return np.exp(-time) * (points[:, 0] ** 2 + points[:, 1] ** 2)
+
+
+def zero_neumann(points: np.ndarray) -> np.ndarray:
+    return np.zeros(points.shape[0])
+
+
+def unit_dirichlet(points: np.ndarray) -> np.ndarray:
+    return np.ones(points.shape[0])
@@ -0,0 +1,45 @@
+"""Advance the same manufactured diffusion problem with FD and PU solvers."""
+
+import numpy as np
+
+from kernelpack import solvers
+
+from _common import build_square_domain, exact_solution, unit_dirichlet, zero_neumann
+
+
+def run_solver(solver, domain) -> np.ndarray:
+    dt = 0.02
+    diffusivity = 0.25
+    points = domain.get_int_bdry_nodes()
+    forcing = lambda nu, time, x: (
+        -np.exp(-time) * (x[:, 0] ** 2 + x[:, 1] ** 2)
+        - 4.0 * nu * np.exp(-time)
+    )
+    boundary = lambda alpha, beta, normals, time, xb: exact_solution(time, xb)
+    solver.init(domain, 3, dt, diffusivity)
+    solver.set_initial_state(exact_solution(0.0, points))
+    solver.bdf1_step(dt, forcing, zero_neumann, unit_dirichlet, boundary)
+    solver.bdf2_step(2.0 * dt, forcing, zero_neumann, unit_dirichlet, boundary)
+    return solver.bdf3_step(3.0 * dt, forcing, zero_neumann, unit_dirichlet, boundary)
+
+
+def main() -> None:
+    domain = build_square_domain()
+    points = domain.get_int_bdry_nodes()
+    fd_solution = run_solver(
+        solvers.DiffusionSolver(
+            lap_assembler="fd",
+            bc_assembler="fd",
+            lap_stencil="wls",
+            bc_stencil="wls",
+        ),
+        domain,
+    )
+    pu_solution = run_solver(solvers.PUDiffusionSolver(), domain)
+    truth = exact_solution(0.06, points)
+    print(f"FD relative l2 error: {np.linalg.norm(fd_solution - truth) / np.linalg.norm(truth):.6e}")
+    print(f"PU relative l2 error: {np.linalg.norm(pu_solution - truth) / np.linalg.norm(truth):.6e}")
+
+
+if __name__ == "__main__":
+    main()
@@ -0,0 +1,36 @@
+"""Smooth an ellipsoid with semi-implicit mean-curvature flow."""
+
+import numpy as np
+
+from kernelpack import geometry, manifold, solvers
+
+
+def main() -> None:
+    points = geometry.fibonacci_sphere(256) * np.array([1.35, 0.85, 0.65])
+    xi = 2
+    dt = 2.0e-4
+    step_count = 5
+    ell = manifold.surface_polynomial_degree(xi)
+    neighbors = manifold.build_surface_stencil_graph(
+        points, manifold.surface_stencil_size(ell)
+    )
+    normals = manifold.estimate_pca_normals(points, stencil_size=24)
+    initial_radius = np.linalg.norm(points - np.mean(points, axis=0), axis=1)
+    for _ in range(step_count):
+        operators, _ = manifold.assemble_tangent_plane_operators(
+            points, normals, neighbors, xi=xi
+        )
+        points, info = solvers.mean_curvature_flow_step(
+            points, normals, operators, dt
+        )
+        normals = manifold.estimate_pca_normals(
+            points, previous_normals=normals, stencil_size=24
+        )
+    final_radius = np.linalg.norm(points - np.mean(points, axis=0), axis=1)
+    print(f"initial radius range: {np.ptp(initial_radius):.6e}")
+    print(f"final radius range: {np.ptp(final_radius):.6e}")
+    print(f"linear relative residual: {info.relative_residual:.6e}")
+
+
+if __name__ == "__main__":
+    main()
@@ -0,0 +1,41 @@
+"""Evolve a sphere under ``X_t = Delta_Gamma X``."""
+
+import numpy as np
+
+from kernelpack import geometry, manifold, solvers
+
+
+def main() -> None:
+    node_count = 256
+    xi = 2
+    radius_zero = 1.4
+    dt = 2.0e-4
+    step_count = 10
+    directions = geometry.fibonacci_sphere(node_count)
+    points = radius_zero * directions
+    ell = manifold.surface_polynomial_degree(xi)
+    neighbors = manifold.build_surface_stencil_graph(
+        points, manifold.surface_stencil_size(ell)
+    )
+
+    for _ in range(step_count):
+        normals = points / np.linalg.norm(points, axis=1, keepdims=True)
+        operators, _ = manifold.assemble_tangent_plane_operators(
+            points, normals, neighbors, xi=xi
+        )
+        points, info = solvers.mean_curvature_flow_step(
+            points, normals, operators, dt
+        )
+
+    final_time = step_count * dt
+    exact_radius = np.sqrt(radius_zero**2 - 4.0 * final_time)
+    radius = np.linalg.norm(points, axis=1)
+    relative_error = np.linalg.norm(radius - exact_radius) / np.linalg.norm(
+        np.full_like(radius, exact_radius)
+    )
+    print(f"relative radius error: {relative_error:.6e}")
+    print(f"linear relative residual: {info.relative_residual:.6e}")
+
+
+if __name__ == "__main__":
+    main()