Manifold-mesh PDE support + extract_surface#202
Open
lmoresi wants to merge 8 commits into
Open
Conversation
KDTree.query reshaped the index array to 1-D when ``k=1`` but left the distance array 2-D. The docstring promises both as shape ``(n,)`` for ``k=1``. The asymmetry tripped naive ``indices[dists < tol]`` patterns in two consumers: ``_build_vertex_map`` and ``_build_dof_map`` (UW3 issue #197 — the second IndexError once the AttributeError there is worked around). One-line fix: also reshape ``d`` to match. Underworld development team with AI support from Claude Code
…ive stack UW3 was implicitly assuming ``mesh.dim == mesh.cdim`` in places that should size by the embedded coordinate space (gradients, flux vectors, constitutive tensors, vector field components). For volume meshes ``dim == cdim`` so every change here is a no-op; for manifold meshes (``SphericalManifold``: ``dim=2, cdim=3``) the gradient, flux, and constitutive tensors are genuinely cdim-sized in the embedding and need the right shape. Changes: - ``constitutive_models.py``: ``self.dim = u.mesh.cdim``. The conductivity / viscosity tensor multiplies the gradient (a cdim-vector from ``u.sym.jacobian(mesh.N)``), so the tensor sizes by cdim. - ``cython/petsc_generic_snes_solvers.pyx``: ``F1.reshape(cdim)`` in ``SNES_Scalar._setup_pointwise_functions`` (the flux vector is cdim-sized). - ``utilities/_jitextension.py``: gradient ``_ccodestr`` patches iterate to ``mesh.cdim`` (the JIT pointwise C code indexes ``petsc_u_x[fid * cdim + i]``). - ``discretisation_mesh_variables.py``: vector / tensor / sym_tensor MeshVariable shape and sym matrix size by cdim. Vtype inference accepts either dim or cdim as a vector match (since on a manifold the natural vector size is cdim). - ``swarm.py``: matching vtype inference change for SwarmVariable. - ``systems/ddt.py``: ``Symbolic_DDt``'s zero flux placeholder uses cdim. Vestigial NodalPointSwarm in ``SemiLagrangian_DDt`` removed — it was constructed but never read anywhere in the codebase (the actual SL trace-back uses ``uw.function.global_evaluate``). - ``systems/solvers.py``: ``SNES_AdvectionDiffusion``'s flux-DDt placeholder uses cdim. Each individual change preserves volume-mesh behaviour bit-exact. Underworld development team with AI support from Claude Code
Opt-in flag on ``SNES_Scalar`` that attaches a constant ``MatNullSpace`` to the Jacobian operator (and its transpose, and the preconditioner) at setup. Needed for scalar problems whose linear operator has a constant kernel and no Dirichlet BC to break it — pure-Neumann domains, closed manifolds (the spherical-surface Poisson case), periodic boxes. Without it the operator is singular and PETSc either diverges or returns an arbitrary solution within the constant null space. With it PETSc projects each KSP RHS onto the orthogonal complement of the nullspace and returns the minimum-norm (near-zero-mean) solution. Default ``False`` — existing scalar solvers untouched. Mirrors the standalone PR #201; once that merges to development the duplicate will be resolved on rebase. Bundled here so the surface- submesh branch is self-contained for testing. Implementation: - ``_petsc_use_constant_nullspace = False`` field in ``SNES_Scalar.__init__`` - ``petsc_use_constant_nullspace`` property + setter (setter invalidates ``is_setup`` so the rebuild attaches the nullspace) - ``_attach_constant_nullspace()`` helper: ``snes.setUp()``, then ``setNullSpace`` on operator + preconditioner + transposes - Single call site at the end of ``_setup_solver`` Underworld development team with AI support from Claude Code
Make the mesh layer's coordinate / navigation surfaces aware of the manifold case where ``dim < cdim`` (e.g. ``SphericalManifold``: ``dim=2, cdim=3``). All branches preserve volume-mesh behaviour (``dim == cdim``). ``coordinates.py``: - SPHERICAL CoordinateSystem dispatches on ``mesh.cdim == 3`` instead of ``mesh.dim == 3``. Spherical ``(r, θ, φ)`` is a function of the embedded Cartesian ``(x, y, z)`` — the right test is cdim. Fires for both 3-D volume shells and 2-manifold spheres. - ``_xRotN`` (the Cartesian-basis rotation matrix) sized by cdim. - ``SphericalCoordinateAccessor._dim = mesh.cdim`` so the ``r, θ, φ`` branch (vs ``r, θ`` polar) fires on any cdim=3 mesh. ``discretisation/discretisation_mesh.py``: - ``_build_kd_tree_index``: early-exit branch when ``dim != cdim`` builds a vertex-only kd-tree placeholder. ``DMPlexCreateSubmesh`` retains phantom depth-3 stratum points inside surface cell closures that break the volume-mesh ``closure[-n:]`` slicing (project workaround). The two-stage ``DMPlexFilter(depth, 2)`` strip in ``extract_surface`` makes this branch a no-op for the prototype path; kept as defence for future cd-1 entry points that don't strip the phantom DAG. - ``_mark_local_boundary_faces_inside_and_out``: short-circuit on cd-1. The 2-D perpendicular control-point trick doesn't apply on a curved cell, and "outside" of a closed manifold has no natural meaning. Sets the boundary-face kd-tree to None. - ``points_in_domain``: when the boundary-face tree is None, falls through to the closest-local-cell test (the right filter for on-manifold queries). - ``_get_closest_local_cells_internal``: skip the ``_test_if_points_in_cells_internal`` rejection on cd-1 (its 2-D perpendicular rule doesn't apply to a curved triangle); trust the centroid kd-tree result on convex manifolds with on-surface queries (the manifold-mesh contract). A proper manifold in-cell test (project query into cell tangent plane, then barycentric — reusing ``geometry_tools.distance_pointcloud_triangle``) is Phase-2 work and would lift the short-circuits. Underworld development team with AI support from Claude Code
``uw.meshing.SphericalManifold(radius, cellSize)`` builds a
triangulated sphere surface with ``dim=2, cdim=3`` (a 2-manifold
embedded in 3-D). The manifold sibling of ``SphericalShell``: same
gmsh sphere primitive, but the volume is removed and only the
surface is meshed.
Used for scalar PDEs on a closed manifold (Laplace–Beltrami, surface
diffusion, future shallow-water / surface-flow). The mesh inherits
``CoordinateSystemType.SPHERICAL`` so ``mesh.X.spherical.{r,theta,phi}``
work exactly as on ``SphericalShell``. The closed sphere's 3 rigid-
body rotation modes are populated on ``mesh._nullspace_rotations``
for future surface-vector-field solvers.
Two implementation notes:
1. ``dm_plex_gmsh_spacedim=3`` is set before ``DMPlexCreateFromFile``
so PETSc preserves the 3-D embedding when loading the 2-D mesh.
The standard ``_from_gmsh`` h5 round-trip drops ``coordDim``, so
the factory bypasses that and hands the DM directly to ``Mesh()``.
2. ``return_coords_to_bounds`` is wired to the closed-form sphere
projection ``x → R · x/|x|``. The standard SemiLagrangian
trace-back machinery already calls this after each Euler / RK2
substep, so on-surface advection lands launch points back on the
manifold automatically — no manifold-specific solver changes
needed at the trace-back-arithmetic layer.
Underworld development team with AI support from Claude Code
Surface submesh is the third UW3 submesh flavour alongside
``extract_region`` (subdomain via ``DMPlexFilter``) and
``coarsened_companion`` (resolution level via ``dm.refine()``).
``extract_surface(parent, label)`` pulls the cd-1 boundary marked
by a face label as a standalone ``uw.Mesh`` with
``dim = parent.dim − 1, cdim = parent.cdim``.
Two-stage PETSc construction:
sub1 = DMPlexCreateSubmesh(parent.dm, label, value,
marked_faces=True)
sub2 = DMPlexFilter(sub1, "depth", 2) # strip phantom DAG
Stage 1 retains an upward-DAG phantom depth-3 stratum (one point per
parent volume cell, celltype 12) so the resulting DM can still be
navigated as a slice of the parent. For a standalone manifold mesh
we don't use that linkage and the phantom points break naive
``closure[-n:]`` slicing across the codebase. Stage 2's
``DMPlexFilter`` on the depth=2 stratum drops the phantom, giving
a clean 3-stratum chart (depth = 0, 1, 2; celltypes [0, 1, 3]).
The two ``subpoint_is`` compose into a single surface → parent point
map.
Files (all under ``docs/examples/submesh_investigation/`` —
investigation tier, not promoted to ``Mesh`` API):
- ``surface_submesh_prototype.py`` — the ``extract_surface()`` impl
- ``example_surface_extraction.py`` — documented round-trip example
- ``probe_surface_extraction.py`` — Phase-0 PETSc sanity probe
- ``probe_surface_strip.py`` — phantom-stratum strip probe
- ``test_surface_submesh_contract.py`` — 6 contract tests covering
bad-label refusal, empty-stratum refusal, geometry sanity, scalar
round-trip, cell navigation returning valid cell IDs, and
``evaluate(expr, surface_pts)`` working to machine precision
The prototype workarounds ``Mesh._build_vertex_map`` (broken on
``development`` — UW3 issue #197) with an inline KDTree match. Once
#197 lands the workaround collapses to ``_build_vertex_map()``.
Underworld development team with AI support from Claude Code
…face) Update the submesh-solver-architecture design doc to describe the three submesh flavours now landed: | Flavour | Constructor | PETSc mechanism | |------------------|------------------------------|----------------------| | Subdomain | mesh.extract_region(label) | DMPlexFilter | | Resolution level | coarsened_companion(...) | dm.refine() hierarchy| | Surface (NEW) | extract_surface(mesh, label) | DMPlexCreateSubmesh | | | | + DMPlexFilter strip | The previous version named ``DMPlexCreateSubmesh`` as "Works but wrong dimension" — wrong for the air/rock subdomain case but exactly right for cd-1 surface extraction, the new flavour. Updated the PETSc-infrastructure table accordingly. Added a "Surface submesh: solver path" section describing the Phase-2 work needed to make solvers run on cd-1 meshes: - JIT audit for ``dim`` vs ``cdim`` - FE assembly with non-square Jacobian (PETSc handles this) - Coordinate-system symbols on a manifold - Surface outward normal symbol - BCs on a 1-D boundary curve (relevant for partial-surface submeshes) Status as of this PR: Tier-1 of that solver-path checklist (scalar Poisson on the closed sphere with constant-nullspace handling) is operational and verified against the Y_lm spherical-harmonic eigenfunctions to FE accuracy. Underworld development team with AI support from Claude Code
Contributor
There was a problem hiding this comment.
Pull request overview
This PR adds first-class support for solving scalar PDEs on manifold meshes (dim < cdim) and introduces new construction paths for those meshes, notably a direct SphericalManifold factory and a surface-extraction prototype workflow. It also extends the scalar SNES solver to optionally attach a constant nullspace for singular operators (e.g., closed-surface Poisson / pure-Neumann problems).
Changes:
- Plumbs
cdim(embedded coordinate dimension) through key variable, JIT, constitutive, and solver paths so flux/gradient/tensor shapes are consistent on manifold meshes. - Adds
uw.meshing.SphericalManifold(...)and updates spherical coordinate branching to work whencdim==3even ifdim==2. - Adds
SNES_Scalar.petsc_use_constant_nullspaceto stabilize singular scalar solves by attaching a constantMatNullSpace.
Reviewed changes
Copilot reviewed 18 out of 18 changed files in this pull request and generated 5 comments.
Show a summary per file
| File | Description |
|---|---|
| src/underworld3/utilities/_jitextension.py | JIT codegen patches gradient component wiring to iterate over mesh.cdim. |
| src/underworld3/systems/solvers.py | Sizes semi-Lagrangian flux placeholders by mesh.cdim. |
| src/underworld3/systems/ddt.py | Removes vestigial NodalPointSwarm allocation (_nswarm_psi). |
| src/underworld3/swarm.py | Extends SwarmVariable vtype inference to accept size==mesh.cdim for vectors/tensors. |
| src/underworld3/meshing/spherical.py | Adds SphericalManifold factory; loads gmsh with dm_plex_gmsh_spacedim=3; adds refinement snap + sphere projection callback. |
| src/underworld3/meshing/init.py | Re-exports SphericalManifold. |
| src/underworld3/discretisation/discretisation_mesh.py | Adds cd-1 navigation/inside-domain short-circuits and skips boundary-face control-point logic on dim!=cdim. |
| src/underworld3/discretisation/discretisation_mesh_variables.py | Infers vector/tensor types using dim or cdim, but stores VECTOR/TENSOR shapes using cdim. |
| src/underworld3/cython/petsc_generic_snes_solvers.pyx | Adds petsc_use_constant_nullspace flag; reshapes flux F1 by cdim; attaches constant nullspace during setup. |
| src/underworld3/coordinates.py | Keys spherical coordinate availability and basis rotation sizing off mesh.cdim (enables spherical accessors for dim=2, cdim=3). |
| src/underworld3/constitutive_models.py | Sizes constitutive tensors using mesh.cdim. |
| src/underworld3/ckdtree.pyx | Makes KDTree.query(k=1) distance output 1-D to match index shape. |
| docs/examples/submesh_investigation/test_surface_submesh_contract.py | Adds contract test script for the surface-submesh prototype. |
| docs/examples/submesh_investigation/surface_submesh_prototype.py | Adds prototype extract_surface() implementation and supporting helpers. |
| docs/examples/submesh_investigation/probe_surface_strip.py | Adds probe script for phantom-stratum stripping and IS composition. |
| docs/examples/submesh_investigation/probe_surface_extraction.py | Adds probe script exploring raw DMPlexCreateSubmesh behavior. |
| docs/examples/submesh_investigation/example_surface_extraction.py | Adds documented example of surface extraction + restrict/prolongate round-trip. |
| docs/developer/design/submesh-solver-architecture.md | Updates design doc to include “surface” as a third submesh flavour and documents the two-stage strip contract. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
Comment on lines
177
to
196
| if vtype == None: | ||
| # Note: on a cd-1 mesh (dim < cdim, e.g. SphericalManifold), | ||
| # vector fields are stored with ``cdim`` components in the | ||
| # embedded coordinate space (tangent-constrained 3-vectors | ||
| # on a 2-manifold) and the SwarmVariable coord cache is | ||
| # built with size == mesh.cdim. We accept either dim or | ||
| # cdim as a vector match; on volume meshes dim == cdim so | ||
| # the two branches coincide. | ||
| if isinstance(size, int) and size == 1: | ||
| vtype = uw.VarType.SCALAR | ||
| elif isinstance(size, int) and size == mesh.dim: | ||
| elif isinstance(size, int) and (size == mesh.dim or size == mesh.cdim): | ||
| vtype = uw.VarType.VECTOR | ||
| elif isinstance(size, tuple): | ||
| if size[0] == mesh.dim and size[1] == mesh.dim: | ||
| if ( | ||
| (size[0] == mesh.dim and size[1] == mesh.dim) | ||
| or (size[0] == mesh.cdim and size[1] == mesh.cdim) | ||
| ): | ||
| vtype = uw.VarType.TENSOR | ||
| else: | ||
| vtype = uw.VarType.MATRIX |
Comment on lines
233
to
+309
| @@ -280,19 +287,26 @@ def __init__( | |||
| else: | |||
| self._units = None | |||
|
|
|||
| # Component and shape handling | |||
| # Component and shape handling. | |||
| # Vector / tensor fields are sized by ``cdim`` (the embedded | |||
| # coordinate space) rather than ``dim`` (topological). For | |||
| # volume meshes dim == cdim so this is unchanged; for manifold | |||
| # meshes (e.g. SphericalManifold: dim=2, cdim=3) vectors are | |||
| # 3-component (tangent-constrained) and rank-2 tensors are | |||
| # 3x3 — matching the gradient / flux dimensions the JIT and | |||
| # solver layers expect. | |||
| if vtype == uw.VarType.SCALAR: | |||
| self.shape = (1, 1) | |||
| self.num_components = 1 | |||
| elif vtype == uw.VarType.VECTOR: | |||
| self.shape = (1, mesh.dim) | |||
| self.num_components = mesh.dim | |||
| self.shape = (1, mesh.cdim) | |||
| self.num_components = mesh.cdim | |||
| elif vtype == uw.VarType.TENSOR: | |||
| self.num_components = mesh.dim * mesh.dim | |||
| self.shape = (mesh.dim, mesh.dim) | |||
| self.num_components = mesh.cdim * mesh.cdim | |||
| self.shape = (mesh.cdim, mesh.cdim) | |||
| elif vtype == uw.VarType.SYM_TENSOR: | |||
| self.num_components = math.comb(mesh.dim + 1, 2) | |||
| self.shape = (mesh.dim, mesh.dim) | |||
| self.num_components = math.comb(mesh.cdim + 1, 2) | |||
| self.shape = (mesh.cdim, mesh.cdim) | |||
Comment on lines
+283
to
+293
| @timing.routine_timer_decorator | ||
| def SphericalManifold( | ||
| radius: float = 1.0, | ||
| cellSize: float = 0.1, | ||
| degree: int = 1, | ||
| qdegree: int = 2, | ||
| filename=None, | ||
| refinement=None, | ||
| gmsh_verbosity=0, | ||
| verbose=False, | ||
| ): |
Comment on lines
+421
to
+424
| PETSc.Options().setValue("dm_plex_gmsh_spacedim", 3) | ||
| surface_dm = PETSc.DMPlex().createFromFile( | ||
| uw_filename, interpolate=True, comm=PETSc.COMM_WORLD, | ||
| ) |
Comment on lines
1649
to
+1686
| @@ -1656,6 +1665,25 @@ class SNES_Scalar(SolverBaseClass): | |||
|
|
|||
| self.is_setup = False | |||
|
|
|||
| @property | |||
| def petsc_use_constant_nullspace(self): | |||
| """Whether to attach a constant MatNullSpace to the Jacobian. | |||
|
|
|||
| Set to ``True`` for scalar problems on closed manifolds (e.g. | |||
| Poisson on a ``SphericalManifold``) or fully-Neumann domains | |||
| where the linear operator has a constant kernel. PETSc projects | |||
| the right-hand side onto the orthogonal complement of the | |||
| nullspace and selects the minimum-norm solution from the | |||
| affine null-affine family, so the system becomes uniquely | |||
| solvable up to that nullspace. | |||
| """ | |||
| return self._petsc_use_constant_nullspace | |||
|
|
|||
| @petsc_use_constant_nullspace.setter | |||
| def petsc_use_constant_nullspace(self, value): | |||
| self._petsc_use_constant_nullspace = bool(value) | |||
| self.is_setup = False | |||
|
|
|||
Two scripts under docs/examples/manifold_pdes/ that demonstrate the manifold solver path end-to-end, with pyvista figures regenerated on each run (figures are .gitignored per repo convention). example_manifold_helmholtz.py: - (A) Spherical-harmonic spectrum recovery — solve (-Δ_S + I) T = Y_lm for 7 modes and verify each recovers T = Y_lm / (l(l+1) + 1) to FE accuracy. Relative L2 error 1e-3 to 2e-2 across the spectrum. - (B) h-convergence sweep on Y_10. cellSize 0.4 → 0.05 gives clean ~2nd-order convergence: 3.3e-3 → 8.4e-5 in L2. - (C) Closed-manifold Poisson with constant-nullspace handling (petsc_use_constant_nullspace=True). -Δ_S T = z recovers T = z/2 + const to FE accuracy (L_inf 3.7e-3 at cellSize 0.15). - pyvista rendering of (C) showing T = z/2 on the unit sphere. example_manifold_diffusion.py: - Pole-localised Gaussian → relaxation toward zero mean, 40 steps of forward time-stepping with Diffusion(theta=1.0). Tracks ||T||_2 vs the analytic Y_10 decay envelope exp(-l(l+1) κ t). - pyvista composite figure of 6 snapshots showing the Gaussian spreading and flattening. Configures the time-stepper with ``snes_type = "ksponly"`` so the linear time-step solve runs as a single KSP per step — no Newton line-search noise on a residual already at quadrature precision. Recommended pattern for linear time-dependent diffusion on manifold meshes. Both examples use SphericalManifold directly (the cleanest construction path — no submesh-extraction complexity). Underworld development team with AI support from Claude Code
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Adds first-class support for solving scalar PDEs on manifold meshes (
dim < cdim) and introduces two construction entry points:uw.meshing.SphericalManifold(direct 2-manifold sphere via gmsh) andextract_surface(parent, label)(codimension-1 submesh from a parent volume mesh). Tier-1 of the manifold-solver-path checklist is operational — scalar Poisson on a closed unit sphere recovers the analyticY_lmsolutions to FE accuracy.Status of related work
petsc_use_constant_nullspace) — duplicated here as commitdc6dacdso this branch is self-contained for testing. Once Add SNES_Scalar.petsc_use_constant_nullspace for singular Poisson #201 merges to development, the duplicate is resolved automatically on rebase.Mesh._build_vertex_mapbroken on development). Worked around in theextract_surfaceprototype via an inline KDTree match.Story (in commit order)
fix(kdtree)— 1-line shape consistency fix.KDTree.queryreshapesito 1-D fork=1but leavesd2-D; the asymmetry trippeddists < tolmasking in_build_vertex_mapand_build_dof_map. Half of [BUG] - extract_region's _build_vertex_map calls _get_kdtree() on CoordinateSystem, which has no such method #197.feat(cdim) plumbing— make the solver / JIT / constitutive / discretisation stack size bymesh.cdim(the embedded-coord dim) where it should, instead ofmesh.dim. Touches constitutive tensor sizing, F1 reshape in SNES_Scalar, JIT gradient_ccodestrpatches, MeshVariable + SwarmVariable shape inference, DDt flux placeholders. Every change is a no-op on volume meshes (dim == cdim); together they let cd-1 meshes' gradient / flux / tensor shapes thread through the assembly correctly. Also removes a vestigial NodalPointSwarm constructor inSemiLagrangian_DDtthat was assigned to_nswarm_psibut never read anywhere.feat(snes_scalar) petsc_use_constant_nullspace— mirrors PR Add SNES_Scalar.petsc_use_constant_nullspace for singular Poisson #201. Opt-in flag that attachesMatNullSpace().create(constant=True)to the Jacobian for singular Poisson on closed manifolds / fully-Neumann domains. DefaultFalseso existing solvers are untouched.feat(cdim) navigation + SPHERICAL— the mesh-level cd-1 awareness. Coord system fires SPHERICAL oncdim == 3(so the spherical accessor works on a 2-manifold sphere as it does on the 3-D shell). The four navigation paths (_build_kd_tree_index,_mark_local_boundary_faces_inside_and_out,_get_closest_local_cells_internal,points_in_domain) gain cd-1 short-circuits that avoid the 2-D perpendicular control-point trick (which doesn't apply on a curved cell) and trust the centroid kd-tree for on-manifold queries.feat(meshing) SphericalManifold— new public factory. Produces auw.Meshwithdim=2, cdim=3lying on the sphere to machine precision. InheritsCoordinateSystemType.SPHERICALsomesh.X.spherical.{r, theta, phi}work. Setsreturn_coords_to_boundsto the closed-form sphere projectionx → R·x/|x|, which the standard SemiLagrangian trace-back machinery calls after each Euler/RK2 substep — on-surface advection lands launch points back on the manifold automatically. (Subject to the DMSwarm-on-manifold question — see Phase-2 plan.)feat(submesh) extract_surface— investigation tier (docs/examples/submesh_investigation/). Two-stage construction:DMPlexCreateSubmesh(parent, label, value, marked_faces=True)produces a cd-1 DM that retains an upward-DAG phantom depth-3 stratum (one point per parent volume cell, celltype 12). The phantom points sit inside surface-cell closures and break naiveclosure[-n:]slicing. A second-stageDMPlexFilter(depth, 2)drops them, giving a clean 3-stratum chart. The twosubpoint_iscompose into a single surface → parent point map.Files: prototype, two probes, documented example, 6-contract test suite covering bad-label refusal, empty-stratum refusal, geometry sanity (
max|r-R| < 1e-10), bit-exact scalar round-trip, cell navigation returning valid cell IDs, andevaluate(expr, surface_pts)working to machine precision.docs(submesh)— design doc updated: two-flavour table becomes three,DMPlexCreateSubmeshrow marked as Works instead of "wrong dimension", new section on the surface-flavour design contract and Phase-2 solver-path blockers.Verification
(-Δ_S + α)T = zon the unit sphere recoversT = z/(2 + α)to L2 error ≈ 3e-4 at cellSize=0.1 (P2 elements).-Δ_S T = zwithpetsc_use_constant_nullspace=TruerecoversT = z/2 + const, monotone h-convergence:1.3e-2 → 5.4e-3 → 1.3e-3 → 3.5e-4halving cellSize from 0.4 → 0.05.extract_surface(SphericalShell(cellSize=0.2), "Upper")produces an equivalent mesh; same Helmholtz problem L2 error 1.4e-2 (10× worse than direct factory at matched cell count — element-shape quality, not assembly).Test plan
test_surface_submesh_contract.py(6/6 contracts pass)example_surface_extraction.py(4/4 sections pass)test_0001_meshes.py(19/19),test_0506_tensor_evaluate.py,test_1010_stokesCart.py(25/25 Stokes tests)What's NOT in this PR (Phase 2)
DMSwarmPIC_cooris allocated dim-sized;global_evaluateinternally adds particles atcdim-sized coords, which the swarm rejects). Two paths under discussion: a "swarm-lite" using the PETSc star-forest directly to scatter the minimal information, or fixing DMSwarm itself to carrycdim-sized coords on a cd-1 cell DM.utilities/geometry_tools.py(distance_pointcloud_triangle).Underworld development team with AI support from Claude Code