diff --git a/src/underworld3/cython/petsc_maths.pyx b/src/underworld3/cython/petsc_maths.pyx
index 4061bf298..effaa9511 100644
--- a/src/underworld3/cython/petsc_maths.pyx
+++ b/src/underworld3/cython/petsc_maths.pyx
@@ -259,14 +259,14 @@ class CellWiseIntegral:
         super().__init__()
 
     @timing.routine_timer_decorator
-    def evaluate(self) -> float:
+    def evaluate(self):
         """
         Evaluate the cell-wise integral and return results per cell.
 
         Returns
         -------
         ndarray
-            Array of integral values, one per mesh cell.
+            Array of integral values, one per (rank-local) mesh cell.
 
         Raises
         ------
@@ -300,27 +300,39 @@ class CellWiseIntegral:
         cdef Vec cgvec
         cgvec = a_global
 
-        # TODO(BUG): clone+createDefault+createDS gives dmc a single P1 field,
-        # but cgvec is packed for mesh.dm's multi-field layout — the integral
-        # reads the wrong DOFs and over-counts by ~2x on the unit square.
-        # Same bug as PR #172 (reverted in PR #173-followup). Tests in
-        # tests/test_0501_integrals.py::test_cellwise_integrate_* are xfail
-        # until this is rewritten to integrate against mesh.dm + getDS()
-        # directly (the pre-PR-172 Integral pattern).
-        cdef DM dmc = self.mesh.dm.clone()
-        cdef FE fec = FE().createDefault(self.mesh.dim, 1, False, -1)
-        dmc.setField(0, fec)
-        dmc.createDS()
-
-        cdef DS ds = dmc.getDS()
+        # Integrate on the mesh DM with its own DS (the `Integral` pattern
+        # above).  The JIT objective is compiled against the mesh's
+        # multi-field packing, so it must be evaluated with the mesh DS and
+        # the mesh-packed solution vector.  An earlier clone + createDefault
+        # + createDS variant attached a single P1 field to the cloned DM
+        # while cgvec remained packed for mesh.dm — the layout mismatch made
+        # the integral read the wrong DOFs and over-count (~2x for fn=1 on
+        # the unit square; the same defect that forced the PR #172 revert
+        # for Integral, see 88807c26).
+        #
+        # Note: setting the objective on the shared mesh DS re-inherits the
+        # issue-#171 repeated-call cost-growth behaviour that Integral has —
+        # accepted here for correctness parity.
+        cdef DM dm = self.mesh.dm
+        cdef DS ds = self.mesh.dm.getDS()
         CHKERRQ( PetscDSSetObjective(ds.ds, 0, ext.fns_residual[0]) )
 
-        cdef Vec rvec = dmc.createGlobalVec()
-        CHKERRQ( DMPlexComputeCellwiseIntegralFEM(dmc.dm, cgvec.vec, rvec.vec, NULL) )
+        # DMPlexComputeCellwiseIntegralFEM writes Nf scalars per cell into a
+        # flat [cell*Nf + field] layout when the output vector carries no
+        # DM/section.  Only field 0 has an objective; the other field slots
+        # are left at zero.
+        num_fields = self.mesh.dm.getNumFields()
+        cStart, cEnd = self.mesh.dm.getHeightStratum(0)
+        num_cells = cEnd - cStart
+
+        rvec_py = PETSc.Vec().createMPI((num_cells * num_fields, None),
+                                        comm=self.mesh.dm.comm)
+        cdef Vec rvec = rvec_py
+        CHKERRQ( DMPlexComputeCellwiseIntegralFEM(dm.dm, cgvec.vec, rvec.vec, NULL) )
         self.mesh.dm.restoreGlobalVec(a_global)
 
-        results = rvec.array.copy()
-        rvec.destroy()
+        results = rvec_py.array.reshape(-1, num_fields)[:, 0].copy()
+        rvec_py.destroy()
 
         return results
 
diff --git a/tests/test_0501_integrals.py b/tests/test_0501_integrals.py
index 309e22cae..4c27108b8 100644
--- a/tests/test_0501_integrals.py
+++ b/tests/test_0501_integrals.py
@@ -169,31 +169,26 @@ def test_integrate_swarmvar_deriv_00():
 
 
 # CellWiseIntegral parity tests — the sister class of Integral that returns a
-# per-cell array.  These are marked xfail because CellWiseIntegral.evaluate()
-# clones mesh.dm and attaches a fresh single-field discretisation
-# (FE.createDefault + setField + createDS), but then integrates against
-# mesh.dm.getGlobalVec() which is packed for the original multi-field layout.
-# The mismatch produces a ~2× over-count on simple cases — the same bug that
-# motivated the revert of PR #172 for Integral.
-#
-# These tests will start passing once CellWiseIntegral.evaluate() is rewritten
-# to integrate against mesh.dm + mesh.dm.getDS() directly (matching the
-# pre-PR-172 Integral pattern).
-@pytest.mark.xfail(reason="CellWiseIntegral has the same field-cloning bug "
-                          "that motivated the PR #172 revert")
+# per-cell array.  Regression tests for the field-cloning bug (LE-02 / BF-03):
+# evaluate() used to clone mesh.dm with a fresh single-P1-field DS while the
+# solution vector remained packed for mesh.dm's multi-field layout, reading
+# the wrong DOFs and over-counting (~2x for fn=1 on the unit square — the
+# same defect that motivated the revert of PR #172 for Integral).
 def test_cellwise_integrate_constants():
 
     calculator = uw.maths.CellWiseIntegral(mesh, fn=1.0)
     cell_values = calculator.evaluate()
 
+    # evaluate() returns one value per rank-local cell; reduce for the
+    # global total (no-op in serial).
+    total = uw.mpi.comm.allreduce(np.sum(cell_values))
+
     # Sum of per-cell volumes equals the total mesh volume (unit square = 1.0).
-    assert abs(np.sum(cell_values) - 1.0) < 0.001
+    assert abs(total - 1.0) < 0.001
 
     return
 
 
-@pytest.mark.xfail(reason="CellWiseIntegral has the same field-cloning bug "
-                          "that motivated the PR #172 revert")
 def test_cellwise_integrate_meshvar():
 
     s_soln.data[:, 0] = np.sin(np.pi * s_soln.coords[:, 0])
@@ -201,7 +196,11 @@ def test_cellwise_integrate_meshvar():
     calculator = uw.maths.CellWiseIntegral(mesh, fn=s_soln.sym[0])
     cell_values = calculator.evaluate()
 
+    # evaluate() returns one value per rank-local cell; reduce for the
+    # global total (no-op in serial).
+    total = uw.mpi.comm.allreduce(np.sum(cell_values))
+
     # Sum matches the global Integral: ∫∫ sin(πx) dx dy = 2/π over the unit square.
-    assert abs(np.sum(cell_values) - 2 / np.pi) < 0.001
+    assert abs(total - 2 / np.pi) < 0.001
 
     return