From 81808dfa6d1f4d741074d38fa779f0f8eaaa6f05 Mon Sep 17 00:00:00 2001
From: Tyagi <gthyagi@gmail.com>
Date: Fri, 24 Jul 2026 01:47:04 +0530
Subject: [PATCH 1/3] Fix 3D boundary reaction flux recovery

Recover pointwise traction from 3D rotated free-slip reactions by assembling a lumped triangular boundary mass and dividing nodal reactions by that mass. Support P1 traces directly and P2 traces through four-triangle facet subdivision, with global coordinate-based assembly for partition independence and a surface-mass-weighted mean gauge.

Add serial and MPI regressions using the Zhong et al. l=2 internal-load shell. The corrected surface and CMB topography coefficients reproduce the reference response and remain invariant across MPI partitions. Retain focused 2D traction and flux coverage.
---
 src/underworld3/utilities/boundary_flux.py    | 108 +++++++++++++++++-
 .../test_1064_rotated_freeslip_parallel.py    |  65 +++++++++++
 tests/test_1018_rotated_freeslip.py           |  69 +++++++++++
 3 files changed, 237 insertions(+), 5 deletions(-)

diff --git a/src/underworld3/utilities/boundary_flux.py b/src/underworld3/utilities/boundary_flux.py
index a36c54662..5c05fdcb1 100644
--- a/src/underworld3/utilities/boundary_flux.py
+++ b/src/underworld3/utilities/boundary_flux.py
@@ -140,12 +140,110 @@ def _desmear(solver, boundary, xs, R, mass, remove_mean, partial_reaction=True):
     csec = dm.getCoordinateSection()
     cvec = np.asarray(dm.getCoordinatesLocal().array).reshape(-1, dim)
     v0, v1 = dm.getDepthStratum(0)
+    if dim == 3:
+        if mass != "lumped":
+            raise NotImplementedError(
+                "3D boundary-flux recovery currently supports mass='lumped' only."
+            )
+
+        lsec = dm.getLocalSection()
+        f0, f1 = dm.getHeightStratum(1)
+        e0, e1 = dm.getDepthStratum(1)
+
+        def coord(q):
+            return _point_coord(dm, dim, cvec, csec, v0, v1, q)
+
+        nodeR = {_key(x, dim): float(r) for x, r in zip(xs, R)}
+        sis = _boundary_stratum_is(dm, solver.mesh, boundary)
+        facets = [] if not (sis and sis.getSize() > 0) else [
+            int(q) for q in sis.getIndices() if f0 <= int(q) < f1
+        ]
+        local_triangles = []
+
+        for facet in facets:
+            closure = [int(q) for q in dm.getTransitiveClosure(facet)[0]]
+            vertices = [q for q in closure if v0 <= q < v1]
+            edges = [q for q in closure if e0 <= q < e1]
+            if len(vertices) != 3:
+                raise NotImplementedError(
+                    "3D boundary-flux recovery currently requires triangular facets."
+                )
+
+            vertex_keys = {_key(coord(q), dim): q for q in vertices}
+            edge_midpoints = {}
+            for edge in edges:
+                if lsec.getFieldDof(edge, 0) <= 0:
+                    continue
+                edge_vertices = [
+                    int(q)
+                    for q in dm.getTransitiveClosure(edge)[0]
+                    if v0 <= int(q) < v1
+                ]
+                if len(edge_vertices) == 2:
+                    edge_key = frozenset(_key(coord(q), dim) for q in edge_vertices)
+                    edge_midpoints[edge_key] = _key(coord(edge), dim)
+
+            vk = list(vertex_keys)
+            if not edge_midpoints:
+                local_triangles.append(tuple(vk))
+                continue
+            if len(edge_midpoints) != 3:
+                raise NotImplementedError(
+                    "3D lumped recovery supports P1 or complete P2 triangular traces."
+                )
+
+            m01 = edge_midpoints[frozenset((vk[0], vk[1]))]
+            m12 = edge_midpoints[frozenset((vk[1], vk[2]))]
+            m20 = edge_midpoints[frozenset((vk[2], vk[0]))]
+            local_triangles.extend(
+                (
+                    (vk[0], m01, m20),
+                    (vk[1], m12, m01),
+                    (vk[2], m20, m12),
+                    (m01, m12, m20),
+                )
+            )
+
+        R_by = {}
+        for rank_values in comm.allgather(nodeR):
+            for key, value in rank_values.items():
+                R_by[key] = (
+                    R_by.get(key, 0.0) + value if partial_reaction else value
+                )
+
+        triangles = {}
+        for rank_triangles in comm.allgather(local_triangles):
+            for triangle in rank_triangles:
+                triangles[tuple(sorted(triangle))] = triangle
+
+        keys = sorted(R_by)
+        global_index = {key: i for i, key in enumerate(keys)}
+        reaction = np.array([R_by[key] for key in keys], dtype=float)
+        lumped_mass = np.zeros(len(keys), dtype=float)
+
+        for triangle in triangles.values():
+            a, b, c = (np.asarray(key, dtype=float) for key in triangle)
+            area = 0.5 * float(np.linalg.norm(np.cross(b - a, c - a)))
+            for key in triangle:
+                if key in global_index:
+                    lumped_mass[global_index[key]] += area / 3.0
+
+        missing = np.flatnonzero(lumped_mass <= 0.0)
+        if missing.size:
+            raise RuntimeError(
+                f"Boundary mass is zero at {missing.size} nodes on {boundary!r}."
+            )
+
+        flux = reaction / lumped_mass
+        if remove_mean:
+            mean = float(np.dot(flux, lumped_mass) / np.sum(lumped_mass))
+            flux -= mean
+        return np.array([flux[global_index[_key(x, dim)]] for x in xs])
+
     if dim != 2:
-        # no line-mass geometry yet in 3D → global-mean lumped fallback
-        tot = comm.allreduce(float(np.sum(R)), op=MPI.SUM)
-        cnt = comm.allreduce(int(len(R)), op=MPI.SUM)
-        m = tot / max(cnt, 1)
-        return np.asarray(R) - (m if remove_mean else 0.0)
+        raise NotImplementedError(
+            f"Boundary-flux recovery is not implemented for mesh dimension {dim}."
+        )
 
     e0, e1 = dm.getDepthStratum(1)
     def vcoord(q): return cvec[csec.getOffset(q) // dim]
diff --git a/tests/parallel/test_1064_rotated_freeslip_parallel.py b/tests/parallel/test_1064_rotated_freeslip_parallel.py
index 1dd9ff35f..43427cead 100644
--- a/tests/parallel/test_1064_rotated_freeslip_parallel.py
+++ b/tests/parallel/test_1064_rotated_freeslip_parallel.py
@@ -51,6 +51,9 @@
 # 3D spherical shell (free-slip both boundaries, all 3 rotation nullspace modes):
 # velocity L2. Recompute with `python <thisfile> spherical3d`.
 GOLDEN_SPHERICAL3D = 4.069689334228e-03
+# Zhong l=2 surface and CMB topography coefficients recovered from the 3D
+# rotated-constraint reaction. Recompute with `python <thisfile> spherical3d_topo`.
+GOLDEN_SPHERICAL3D_TOPO = (4.205055186099e-01, 7.704049016458e-01)
 # NONLINEAR (power-law) box with rotated free-slip through the manual Newton loop
 # (consistent tangent): (velocity L2, nonlinear iteration count — the number of
 # Newton increments solved, == len(ksp_its); this solve exits on the step-norm
@@ -158,6 +161,52 @@ def _spherical3d_diagnostics():
     return L2, int(info["ksp_its"]), int(info["ksp_reason"])
 
 
+def _spherical3d_topography_diagnostics():
+    """Zhong l=2 topography from the 3D rotated-constraint reaction."""
+    RI, RO, RINT = 0.55, 1.0, 0.775
+    mesh = uw.meshing.SphericalShellInternalBoundary(
+        radiusOuter=RO, radiusInternal=RINT, radiusInner=RI,
+        cellSize=0.25, qdegree=2, degree=1)
+    v = uw.discretisation.MeshVariable(
+        "Vst", mesh, mesh.dim, degree=2, continuous=True)
+    p = uw.discretisation.MeshVariable(
+        "Pst", mesh, 1, degree=1, continuous=True)
+    s = uw.systems.Stokes(mesh, velocityField=v, pressureField=p)
+    s.constitutive_model = uw.constitutive_models.ViscousFlowModel
+    s.constitutive_model.Parameters.shear_viscosity_0 = 1.0
+    theta = mesh.CoordinateSystem.xR[1]
+    unit_r = mesh.CoordinateSystem.unit_e_0
+    harmonic = sympy.assoc_legendre(2, 0, sympy.cos(theta))
+    s.add_natural_bc(harmonic * unit_r, "Internal")
+    s.add_rotated_freeslip_bc(0, "Upper", normal=unit_r)
+    s.add_rotated_freeslip_bc(0, "Lower", normal=-unit_r)
+    s.petsc_use_pressure_nullspace = True
+    s.petsc_options["snes_type"] = "ksponly"
+    s.tolerance = 1.0e-5
+    s.solve()
+
+    def harmonic_coefficient(boundary, response_sign):
+        xs, sigma_nn = s.boundary_normal_traction(boundary)
+        local = {
+            tuple(np.round(x, 12)): -float(value)
+            for x, value in zip(xs, sigma_nn)
+        }
+        samples = {}
+        for rank_values in uw.mpi.comm.allgather(local):
+            samples.update(rank_values)
+        coords = np.asarray(list(samples))
+        topography = np.asarray(list(samples.values()))
+        radii = np.linalg.norm(coords, axis=1)
+        harmonic_values = 0.5 * (3.0 * (coords[:, 2] / radii) ** 2 - 1.0)
+        return float(
+            response_sign
+            * np.dot(topography, harmonic_values)
+            / np.dot(harmonic_values, harmonic_values)
+        )
+
+    return harmonic_coefficient("Upper", 1.0), harmonic_coefficient("Lower", -1.0)
+
+
 def _annulus_fmg_diagnostics():
     """Annulus radial free-slip whose velocity block is CUSTOM GEOMETRIC FMG on a
     nested hierarchy (coarse annulus -> refine -> refine), via set_custom_fmg — no
@@ -360,6 +409,18 @@ def test_rotated_freeslip_spherical3d_partition_independent():
         f"{L2_ref} vs {L2}")
 
 
+def test_rotated_freeslip_spherical3d_topography_partition_independent():
+    """3D boundary-mass recovery gives partition-independent topography coefficients."""
+    surface, cmb = _spherical3d_topography_diagnostics()
+    surface_ref, cmb_ref = GOLDEN_SPHERICAL3D_TOPO
+    assert np.isclose(surface, surface_ref, rtol=1e-6, atol=0), (
+        f"3D surface topography differs serial vs np={uw.mpi.size}: "
+        f"{surface_ref} vs {surface}")
+    assert np.isclose(cmb, cmb_ref, rtol=1e-6, atol=0), (
+        f"3D CMB topography differs serial vs np={uw.mpi.size}: "
+        f"{cmb_ref} vs {cmb}")
+
+
 def test_rotated_freeslip_box_nonlinear_partition_independent():
     """NONLINEAR rotated free-slip is partition-independent: a power-law box solved by
     the manual Newton/Picard loop reproduces the serial velocity L2 and iteration count
@@ -428,6 +489,10 @@ def test_rotated_freeslip_dynamic_topography_partition_independent():
         _L2, _its, _reason = _spherical3d_diagnostics()
         if uw.mpi.rank == 0:
             print(f"DIAG_SPHERICAL3D {_L2:.12e} its={_its} reason={_reason}")
+    elif _kind == "spherical3d_topo":
+        _surface, _cmb = _spherical3d_topography_diagnostics()
+        if uw.mpi.rank == 0:
+            print(f"DIAG_SPHERICAL3D_TOPO {_surface:.12e} {_cmb:.12e}")
     else:
         _L2, _verr = _box_diagnostics()
         if uw.mpi.rank == 0:
diff --git a/tests/test_1018_rotated_freeslip.py b/tests/test_1018_rotated_freeslip.py
index 35f0da4be..27be50a53 100644
--- a/tests/test_1018_rotated_freeslip.py
+++ b/tests/test_1018_rotated_freeslip.py
@@ -120,6 +120,75 @@ def test_rotated_freeslip_spherical_shell_3d():
         assert rotfrac < 1e-8, f"rotation mode {k} gauge {rotfrac:.2e} not removed"
 
 
+def _spherical3d_reaction_topography():
+    """Zhong l=2 topography recovered directly from rotated constraint reactions."""
+
+    radius_inner = 0.55
+    radius_outer = 1.0
+    radius_internal = 0.775
+    mesh = uw.meshing.SphericalShellInternalBoundary(
+        radiusOuter=radius_outer,
+        radiusInternal=radius_internal,
+        radiusInner=radius_inner,
+        cellSize=0.25,
+        qdegree=2,
+        degree=1,
+    )
+    velocity = uw.discretisation.MeshVariable(
+        "U_topography_3d", mesh, mesh.dim, degree=2, continuous=True
+    )
+    pressure = uw.discretisation.MeshVariable(
+        "P_topography_3d", mesh, 1, degree=1, continuous=True
+    )
+    stokes = uw.systems.Stokes(
+        mesh, velocityField=velocity, pressureField=pressure
+    )
+    stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel
+    stokes.constitutive_model.Parameters.shear_viscosity_0 = 1.0
+    theta = mesh.CoordinateSystem.xR[1]
+    unit_r = mesh.CoordinateSystem.unit_e_0
+    harmonic = sympy.assoc_legendre(2, 0, sympy.cos(theta))
+    stokes.add_natural_bc(harmonic * unit_r, "Internal")
+    stokes.add_rotated_freeslip_bc(0, "Upper", normal=unit_r)
+    stokes.add_rotated_freeslip_bc(0, "Lower", normal=-unit_r)
+    stokes.petsc_use_pressure_nullspace = True
+    stokes.petsc_options["snes_type"] = "ksponly"
+    stokes.tolerance = 1.0e-5
+    stokes.solve()
+
+    def harmonic_coefficient(boundary, response_sign):
+        xs, sigma_nn = stokes.boundary_normal_traction(boundary)
+        local = {
+            tuple(np.round(x, 12)): -float(value)
+            for x, value in zip(xs, sigma_nn)
+        }
+        samples = {}
+        for rank_values in uw.mpi.comm.allgather(local):
+            samples.update(rank_values)
+        coords = np.asarray(list(samples))
+        topography = np.asarray(list(samples.values()))
+        radii = np.linalg.norm(coords, axis=1)
+        harmonic_values = 0.5 * (3.0 * (coords[:, 2] / radii) ** 2 - 1.0)
+        return float(
+            response_sign
+            * np.dot(topography, harmonic_values)
+            / np.dot(harmonic_values, harmonic_values)
+        )
+
+    surface = harmonic_coefficient("Upper", 1.0)
+    cmb = harmonic_coefficient("Lower", -1.0)
+    return surface, cmb
+
+
+@pytest.mark.level_2
+def test_rotated_freeslip_spherical3d_reaction_topography():
+    """3D reaction loads must be divided by boundary mass to recover pointwise stress."""
+
+    surface, cmb = _spherical3d_reaction_topography()
+    assert np.isclose(surface, 0.41920, rtol=0.10)
+    assert np.isclose(cmb, 0.77060, rtol=0.10)
+
+
 def test_rotated_freeslip_annulus_zero_leakage():
     """Annulus: per-node radial free-slip on both arcs → machine-zero v_r leakage,
     with the analytic radial normal."""

From bc62b8705292a24557ad455d4db63f96c34e110e Mon Sep 17 00:00:00 2001
From: Tyagi <gthyagi@gmail.com>
Date: Fri, 24 Jul 2026 08:28:58 +0530
Subject: [PATCH 2/3] Fix pointwise P2 surface reaction recovery

Replace the invalid four-subtriangle P2 lumping with the exact six-node triangular surface mass matrix and a sparse consistent solve. Add mass=auto so existing 2D and 3D P1 traces retain lumped recovery while 3D P2 selects the mathematically required consistent operator; reject explicit P2 lumping and unsupported 3D traces.

Add pointwise constant-flux P1/P2 tetrahedral-box regressions in serial and MPI, including two- and four-rank partitions. Extend the Zhong l=2 validation to compare all nodes, vertices, and edge midpoints independently, removing the prior all-node cancellation blind spot.

Document trace-order behavior, non-triangular limitations, analytic-normal partition guarantees, and the corrected 3D CBF behavior.
---
 docs/developer/CHANGELOG.md                   |   7 +
 .../cython/petsc_generic_snes_solvers.pyx     |  40 +++--
 src/underworld3/utilities/boundary_flux.py    | 161 +++++++++++++-----
 src/underworld3/utilities/rotated_bc.py       |  27 ++-
 .../test_1064_rotated_freeslip_parallel.py    |  79 ++++++---
 .../test_1065_boundary_flux_parallel.py       |  36 ++++
 tests/test_1018_rotated_freeslip.py           |  44 +++--
 tests/test_1019_boundary_flux.py              |  41 +++++
 8 files changed, 330 insertions(+), 105 deletions(-)

diff --git a/docs/developer/CHANGELOG.md b/docs/developer/CHANGELOG.md
index 86d694ef2..761513dcd 100644
--- a/docs/developer/CHANGELOG.md
+++ b/docs/developer/CHANGELOG.md
@@ -95,6 +95,13 @@ component exactly — correct on curved, tilted, and deformed boundaries (#293).
 - A general **consistent boundary flux (CBF) primitive** recovers boundary
   fluxes for any solver — surface heat flux / Nusselt number for scalar
   diffusion, boundary traction σ·n for Stokes (#294).
+- Three-dimensional CBF recovery now assembles the exact triangular trace mass:
+  P1 supports lumped or consistent recovery, while P2 uses the required
+  consistent six-node surface-mass solve. The default `mass="auto"` selects the
+  valid method; explicit P2 lumping and non-triangular 3D traces raise instead
+  of returning a non-pointwise reaction scaling. Strict MPI invariance of a
+  vector normal projection requires an analytic normal; geometric facet-normal
+  seam sensitivity is unchanged (#404).
 - Recorded as the preferred free-slip BC in the project guidance (#300);
   conda PETSc floor raised to ≥ 3.25 for FMG/rotation API consistency (#304).
 
diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx
index e9283e5af..cd448f458 100644
--- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx
+++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx
@@ -2294,7 +2294,7 @@ class SolverBaseClass(uw_object):
             self.dm.restoreLocalVec(xlocal)
             self.dm.restoreGlobalVec(gvec)
 
-    def boundary_flux(self, boundary, mass="lumped", remove_mean=False, normal=None):
+    def boundary_flux(self, boundary, mass="auto", remove_mean=False, normal=None):
         r"""Consistent boundary flux on ``boundary``, recovered from the essential-BC
         reaction of the last solve (the Consistent Boundary Flux method).
 
@@ -2304,15 +2304,21 @@ class SolverBaseClass(uw_object):
         number); for a **vector** solver the traction :math:`\sigma\cdot\hat n` (pass
         ``normal`` to get the scalar normal component :math:`\hat n\cdot\sigma\cdot\hat n`).
 
-        ``mass`` de-smears the nodal reaction with the ``"lumped"`` (diagonal, monotone —
-        no overshoot at a flux jump) or ``"consistent"`` boundary mass. ``remove_mean``
-        subtracts the boundary mean — leave ``False`` for a physical flux (the mean is
-        the Nusselt number); ``True`` gives a gauge-free field (e.g. dynamic topography).
-        Parallel-safe and partition-independent."""
+        ``mass`` de-smears the nodal reaction with ``"lumped"`` or ``"consistent"``
+        boundary mass. ``"auto"`` (default) selects lumped recovery for 2D traces and
+        3D P1 triangles, and the required consistent solve for 3D P2 triangles.
+        ``remove_mean`` subtracts the boundary mean — leave ``False`` for a physical
+        flux (the mean is the Nusselt number); ``True`` gives a gauge-free field.
+
+        Three-dimensional recovery supports triangular P1/P2 traces; quadrilateral
+        traces raise explicitly. Reaction and mass assembly are partition-independent.
+        For vector fluxes, supply an analytic ``normal`` when strict partition
+        independence of the normal projection is required; geometric facet-normal
+        averaging at partition seams has a small pre-existing partition sensitivity."""
         from underworld3.utilities.boundary_flux import boundary_flux as _bf
         return _bf(self, boundary, mass=mass, remove_mean=remove_mean, normal=normal)
 
-    def boundary_flux_field(self, boundary, field, mass="lumped",
+    def boundary_flux_field(self, boundary, field, mass="auto",
                             remove_mean=False, scale=1.0, normal=None):
         r"""Write the consistent boundary flux (see :meth:`boundary_flux`) onto a scalar
         MeshVariable ``field`` at the boundary nodes (interior untouched), multiplied by
@@ -5364,7 +5370,7 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
             J1.destroy(); J2.destroy()
         return rel > tol
 
-    def boundary_normal_traction(self, boundary, mass="lumped"):
+    def boundary_normal_traction(self, boundary, mass="auto"):
         r"""Return the boundary normal traction :math:`\sigma_{nn}` on a
         rotated-free-slip ``boundary`` as the constraint reaction from the last
         solve — the smooth, bounded quantity used for dynamic topography
@@ -5372,18 +5378,18 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
         prior :meth:`add_rotated_freeslip_bc` on ``boundary`` and a completed
         :meth:`solve`.
 
-        ``mass`` chooses the boundary-mass de-smear of the nodal reaction:
-        ``"lumped"`` (default) is monotone — it cannot overshoot where the traction
-        jumps (e.g. across a viscosity contrast), so it is the safe choice for driving
-        a free surface; ``"consistent"`` uses the full P2 line mass (marginally sharper
-        on smooth tractions, but overshoots at discontinuities)."""
+        ``mass="auto"`` (default) uses lumped recovery for 2D traces and 3D P1
+        triangles, and the required consistent surface-mass solve for 3D P2 triangles.
+        Explicit ``"lumped"`` and ``"consistent"`` choices remain available where
+        mathematically valid. Three-dimensional recovery currently supports triangular
+        P1/P2 traces only."""
         if self._rotated_freeslip_info is None:
             raise RuntimeError(
                 "boundary_normal_traction requires a completed rotated-free-slip solve.")
         from underworld3.utilities.rotated_bc import boundary_normal_traction as _bnt
         return _bnt(self, boundary, self._rotated_freeslip_info, mass=mass)
 
-    def dynamic_topography(self, boundary, field, buoyancy_scale=1.0, mass="lumped"):
+    def dynamic_topography(self, boundary, field, buoyancy_scale=1.0, mass="auto"):
         r"""Write the dynamic topography
         :math:`h = -(\sigma_{nn}-\overline{\sigma_{nn}})/(\Delta\rho\,g)` on a
         rotated-free-slip ``boundary`` onto a scalar MeshVariable ``field``, from the
@@ -5393,9 +5399,9 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
         pass it here after each :meth:`solve`; its boundary nodes are filled and the
         interior left untouched.
 
-        ``buoyancy_scale`` is :math:`\Delta\rho\,g` (traction → length). ``mass`` selects
-        the recovery de-smear (``"lumped"`` default is monotone — no overshoot at a
-        stress jump — and is the safe choice for a free surface). Requires a prior
+        ``buoyancy_scale`` is :math:`\Delta\rho\,g` (traction → length).
+        ``mass="auto"`` selects lumped recovery where valid and the consistent
+        surface-mass solve for 3D P2 triangles. Requires a prior
         :meth:`add_rotated_freeslip_bc` on ``boundary`` and a completed :meth:`solve`."""
         if self._rotated_freeslip_info is None:
             raise RuntimeError(
diff --git a/src/underworld3/utilities/boundary_flux.py b/src/underworld3/utilities/boundary_flux.py
index 5c05fdcb1..9be243592 100644
--- a/src/underworld3/utilities/boundary_flux.py
+++ b/src/underworld3/utilities/boundary_flux.py
@@ -15,8 +15,10 @@
 assembled here in ``_desmear`` by SUMMING each rank's partial contribution by coordinate
 (the same rock-solid gather used for the boundary mass — no hand-rolled global assembly).
 
-``mass="lumped"`` (default) uses the diagonal boundary mass: being an M-matrix it cannot
-overshoot where the flux jumps (no Gibbs wiggle) and is a purely local division.
+``mass="auto"`` (default) uses a diagonal lumped mass where the trace basis admits
+positive row sums (the 2D P2 line trace and the 3D P1 triangle trace), and the consistent
+mass otherwise. A 3D P2 triangle has exactly zero row sum at every vertex, so its
+pointwise recovery requires the consistent surface-mass solve.
 ``remove_mean=False`` (default) keeps the physical mean flux (the Nusselt number);
 set ``remove_mean=True`` for a gauge-free field (e.g. dynamic topography).
 """
@@ -24,6 +26,29 @@
 from mpi4py import MPI
 
 
+# M_e = (area / 12) * _P1_TRIANGLE_MASS.
+_P1_TRIANGLE_MASS = np.array(
+    (
+        (2.0, 1.0, 1.0),
+        (1.0, 2.0, 1.0),
+        (1.0, 1.0, 2.0),
+    )
+)
+
+# M_e = (area / 180) * _P2_TRIANGLE_MASS. Node order: vertices
+# (0, 1, 2), then edge nodes (01, 12, 20).
+_P2_TRIANGLE_MASS = np.array(
+    (
+        (6.0, -1.0, -1.0, 0.0, -4.0, 0.0),
+        (-1.0, 6.0, -1.0, 0.0, 0.0, -4.0),
+        (-1.0, -1.0, 6.0, -4.0, 0.0, 0.0),
+        (0.0, 0.0, -4.0, 32.0, 16.0, 16.0),
+        (-4.0, 0.0, 0.0, 16.0, 32.0, 16.0),
+        (0.0, -4.0, 0.0, 16.0, 16.0, 32.0),
+    )
+)
+
+
 def _key(c, dim):
     return tuple(round(float(t), 9) for t in np.asarray(c).ravel()[:dim])
 
@@ -140,13 +165,11 @@ def _desmear(solver, boundary, xs, R, mass, remove_mean, partial_reaction=True):
     csec = dm.getCoordinateSection()
     cvec = np.asarray(dm.getCoordinatesLocal().array).reshape(-1, dim)
     v0, v1 = dm.getDepthStratum(0)
+    if mass not in ("auto", "lumped", "consistent"):
+        raise ValueError("mass must be 'auto', 'lumped', or 'consistent'.")
     if dim == 3:
-        if mass != "lumped":
-            raise NotImplementedError(
-                "3D boundary-flux recovery currently supports mass='lumped' only."
-            )
-
         lsec = dm.getLocalSection()
+        ncomp = lsec.getFieldComponents(0)
         f0, f1 = dm.getHeightStratum(1)
         e0, e1 = dm.getDepthStratum(1)
 
@@ -158,7 +181,7 @@ def coord(q):
         facets = [] if not (sis and sis.getSize() > 0) else [
             int(q) for q in sis.getIndices() if f0 <= int(q) < f1
         ]
-        local_triangles = []
+        local_elements = []
 
         for facet in facets:
             closure = [int(q) for q in dm.getTransitiveClosure(facet)[0]]
@@ -168,12 +191,22 @@ def coord(q):
                 raise NotImplementedError(
                     "3D boundary-flux recovery currently requires triangular facets."
                 )
+            if lsec.getFieldDof(facet, 0) > 0:
+                raise NotImplementedError(
+                    "3D boundary-flux recovery supports P1 or P2 triangular traces."
+                )
 
-            vertex_keys = {_key(coord(q), dim): q for q in vertices}
+            vertex_coords = [np.asarray(coord(q), dtype=float) for q in vertices]
+            vertex_keys = [_key(value, dim) for value in vertex_coords]
             edge_midpoints = {}
             for edge in edges:
-                if lsec.getFieldDof(edge, 0) <= 0:
+                edge_dof = lsec.getFieldDof(edge, 0)
+                if edge_dof <= 0:
                     continue
+                if edge_dof != ncomp:
+                    raise NotImplementedError(
+                        "3D boundary-flux recovery supports P1 or P2 triangular traces."
+                    )
                 edge_vertices = [
                     int(q)
                     for q in dm.getTransitiveClosure(edge)[0]
@@ -183,24 +216,24 @@ def coord(q):
                     edge_key = frozenset(_key(coord(q), dim) for q in edge_vertices)
                     edge_midpoints[edge_key] = _key(coord(edge), dim)
 
-            vk = list(vertex_keys)
+            a, b, c = vertex_coords
+            area = 0.5 * float(np.linalg.norm(np.cross(b - a, c - a)))
             if not edge_midpoints:
-                local_triangles.append(tuple(vk))
+                local_elements.append((1, tuple(vertex_keys), area))
                 continue
             if len(edge_midpoints) != 3:
                 raise NotImplementedError(
-                    "3D lumped recovery supports P1 or complete P2 triangular traces."
+                    "3D boundary-flux recovery supports P1 or complete P2 triangular traces."
                 )
 
-            m01 = edge_midpoints[frozenset((vk[0], vk[1]))]
-            m12 = edge_midpoints[frozenset((vk[1], vk[2]))]
-            m20 = edge_midpoints[frozenset((vk[2], vk[0]))]
-            local_triangles.extend(
+            m01 = edge_midpoints[frozenset((vertex_keys[0], vertex_keys[1]))]
+            m12 = edge_midpoints[frozenset((vertex_keys[1], vertex_keys[2]))]
+            m20 = edge_midpoints[frozenset((vertex_keys[2], vertex_keys[0]))]
+            local_elements.append(
                 (
-                    (vk[0], m01, m20),
-                    (vk[1], m12, m01),
-                    (vk[2], m20, m12),
-                    (m01, m12, m20),
+                    2,
+                    tuple(vertex_keys) + (m01, m12, m20),
+                    area,
                 )
             )
 
@@ -211,32 +244,72 @@ def coord(q):
                     R_by.get(key, 0.0) + value if partial_reaction else value
                 )
 
-        triangles = {}
-        for rank_triangles in comm.allgather(local_triangles):
-            for triangle in rank_triangles:
-                triangles[tuple(sorted(triangle))] = triangle
+        elements = {}
+        for rank_elements in comm.allgather(local_elements):
+            for order, nodes, area in rank_elements:
+                elements[(order, tuple(sorted(nodes)))] = (order, nodes, area)
+
+        orders = {order for order, _nodes, _area in elements.values()}
+        if len(orders) != 1:
+            raise RuntimeError(
+                f"Expected one trace order on boundary {boundary!r}, found {sorted(orders)}."
+            )
+        order = orders.pop()
+        if mass == "auto":
+            mass = "consistent" if order == 2 else "lumped"
+        if order == 2 and mass == "lumped":
+            raise ValueError(
+                "A 3D P2 triangular trace has zero row-sum mass at its vertices; "
+                "use mass='consistent' for pointwise boundary-flux recovery."
+            )
 
         keys = sorted(R_by)
         global_index = {key: i for i, key in enumerate(keys)}
         reaction = np.array([R_by[key] for key in keys], dtype=float)
-        lumped_mass = np.zeros(len(keys), dtype=float)
-
-        for triangle in triangles.values():
-            a, b, c = (np.asarray(key, dtype=float) for key in triangle)
-            area = 0.5 * float(np.linalg.norm(np.cross(b - a, c - a)))
-            for key in triangle:
-                if key in global_index:
-                    lumped_mass[global_index[key]] += area / 3.0
-
-        missing = np.flatnonzero(lumped_mass <= 0.0)
-        if missing.size:
-            raise RuntimeError(
-                f"Boundary mass is zero at {missing.size} nodes on {boundary!r}."
+        if mass == "lumped":
+            boundary_mass = np.zeros(len(keys), dtype=float)
+            for _order, nodes, area in elements.values():
+                for key in nodes:
+                    boundary_mass[global_index[key]] += area / 3.0
+            missing = np.flatnonzero(boundary_mass <= 0.0)
+            if missing.size:
+                raise RuntimeError(
+                    f"Boundary mass is zero at {missing.size} nodes on {boundary!r}."
+                )
+            flux = reaction / boundary_mass
+        elif mass == "consistent":
+            from scipy.sparse import coo_matrix
+            from scipy.sparse.linalg import spsolve
+
+            rows = []
+            cols = []
+            values = []
+            reference_mass = (
+                _P1_TRIANGLE_MASS if order == 1 else _P2_TRIANGLE_MASS
             )
-
-        flux = reaction / lumped_mass
+            mass_scale = 12.0 if order == 1 else 180.0
+            for _order, nodes, area in elements.values():
+                indices = [global_index[key] for key in nodes]
+                element_mass = (area / mass_scale) * reference_mass
+                for i, row in enumerate(indices):
+                    for j, col in enumerate(indices):
+                        rows.append(row)
+                        cols.append(col)
+                        values.append(element_mass[i, j])
+            surface_mass = coo_matrix(
+                (values, (rows, cols)), shape=(len(keys), len(keys))
+            ).tocsr()
+            surface_mass.sum_duplicates()
+            flux = np.asarray(spsolve(surface_mass, reaction), dtype=float)
+            boundary_mass = np.asarray(
+                surface_mass @ np.ones(len(keys), dtype=float)
+            )
+            if not np.all(np.isfinite(flux)):
+                raise RuntimeError(
+                    f"Consistent boundary-mass solve failed on boundary {boundary!r}."
+                )
         if remove_mean:
-            mean = float(np.dot(flux, lumped_mass) / np.sum(lumped_mass))
+            mean = float(np.dot(flux, boundary_mass) / np.sum(boundary_mass))
             flux -= mean
         return np.array([flux[global_index[_key(x, dim)]] for x in xs])
 
@@ -244,6 +317,8 @@ def coord(q):
         raise NotImplementedError(
             f"Boundary-flux recovery is not implemented for mesh dimension {dim}."
         )
+    if mass == "auto":
+        mass = "lumped"
 
     e0, e1 = dm.getDepthStratum(1)
     def vcoord(q): return cvec[csec.getOffset(q) // dim]
@@ -296,7 +371,7 @@ def vcoord(q): return cvec[csec.getOffset(q) // dim]
     return np.array([sig[gi[_key(x, dim)]] for x in xs])
 
 
-def boundary_flux(solver, boundary, mass="lumped", remove_mean=False, normal=None):
+def boundary_flux(solver, boundary, mass="auto", remove_mean=False, normal=None):
     """See ``SolverBaseClass.boundary_flux``. Returns ``(xs, flux)`` for this rank's
     boundary nodes; scalar solver → normal flux, vector solver → traction (or its normal
     component if ``normal`` is given)."""
@@ -351,7 +426,7 @@ def write_boundary_scalar_field(solver, field, value_by_key, dim):
     return field
 
 
-def boundary_flux_field(solver, boundary, field, mass="lumped",
+def boundary_flux_field(solver, boundary, field, mass="auto",
                         remove_mean=False, scale=1.0, normal=None):
     r"""See ``SolverBaseClass.boundary_flux_field`` (the documented entry point;
     this free function is its implementation and shares its name). Writes
diff --git a/src/underworld3/utilities/rotated_bc.py b/src/underworld3/utilities/rotated_bc.py
index 5ea39acae..1bb9c9c91 100644
--- a/src/underworld3/utilities/rotated_bc.py
+++ b/src/underworld3/utilities/rotated_bc.py
@@ -1042,7 +1042,7 @@ def _mode_satisfies_constraints(solver, Q, normal_rows, tg, tol=1e-8):
     return viol < tol
 
 
-def boundary_normal_traction(solver, boundary, solve_result, mass="lumped"):
+def boundary_normal_traction(solver, boundary, solve_result, mass="auto"):
     """Boundary normal traction σ_nn on `boundary` from the constraint reaction of the
     last rotated-free-slip solve (``solve_result`` is the dict returned by
     ``solve_rotated_freeslip`` / ``solve_rotated_freeslip_nonlinear`` — see their
@@ -1055,22 +1055,19 @@ def boundary_normal_traction(solver, boundary, solve_result, mass="lumped"):
     rotated frame's normal row) is corner-correct — at a node shared with another
     rotated-free-slip boundary the rotated frame's first row is a mix of both walls'
     normals, but n̂·r_c is the true normal traction for this boundary. The pointwise
-    σ_nn is the boundary-mass de-smear of R (2D).
+    σ_nn is the boundary-mass de-smear of R.
 
     ``mass`` selects the de-smear:
-      * ``"lumped"`` (default) — the diagonal (row-sum) boundary mass. Being an M-matrix
-        it CANNOT overshoot at a stress discontinuity (no Gibbs wiggle where the traction
-        jumps, e.g. across a viscosity contrast), it is a purely local division (no global
-        mass solve → trivially parallel), and it is marginally more accurate than the
-        consistent mass on SolCx. Recommended for driving a free surface, where an
-        overshoot at a sharp feature injects a spurious surface-velocity pulse.
-      * ``"consistent"`` — the full consistent P2 line mass. Marginally sharper on smooth
-        tractions but overshoots at discontinuities.
+      * ``"auto"`` (default) — lumped for 2D traces and 3D P1 triangles, consistent for
+        3D P2 triangles.
+      * ``"lumped"`` — the diagonal row-sum mass. It is monotone for supported traces,
+        but invalid for 3D P2 triangles because their vertex row sums are exactly zero.
+      * ``"consistent"`` — the full trace mass. Required for pointwise 3D P2 recovery.
 
     Parallel-safe: r_c is scattered to a local vector (ghosts included) and read by LOCAL
     section offset; the boundary mass is assembled globally by a coordinate-keyed
     allgather of the boundary elements, so every rank produces the same de-smear and the
-    mean-removal gauge is global.
+    mean-removal gauge is global. In 3D, only triangular P1/P2 traces are supported.
     """
     dm = solver.dm
     dim = solver.mesh.dim
@@ -1119,12 +1116,12 @@ def boundary_normal_traction(solver, boundary, solve_result, mass="lumped"):
 
 
 def dynamic_topography_field(solver, boundary, solve_result, field,
-                             buoyancy_scale=1.0, mass="lumped"):
+                             buoyancy_scale=1.0, mass="auto"):
     """Populate a scalar MeshVariable ``field`` with the dynamic topography
     :math:`h = -(\\sigma_{nn}-\\overline{\\sigma_{nn}})/(\\Delta\\rho\\,g)` on ``boundary``,
-    recovered from the rotated-free-slip constraint reaction (lumped by default —
-    monotone, no Gibbs overshoot at a stress jump). Interior nodes are left untouched.
-    Returns ``field``.
+    recovered from the rotated-free-slip constraint reaction. ``mass="auto"`` uses
+    lumped recovery where valid and the consistent surface mass for 3D P2 triangles.
+    Interior nodes are left untouched. Returns ``field``.
 
     This is the hand-off to the free-surface machinery: the 3-number topography
     integrator drives node motion from a surface field, so σ_nn is written onto the
diff --git a/tests/parallel/test_1064_rotated_freeslip_parallel.py b/tests/parallel/test_1064_rotated_freeslip_parallel.py
index 43427cead..4c882c389 100644
--- a/tests/parallel/test_1064_rotated_freeslip_parallel.py
+++ b/tests/parallel/test_1064_rotated_freeslip_parallel.py
@@ -51,9 +51,17 @@
 # 3D spherical shell (free-slip both boundaries, all 3 rotation nullspace modes):
 # velocity L2. Recompute with `python <thisfile> spherical3d`.
 GOLDEN_SPHERICAL3D = 4.069689334228e-03
-# Zhong l=2 surface and CMB topography coefficients recovered from the 3D
-# rotated-constraint reaction. Recompute with `python <thisfile> spherical3d_topo`.
-GOLDEN_SPHERICAL3D_TOPO = (4.205055186099e-01, 7.704049016458e-01)
+# Zhong l=2 topography coefficients recovered from the 3D rotated-constraint
+# reaction: surface (all, vertices, midpoints), CMB (all, vertices, midpoints).
+# Recompute with `python <thisfile> spherical3d_topo`.
+GOLDEN_SPHERICAL3D_TOPO = (
+    4.149689252074e-01,
+    3.952301937705e-01,
+    4.215939953379e-01,
+    7.932177563075e-01,
+    8.426041179682e-01,
+    7.762363224500e-01,
+)
 # NONLINEAR (power-law) box with rotated free-slip through the manual Newton loop
 # (consistent tangent): (velocity L2, nonlinear iteration count — the number of
 # Newton increments solved, == len(ksp_its); this solve exits on the step-norm
@@ -185,7 +193,19 @@ def _spherical3d_topography_diagnostics():
     s.tolerance = 1.0e-5
     s.solve()
 
-    def harmonic_coefficient(boundary, response_sign):
+    dm = mesh.dm
+    csec = dm.getCoordinateSection()
+    cvec = np.asarray(dm.getCoordinatesLocal().array).reshape(-1, mesh.dim)
+    vertex_start, vertex_end = dm.getDepthStratum(0)
+    local_vertex_keys = {
+        tuple(np.round(cvec[csec.getOffset(point) // mesh.dim], 12))
+        for point in range(vertex_start, vertex_end)
+    }
+    vertex_keys = set()
+    for rank_keys in uw.mpi.comm.allgather(local_vertex_keys):
+        vertex_keys.update(rank_keys)
+
+    def harmonic_coefficients(boundary, response_sign):
         xs, sigma_nn = s.boundary_normal_traction(boundary)
         local = {
             tuple(np.round(x, 12)): -float(value)
@@ -198,13 +218,21 @@ def harmonic_coefficient(boundary, response_sign):
         topography = np.asarray(list(samples.values()))
         radii = np.linalg.norm(coords, axis=1)
         harmonic_values = 0.5 * (3.0 * (coords[:, 2] / radii) ** 2 - 1.0)
-        return float(
-            response_sign
-            * np.dot(topography, harmonic_values)
-            / np.dot(harmonic_values, harmonic_values)
-        )
+        is_vertex = np.array([key in vertex_keys for key in samples], dtype=bool)
+
+        def fit(mask):
+            return float(
+                response_sign
+                * np.dot(topography[mask], harmonic_values[mask])
+                / np.dot(harmonic_values[mask], harmonic_values[mask])
+            )
 
-    return harmonic_coefficient("Upper", 1.0), harmonic_coefficient("Lower", -1.0)
+        return fit(np.ones(len(coords), dtype=bool)), fit(is_vertex), fit(~is_vertex)
+
+    return (
+        *harmonic_coefficients("Upper", 1.0),
+        *harmonic_coefficients("Lower", -1.0),
+    )
 
 
 def _annulus_fmg_diagnostics():
@@ -411,14 +439,22 @@ def test_rotated_freeslip_spherical3d_partition_independent():
 
 def test_rotated_freeslip_spherical3d_topography_partition_independent():
     """3D boundary-mass recovery gives partition-independent topography coefficients."""
-    surface, cmb = _spherical3d_topography_diagnostics()
-    surface_ref, cmb_ref = GOLDEN_SPHERICAL3D_TOPO
-    assert np.isclose(surface, surface_ref, rtol=1e-6, atol=0), (
-        f"3D surface topography differs serial vs np={uw.mpi.size}: "
-        f"{surface_ref} vs {surface}")
-    assert np.isclose(cmb, cmb_ref, rtol=1e-6, atol=0), (
-        f"3D CMB topography differs serial vs np={uw.mpi.size}: "
-        f"{cmb_ref} vs {cmb}")
+    coefficients = _spherical3d_topography_diagnostics()
+    labels = (
+        "surface all",
+        "surface vertices",
+        "surface midpoints",
+        "CMB all",
+        "CMB vertices",
+        "CMB midpoints",
+    )
+    for label, value, reference in zip(
+        labels, coefficients, GOLDEN_SPHERICAL3D_TOPO
+    ):
+        assert np.isclose(value, reference, rtol=1e-6, atol=0), (
+            f"3D {label} differs serial vs np={uw.mpi.size}: "
+            f"{reference} vs {value}"
+        )
 
 
 def test_rotated_freeslip_box_nonlinear_partition_independent():
@@ -490,9 +526,12 @@ def test_rotated_freeslip_dynamic_topography_partition_independent():
         if uw.mpi.rank == 0:
             print(f"DIAG_SPHERICAL3D {_L2:.12e} its={_its} reason={_reason}")
     elif _kind == "spherical3d_topo":
-        _surface, _cmb = _spherical3d_topography_diagnostics()
+        _coefficients = _spherical3d_topography_diagnostics()
         if uw.mpi.rank == 0:
-            print(f"DIAG_SPHERICAL3D_TOPO {_surface:.12e} {_cmb:.12e}")
+            print(
+                "DIAG_SPHERICAL3D_TOPO "
+                + " ".join(f"{value:.12e}" for value in _coefficients)
+            )
     else:
         _L2, _verr = _box_diagnostics()
         if uw.mpi.rank == 0:
diff --git a/tests/parallel/test_1065_boundary_flux_parallel.py b/tests/parallel/test_1065_boundary_flux_parallel.py
index f3cb42a89..ce21e1cea 100644
--- a/tests/parallel/test_1065_boundary_flux_parallel.py
+++ b/tests/parallel/test_1065_boundary_flux_parallel.py
@@ -15,6 +15,7 @@
 import sympy
 import pytest
 import underworld3 as uw
+from mpi4py import MPI
 
 pytestmark = [pytest.mark.mpi(min_size=2), pytest.mark.timeout(180)]
 
@@ -73,6 +74,41 @@ def test_boundary_flux_partition_independent():
     assert relL2 < 0.01, f"heat flux relL2 vs analytic {relL2:.4f} too large at np={uw.mpi.size}"
 
 
+def _uniform_flux_3d_error(degree, mass):
+    """Maximum pointwise error for unit-cube conduction on the Bottom trace."""
+    mesh = uw.meshing.UnstructuredSimplexBox(
+        minCoords=(0.0, 0.0, 0.0),
+        maxCoords=(1.0, 1.0, 1.0),
+        cellSize=0.45,
+        regular=True,
+        qdegree=3,
+    )
+    temperature = uw.discretisation.MeshVariable(
+        f"Tbf3dp_p{degree}", mesh, 1, degree=degree
+    )
+    poisson = uw.systems.Poisson(mesh, u_Field=temperature)
+    poisson.constitutive_model = uw.constitutive_models.DiffusionModel
+    poisson.constitutive_model.Parameters.diffusivity = 1.0
+    poisson.f = 0.0
+    poisson.add_dirichlet_bc(0.0, "Bottom")
+    poisson.add_dirichlet_bc(1.0, "Top")
+    poisson.tolerance = 1.0e-11
+    poisson.petsc_options["snes_type"] = "ksponly"
+    poisson.solve()
+    _xs, flux = poisson.boundary_flux("Bottom", mass=mass)
+    local_error = float(np.max(np.abs(np.asarray(flux) + 1.0))) if len(flux) else 0.0
+    return uw.mpi.comm.allreduce(local_error, op=MPI.MAX)
+
+
+@pytest.mark.parametrize(("degree", "mass"), ((1, "lumped"), (2, "auto")))
+def test_boundary_flux_3d_pointwise_uniform_partition_independent(degree, mass):
+    """P1 and P2 constant-flux recovery is pointwise exact on every MPI partition."""
+    max_error = _uniform_flux_3d_error(degree, mass)
+    assert max_error < 1.0e-8, (
+        f"P{degree} pointwise flux error {max_error:.3e} at np={uw.mpi.size}"
+    )
+
+
 if __name__ == "__main__":
     _b, _r = _flux_diagnostics()
     if uw.mpi.rank == 0:
diff --git a/tests/test_1018_rotated_freeslip.py b/tests/test_1018_rotated_freeslip.py
index 27be50a53..2035d827f 100644
--- a/tests/test_1018_rotated_freeslip.py
+++ b/tests/test_1018_rotated_freeslip.py
@@ -156,7 +156,19 @@ def _spherical3d_reaction_topography():
     stokes.tolerance = 1.0e-5
     stokes.solve()
 
-    def harmonic_coefficient(boundary, response_sign):
+    dm = mesh.dm
+    csec = dm.getCoordinateSection()
+    cvec = np.asarray(dm.getCoordinatesLocal().array).reshape(-1, mesh.dim)
+    vertex_start, vertex_end = dm.getDepthStratum(0)
+    local_vertex_keys = {
+        tuple(np.round(cvec[csec.getOffset(point) // mesh.dim], 12))
+        for point in range(vertex_start, vertex_end)
+    }
+    vertex_keys = set()
+    for rank_keys in uw.mpi.comm.allgather(local_vertex_keys):
+        vertex_keys.update(rank_keys)
+
+    def harmonic_coefficients(boundary, response_sign):
         xs, sigma_nn = stokes.boundary_normal_traction(boundary)
         local = {
             tuple(np.round(x, 12)): -float(value)
@@ -169,24 +181,36 @@ def harmonic_coefficient(boundary, response_sign):
         topography = np.asarray(list(samples.values()))
         radii = np.linalg.norm(coords, axis=1)
         harmonic_values = 0.5 * (3.0 * (coords[:, 2] / radii) ** 2 - 1.0)
-        return float(
-            response_sign
-            * np.dot(topography, harmonic_values)
-            / np.dot(harmonic_values, harmonic_values)
-        )
+        is_vertex = np.array([key in vertex_keys for key in samples], dtype=bool)
 
-    surface = harmonic_coefficient("Upper", 1.0)
-    cmb = harmonic_coefficient("Lower", -1.0)
-    return surface, cmb
+        def fit(mask):
+            return float(
+                response_sign
+                * np.dot(topography[mask], harmonic_values[mask])
+                / np.dot(harmonic_values[mask], harmonic_values[mask])
+            )
+
+        return fit(np.ones(len(coords), dtype=bool)), fit(is_vertex), fit(~is_vertex)
+
+    return (
+        *harmonic_coefficients("Upper", 1.0),
+        *harmonic_coefficients("Lower", -1.0),
+    )
 
 
 @pytest.mark.level_2
 def test_rotated_freeslip_spherical3d_reaction_topography():
     """3D reaction loads must be divided by boundary mass to recover pointwise stress."""
 
-    surface, cmb = _spherical3d_reaction_topography()
+    surface, surface_vertices, surface_midpoints, cmb, cmb_vertices, cmb_midpoints = (
+        _spherical3d_reaction_topography()
+    )
     assert np.isclose(surface, 0.41920, rtol=0.10)
     assert np.isclose(cmb, 0.77060, rtol=0.10)
+    assert np.isclose(surface_vertices, 0.41920, rtol=0.12)
+    assert np.isclose(surface_midpoints, 0.41920, rtol=0.12)
+    assert np.isclose(cmb_vertices, 0.77060, rtol=0.12)
+    assert np.isclose(cmb_midpoints, 0.77060, rtol=0.12)
 
 
 def test_rotated_freeslip_annulus_zero_leakage():
diff --git a/tests/test_1019_boundary_flux.py b/tests/test_1019_boundary_flux.py
index 644a8af70..cf5f5e01b 100644
--- a/tests/test_1019_boundary_flux.py
+++ b/tests/test_1019_boundary_flux.py
@@ -59,6 +59,47 @@ def test_boundary_flux_scalar_heatflux_serial():
     assert abs(bd_q) > 0.0                                    # field populated + usable
 
 
+def _uniform_flux_3d(degree, mass):
+    """Unit-cube conduction with exact pointwise outward flux -1 on Bottom."""
+    mesh = uw.meshing.UnstructuredSimplexBox(
+        minCoords=(0.0, 0.0, 0.0),
+        maxCoords=(1.0, 1.0, 1.0),
+        cellSize=0.45,
+        regular=True,
+        qdegree=3,
+    )
+    temperature = uw.discretisation.MeshVariable(
+        f"Tbf3d_p{degree}", mesh, 1, degree=degree
+    )
+    poisson = uw.systems.Poisson(mesh, u_Field=temperature)
+    poisson.constitutive_model = uw.constitutive_models.DiffusionModel
+    poisson.constitutive_model.Parameters.diffusivity = 1.0
+    poisson.f = 0.0
+    poisson.add_dirichlet_bc(0.0, "Bottom")
+    poisson.add_dirichlet_bc(1.0, "Top")
+    poisson.tolerance = 1.0e-11
+    poisson.petsc_options["snes_type"] = "ksponly"
+    poisson.solve()
+    xs, flux = poisson.boundary_flux("Bottom", mass=mass)
+    return poisson, np.asarray(xs), np.asarray(flux)
+
+
+@pytest.mark.level_2
+@pytest.mark.parametrize(("degree", "mass"), ((1, "lumped"), (2, "auto")))
+def test_boundary_flux_3d_pointwise_uniform_serial(degree, mass):
+    """P1 and P2 recovery reproduce constant flux at every triangular-trace node."""
+    _poisson, _xs, flux = _uniform_flux_3d(degree, mass)
+    assert np.allclose(flux, -1.0, rtol=0.0, atol=1.0e-8)
+
+
+@pytest.mark.level_2
+def test_boundary_flux_3d_p2_lumped_rejected():
+    """P2 triangle vertex row sums are zero, so nodal lumping is not pointwise valid."""
+    poisson, _xs, _flux = _uniform_flux_3d(2, "consistent")
+    with pytest.raises(ValueError, match="zero row-sum mass"):
+        poisson.boundary_flux("Bottom", mass="lumped")
+
+
 if __name__ == "__main__":
     _f, _a, _b = _heatflux_diagnostics()
     c = np.dot(_f, _a) / (np.linalg.norm(_f) * np.linalg.norm(_a))

From 1533574a034d774f19d3a9aebb5ab7dd406995de Mon Sep 17 00:00:00 2001
From: Tyagi <gthyagi@gmail.com>
Date: Fri, 24 Jul 2026 08:58:19 +0530
Subject: [PATCH 3/3] Parameterize Zhong topography diagnostic resolution

Allow the serial and MPI Zhong l=2 rotated free-slip diagnostics to accept an optional shell cell size while retaining 0.25 as the regression-test default. Extend the standalone MPI diagnostic output with the selected cell size so manual convergence results are self-describing. This supports serial/MPI comparisons at 1/4, 1/8, and 1/16 resolution without adding costly high-resolution cases to routine CI.
---
 tests/parallel/test_1064_rotated_freeslip_parallel.py | 9 +++++----
 tests/test_1018_rotated_freeslip.py                   | 4 ++--
 2 files changed, 7 insertions(+), 6 deletions(-)

diff --git a/tests/parallel/test_1064_rotated_freeslip_parallel.py b/tests/parallel/test_1064_rotated_freeslip_parallel.py
index 4c882c389..576d5340a 100644
--- a/tests/parallel/test_1064_rotated_freeslip_parallel.py
+++ b/tests/parallel/test_1064_rotated_freeslip_parallel.py
@@ -169,12 +169,12 @@ def _spherical3d_diagnostics():
     return L2, int(info["ksp_its"]), int(info["ksp_reason"])
 
 
-def _spherical3d_topography_diagnostics():
+def _spherical3d_topography_diagnostics(cell_size=0.25):
     """Zhong l=2 topography from the 3D rotated-constraint reaction."""
     RI, RO, RINT = 0.55, 1.0, 0.775
     mesh = uw.meshing.SphericalShellInternalBoundary(
         radiusOuter=RO, radiusInternal=RINT, radiusInner=RI,
-        cellSize=0.25, qdegree=2, degree=1)
+        cellSize=cell_size, qdegree=2, degree=1)
     v = uw.discretisation.MeshVariable(
         "Vst", mesh, mesh.dim, degree=2, continuous=True)
     p = uw.discretisation.MeshVariable(
@@ -526,10 +526,11 @@ def test_rotated_freeslip_dynamic_topography_partition_independent():
         if uw.mpi.rank == 0:
             print(f"DIAG_SPHERICAL3D {_L2:.12e} its={_its} reason={_reason}")
     elif _kind == "spherical3d_topo":
-        _coefficients = _spherical3d_topography_diagnostics()
+        _cell_size = float(sys.argv[2]) if len(sys.argv) > 2 else 0.25
+        _coefficients = _spherical3d_topography_diagnostics(_cell_size)
         if uw.mpi.rank == 0:
             print(
-                "DIAG_SPHERICAL3D_TOPO "
+                f"DIAG_SPHERICAL3D_TOPO cell_size={_cell_size:.8f} "
                 + " ".join(f"{value:.12e}" for value in _coefficients)
             )
     else:
diff --git a/tests/test_1018_rotated_freeslip.py b/tests/test_1018_rotated_freeslip.py
index 2035d827f..f6db5294e 100644
--- a/tests/test_1018_rotated_freeslip.py
+++ b/tests/test_1018_rotated_freeslip.py
@@ -120,7 +120,7 @@ def test_rotated_freeslip_spherical_shell_3d():
         assert rotfrac < 1e-8, f"rotation mode {k} gauge {rotfrac:.2e} not removed"
 
 
-def _spherical3d_reaction_topography():
+def _spherical3d_reaction_topography(cell_size=0.25):
     """Zhong l=2 topography recovered directly from rotated constraint reactions."""
 
     radius_inner = 0.55
@@ -130,7 +130,7 @@ def _spherical3d_reaction_topography():
         radiusOuter=radius_outer,
         radiusInternal=radius_internal,
         radiusInner=radius_inner,
-        cellSize=0.25,
+        cellSize=cell_size,
         qdegree=2,
         degree=1,
     )