From 443ffb58fb8acfab2b2a0737225626dc6f985c4a Mon Sep 17 00:00:00 2001 From: Tyagi Date: Sun, 6 Sep 2026 01:17:47 +1000 Subject: [PATCH 1/3] fix: make stabilization cell sizes local to each cell (#687) Adapt only the mesh-size correction from lmoresi's 68e545fd on feature/navier-stokes-supg; do not import Navier-Stokes or other branch changes. Cache _radii_own from current DM vertex coordinates and use it for mesh.cell_size(). Preserve the legacy kd-tree radius arrays and global timestep/mesh-motion consumers. Use coordinate-section offsets and the full vertex stratum so the own-cell RMS definition also handles hexahedra, which have eight vertices but six faces. Correct the field documentation and Nitsche mechanism tests for the new definition; retain physical solve tolerances and use the exact nearest-centroid <= own-centroid ordering instead of an arbitrary approximate-equality tolerance. Add a first-failing independent geometry/deformation regression for triangles, tetrahedra, quadrilaterals and hexahedra plus a regular-square analytical control. Before: four failures in serial and on eight ranks. After rebuild: 21 passed/one expected skip serial (22.90 s), 22 passed on eight ranks (40.45 s), covering Nitsche solves, radius accessors, frozen PC2 migration and memory/disk snapshots. Own-cell geometry error is zero in these tests; style and whitespace gates pass. --- .../discretisation/discretisation_mesh.py | 48 +++++++------- tests/test_0010_cell_size_geometry.py | 66 +++++++++++++++++++ tests/test_1065_nitsche_local_h.py | 23 +++---- 3 files changed, 103 insertions(+), 34 deletions(-) create mode 100644 tests/test_0010_cell_size_geometry.py diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index d22d8b190..e48d3a56f 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -3224,7 +3224,8 @@ def cell_size(self): Returns the ``.sym`` of a cell-constant (degree-0, discontinuous) scalar MeshVariable holding each cell's characteristic length (the - ``volume**(1/dim)`` equivalent radius, i.e. ``self._radii``). Unlike + RMS distance of its vertices from their own centroid). This is a + purely cell-local quantity, independent of the MPI partition. Unlike the single *global* scalar from :meth:`get_min_radius` (the smallest cell anywhere), this varies cell to cell, so a stabilisation that scales as :math:`1/h` — e.g. the Nitsche free-slip penalty @@ -3290,9 +3291,9 @@ def _refresh(): def _assemble_cell_size(self, var): """Fill ``var`` (degree-0 scalar) with each cell's characteristic size. - Uses the per-cell characteristic lengths ``self._radii`` computed by + Uses the own-cell characteristic lengths ``self._radii_own`` computed by :meth:`_get_mesh_sizes` on the *current* geometry. A degree-0 - discontinuous variable's local DOFs and ``self._radii`` are BOTH + discontinuous variable's local DOFs and ``self._radii_own`` are BOTH indexed by this rank's cell-stratum order, so a direct assignment is correct on every rank. @@ -3300,24 +3301,9 @@ def _assemble_cell_size(self, var): access, no collective): mixing a rank-local fast path with a collective fallback would diverge across ranks and deadlock, because ``var.coords`` triggers the collective ``_get_coords_for_basis``.""" - # TODO(BUG): this field is PARTITION-DEPENDENT, and so therefore is the - # Nitsche penalty gamma*mu/h that consumes it (local_h=True, the default). - # Not the indexing here — the values. `_get_mesh_sizes` measures a cell by - # the distance from its vertices to the NEAREST CENTROID in a kd-tree built - # from THIS RANK's centroids, so near a partition seam the nearest centroid - # may simply be absent. Measured on Annulus(cellSize=0.12): the field's sum - # is 26.0822 at np=1, 26.1211 at np=2 and 26.1386 at np=4, and its max moves - # at np=4. End to end that is 6.6e-03 in the velocity of a Nitsche free-slip - # annulus and it does NOT shrink with solver tolerance. - # This is a DIFFERENT defect from the boundary normal fixed for #564 (which - # is now clean: the same solve with local_h=False agrees to 3.6e-10 at - # np=1..4). It is the local h that is left, and it also reaches every other - # consumer of `cell_size()`. Not fixed here because `_get_mesh_sizes` also - # feeds `get_min_radius`, the adaptivity metrics and the free-surface - # relaxation, and it needs its own benchmarking. - # Guard/measurement: tests/parallel/test_1069_boundary_normal_parallel.py - # (_nitsche_annulus_diagnostics docstring records the numbers). - radii = numpy.asarray(self._radii).reshape(-1) + # Own-cell radii fix #687 without changing the legacy kd-tree radii + # used by global timestep estimates, adaptivity, and mesh relaxation. + radii = numpy.asarray(self._radii_own).reshape(-1) # Empty partition (no local cells): nothing to fill on this rank. if radii.size == 0 or var.data.shape[0] == 0: return @@ -6895,8 +6881,11 @@ def _eval_use_robust_location(self) -> bool: def _get_mesh_sizes(self, verbose=False): """ - Obtain the (local) mesh radii and centroids using kdtree distances - This routine is called when the mesh is built / rebuilt + Cache own-cell radii for cell_size and return legacy kd-tree radii. + + Own-cell sizes use current DM vertices, so neither partition-local + neighbours nor stale coordinate views affect stabilization (#687). + Legacy radii remain unchanged for their other consumers. """ centroids = self._get_coords_for_basis(0, False) @@ -6909,6 +6898,9 @@ def _get_mesh_sizes(self, verbose=False): cell_length = np.empty(centroids.shape[0]) cell_min_r = np.empty(centroids.shape[0]) cell_r = np.empty(centroids.shape[0]) + cell_r_own = np.empty(centroids.shape[0]) + coordinate_section = self.dm.getCoordinateDM().getLocalSection() + vertex_coordinates = self.dm.getCoordinatesLocal().array for cell in range(cEnd - cStart): cell_num_points = self.dm.getConeSize(cell) @@ -6922,6 +6914,16 @@ def _get_mesh_sizes(self, verbose=False): cell_r[cell] = np.sqrt(distsq.mean()) cell_min_r[cell] = np.sqrt(distsq.min()) + # A hex has six faces but eight vertices: select the vertex + # stratum, not a cone-sized suffix of its transitive closure. + closure = self.dm.getTransitiveClosure(cStart + cell)[0] + vertices = closure[(closure >= pStart) & (closure < pEnd)] + offsets = np.array([coordinate_section.getOffset(int(v)) for v in vertices]) + own_coords = vertex_coordinates[offsets[:, None] + np.arange(self.cdim)] + delta = own_coords - own_coords.mean(axis=0) + cell_r_own[cell] = np.sqrt(np.mean(np.sum(delta ** 2, axis=1))) + + self._radii_own = cell_r_own return cell_min_r, cell_r, centroids, cell_length # ========== diff --git a/tests/test_0010_cell_size_geometry.py b/tests/test_0010_cell_size_geometry.py new file mode 100644 index 000000000..9e955328e --- /dev/null +++ b/tests/test_0010_cell_size_geometry.py @@ -0,0 +1,66 @@ +"""Issue #687: cell_size is an own-cell geometric quantity, including after deform. + +The independent oracle reads vertex coordinates through the coordinate section; +it does not use the mesh's cached radii or centroid kd-tree. Run serial and MPI. +""" + +import numpy as np +import pytest + +import underworld3 as uw + +pytestmark = [pytest.mark.level_1, pytest.mark.tier_b] + + +def _vertex_rms(mesh): + dm = mesh.dm + section = dm.getCoordinateDM().getLocalSection() + coordinates = dm.getCoordinatesLocal().array + start, end = dm.getHeightStratum(0) + first_vertex, last_vertex = dm.getDepthStratum(0) + radii = [] + for cell in range(start, end): + vertices = [int(point) for point in dm.getTransitiveClosure(cell)[0] + if first_vertex <= point < last_vertex] + points = np.array([coordinates[section.getOffset(v):section.getOffset(v) + mesh.cdim] + for v in vertices]) + radii.append(np.sqrt(np.mean(np.sum((points - points.mean(axis=0)) ** 2, axis=1)))) + return np.asarray(radii) + + +@pytest.mark.parametrize("dim", [2, 3]) +@pytest.mark.parametrize("simplex", [True, False], ids=["simplex", "tensor"]) +def test_cell_size_matches_own_vertices_and_tracks_deform(dim, simplex): + geometry = dict(minCoords=(0.0,) * dim, maxCoords=(1.0,) * dim, qdegree=3) + mesh = (uw.meshing.UnstructuredSimplexBox(**geometry, cellSize=0.25, regular=False) + if simplex else uw.meshing.StructuredQuadBox(**geometry, elementRes=(4,) * dim)) + mesh.cell_size() + field = mesh._cell_size_variable + errors = [] + for phase in ("initial", "deformed"): + if phase == "deformed": + coordinates = np.array(mesh.X.coords) + coordinates[:, 0] = 1.7 * coordinates[:, 0] + 0.2 * coordinates[:, 1] + mesh.deform(coordinates) + expected = _vertex_rms(mesh) + actual = np.asarray(field.array[:, 0, 0]) + shapes_match = actual.shape == expected.shape + assert all(uw.mpi.comm.allgather(shapes_match)), (actual.shape, expected.shape) + local_error = float(np.abs(actual - expected).max(initial=0.0)) + error = max(uw.mpi.comm.allgather(local_error)) + errors.append(error) + uw.pprint(f"CELL_SIZE_GEOMETRY dim={dim} simplex={simplex} phase={phase} " + f"ranks={uw.mpi.size} max_error={error:.12g}") + assert max(errors) < 1e-12, errors + + +def test_regular_square_cell_size_keeps_global_radius(): + mesh = uw.meshing.StructuredQuadBox(elementRes=(4, 4), qdegree=2) + legacy = np.array(mesh._radii) + global_radius = mesh.get_min_radius() + mesh.cell_size() + expected = np.sqrt(2.0) / 8.0 + error = float(np.abs(np.asarray(mesh._cell_size_variable.array) - expected).max(initial=0.0)) + assert max(uw.mpi.comm.allgather(error)) < 1e-12 + assert global_radius == pytest.approx(expected, rel=1e-12) + assert all(uw.mpi.comm.allgather(np.array_equal(mesh._radii, legacy))) diff --git a/tests/test_1065_nitsche_local_h.py b/tests/test_1065_nitsche_local_h.py index 15c0bcb3b..971130a0f 100644 --- a/tests/test_1065_nitsche_local_h.py +++ b/tests/test_1065_nitsche_local_h.py @@ -152,11 +152,11 @@ def _box_wobble(X0, amp): # -------------------------------------------------------------------------- def test_cell_size_is_local_per_cell(): """``mesh.cell_size()`` is a per-cell field equal to each cell's - characteristic size (``mesh._radii``), not the single global minimum.""" + own-vertex RMS size (``mesh._radii_own``), not the single global minimum.""" mesh = _graded_box() h = mesh.cell_size() # sympy symbol -> backed by a P0 field - field = np.asarray(mesh._cell_size_variable.data[:, 0]).reshape(-1) - radii = np.asarray(mesh._radii).reshape(-1) + field = np.asarray(mesh._cell_size_variable.array[:, 0, 0]).reshape(-1) + radii = np.asarray(mesh._radii_own).reshape(-1) # field exactly mirrors the per-cell characteristic size (rank-local check, # reduced to a single global pass/fail so all ranks agree) @@ -168,8 +168,9 @@ def test_cell_size_is_local_per_cell(): gfmin, gfmax = _gmin(field), _gmax(field) assert gfmax / gfmin > 3.0 - # the global scalar that global-h would use is just the minimum cell size - assert np.isclose(mesh.get_min_radius(), gfmin, rtol=1e-6) + # The unchanged nearest-centroid minimum cannot exceed the own-cell + # minimum; these are no longer the same definition on an irregular mesh. + assert 0.0 < mesh.get_min_radius() <= gfmin * (1.0 + 1e-12) def test_local_h_at_coarse_freeslip_boundary_exceeds_global_min(): @@ -177,13 +178,13 @@ def test_local_h_at_coarse_freeslip_boundary_exceeds_global_min(): times the global minimum. global-h would over-stiffen that penalty by exactly this factor; local-h scales it correctly.""" mesh = _graded_box(h_fine=0.04, h_coarse=0.12) - # build/exercise the field; its data equals mesh._radii (asserted in + # build/exercise the field; its data equals mesh._radii_own (asserted in # test_cell_size_is_local_per_cell), so we read the per-cell sizes directly - # from _radii / _centroids — a rank-local lookup, avoiding the collective + # from the field / _centroids — a rank-local lookup, avoiding the collective # arbitrary-point uw.function.evaluate (which deadlocks in parallel). _ = mesh.cell_size() cen = np.asarray(mesh._centroids) - radii = np.asarray(mesh._radii).reshape(-1) + radii = np.asarray(mesh._cell_size_variable.array[:, 0, 0]).reshape(-1) near_top = cen[:, 1] > 0.85 # cells adjacent to the Top free-slip edge h_top = radii[near_top] @@ -204,14 +205,14 @@ def test_cell_size_tracks_deformation(): the Nitsche mis-scaling on the free surface.""" mesh = _graded_box() _ = mesh.cell_size() - h_before = mesh._cell_size_variable.data[:, 0].copy() + h_before = np.array(mesh._cell_size_variable.array[:, 0, 0]) X = np.asarray(mesh.X.coords).copy() moved = mesh.deform(_box_wobble(X, amp=0.04)) assert moved # geometry actually changed - h_after = mesh._cell_size_variable.data[:, 0].copy() - radii_after = np.asarray(mesh._radii).reshape(-1) + h_after = np.array(mesh._cell_size_variable.array[:, 0, 0]) + radii_after = np.asarray(mesh._radii_own).reshape(-1) # not stale: the field changed with the geometry SOMEWHERE (global OR) ... nb = min(h_after.shape[0], h_before.shape[0]) From 6fde1ac931c39632a930418a2f461c6ae112a4b5 Mon Sep 17 00:00:00 2001 From: Tyagi Date: Sun, 6 Sep 2026 21:37:44 +1000 Subject: [PATCH 2/3] test: prove cell-size partition independence Rename the new per-cell geometric radius cache from _radii_own to _cell_radii so the name describes cell geometry rather than rank ownership. Update the focused Nitsche and deformation checks accordingly.\n\nAdd an enumerated parallel regression that gathers owned-cell centroid/radius pairs and compares the complete sorted table with a fresh single-rank run on the same cached Gmsh mesh. This directly guards the rank-count-independence claim at np=2, np=4 and np=8 instead of relying only on within-rank geometric identities.\n\nValidated locally with 9 focused serial tests and the new MPI test at 2, 4 and 8 ranks. --- .../discretisation/discretisation_mesh.py | 12 ++--- ...t_1077_cell_size_partition_independence.py | 47 +++++++++++++++++++ tests/test_1065_nitsche_local_h.py | 8 ++-- 3 files changed, 57 insertions(+), 10 deletions(-) create mode 100644 tests/parallel/test_1077_cell_size_partition_independence.py diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index e48d3a56f..51f260e61 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -3291,9 +3291,9 @@ def _refresh(): def _assemble_cell_size(self, var): """Fill ``var`` (degree-0 scalar) with each cell's characteristic size. - Uses the own-cell characteristic lengths ``self._radii_own`` computed by + Uses the cell-geometry characteristic lengths ``self._cell_radii`` computed by :meth:`_get_mesh_sizes` on the *current* geometry. A degree-0 - discontinuous variable's local DOFs and ``self._radii_own`` are BOTH + discontinuous variable's local DOFs and ``self._cell_radii`` are BOTH indexed by this rank's cell-stratum order, so a direct assignment is correct on every rank. @@ -3303,7 +3303,7 @@ def _assemble_cell_size(self, var): ``var.coords`` triggers the collective ``_get_coords_for_basis``.""" # Own-cell radii fix #687 without changing the legacy kd-tree radii # used by global timestep estimates, adaptivity, and mesh relaxation. - radii = numpy.asarray(self._radii_own).reshape(-1) + radii = numpy.asarray(self._cell_radii).reshape(-1) # Empty partition (no local cells): nothing to fill on this rank. if radii.size == 0 or var.data.shape[0] == 0: return @@ -6898,7 +6898,7 @@ def _get_mesh_sizes(self, verbose=False): cell_length = np.empty(centroids.shape[0]) cell_min_r = np.empty(centroids.shape[0]) cell_r = np.empty(centroids.shape[0]) - cell_r_own = np.empty(centroids.shape[0]) + cell_radii = np.empty(centroids.shape[0]) coordinate_section = self.dm.getCoordinateDM().getLocalSection() vertex_coordinates = self.dm.getCoordinatesLocal().array @@ -6921,9 +6921,9 @@ def _get_mesh_sizes(self, verbose=False): offsets = np.array([coordinate_section.getOffset(int(v)) for v in vertices]) own_coords = vertex_coordinates[offsets[:, None] + np.arange(self.cdim)] delta = own_coords - own_coords.mean(axis=0) - cell_r_own[cell] = np.sqrt(np.mean(np.sum(delta ** 2, axis=1))) + cell_radii[cell] = np.sqrt(np.mean(np.sum(delta ** 2, axis=1))) - self._radii_own = cell_r_own + self._cell_radii = cell_radii return cell_min_r, cell_r, centroids, cell_length # ========== diff --git a/tests/parallel/test_1077_cell_size_partition_independence.py b/tests/parallel/test_1077_cell_size_partition_independence.py new file mode 100644 index 000000000..4348ece63 --- /dev/null +++ b/tests/parallel/test_1077_cell_size_partition_independence.py @@ -0,0 +1,47 @@ +"""Rank-count regression for the cell-local stabilization length. + +The parallel result is compared with a fresh single-rank run on the same Gmsh +mesh. This checks the complete cell geometry table, not merely a reduction or +a within-rank geometric identity. +""" + +import numpy as np +import pytest + +import underworld3 as uw + +from serial_reference import emit, mesh_fingerprint, serial_reference + + +pytestmark = [pytest.mark.mpi(min_size=2), pytest.mark.timeout(300)] + + +def _cell_geometry_table(): + mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.12, qdegree=2) + mesh.cell_size() + + local = np.column_stack((mesh._centroids, mesh._cell_radii)) + local = local[mesh._get_owned_cells_mask()] + gathered = uw.mpi.comm.allgather(local) + table = np.vstack(gathered) + order = np.lexsort(tuple(table[:, axis] for axis in reversed(range(mesh.dim)))) + return table[order].reshape(-1), mesh_fingerprint(mesh) + + +def test_cell_size_matches_single_rank_cell_by_cell(): + values, fingerprint = _cell_geometry_table() + reference = serial_reference(__file__, "cell_size") + + assert int(fingerprint[0]) == int(reference["fingerprint"][0]) + assert np.isclose( + fingerprint[1], reference["fingerprint"][1], rtol=1.0e-12, atol=0.0 + ) + + expected = np.asarray(reference["values"]) + assert values.shape == expected.shape + np.testing.assert_allclose(values, expected, rtol=0.0, atol=1.0e-14) + + +if __name__ == "__main__": + _values, _fingerprint = _cell_geometry_table() + emit(_values, _fingerprint) diff --git a/tests/test_1065_nitsche_local_h.py b/tests/test_1065_nitsche_local_h.py index 971130a0f..e6c225611 100644 --- a/tests/test_1065_nitsche_local_h.py +++ b/tests/test_1065_nitsche_local_h.py @@ -152,11 +152,11 @@ def _box_wobble(X0, amp): # -------------------------------------------------------------------------- def test_cell_size_is_local_per_cell(): """``mesh.cell_size()`` is a per-cell field equal to each cell's - own-vertex RMS size (``mesh._radii_own``), not the single global minimum.""" + own-vertex RMS size (``mesh._cell_radii``), not the single global minimum.""" mesh = _graded_box() h = mesh.cell_size() # sympy symbol -> backed by a P0 field field = np.asarray(mesh._cell_size_variable.array[:, 0, 0]).reshape(-1) - radii = np.asarray(mesh._radii_own).reshape(-1) + radii = np.asarray(mesh._cell_radii).reshape(-1) # field exactly mirrors the per-cell characteristic size (rank-local check, # reduced to a single global pass/fail so all ranks agree) @@ -178,7 +178,7 @@ def test_local_h_at_coarse_freeslip_boundary_exceeds_global_min(): times the global minimum. global-h would over-stiffen that penalty by exactly this factor; local-h scales it correctly.""" mesh = _graded_box(h_fine=0.04, h_coarse=0.12) - # build/exercise the field; its data equals mesh._radii_own (asserted in + # build/exercise the field; its data equals mesh._cell_radii (asserted in # test_cell_size_is_local_per_cell), so we read the per-cell sizes directly # from the field / _centroids — a rank-local lookup, avoiding the collective # arbitrary-point uw.function.evaluate (which deadlocks in parallel). @@ -212,7 +212,7 @@ def test_cell_size_tracks_deformation(): assert moved # geometry actually changed h_after = np.array(mesh._cell_size_variable.array[:, 0, 0]) - radii_after = np.asarray(mesh._radii_own).reshape(-1) + radii_after = np.asarray(mesh._cell_radii).reshape(-1) # not stale: the field changed with the geometry SOMEWHERE (global OR) ... nb = min(h_after.shape[0], h_before.shape[0]) From e6eaac28c5101e14545b4dcd78de5853ca2677de Mon Sep 17 00:00:00 2001 From: Tyagi Date: Sun, 6 Sep 2026 21:56:35 +1000 Subject: [PATCH 3/3] test: restore default Nitsche local-h path Remove the local_h=False workaround from the boundary-normal MPI regression now that Mesh.cell_size() is partition independent. The test again exercises the public local_h=True default and compares its Nitsche solve with a fresh serial process.\n\nRecord the user-visible consequence in the development changelog: the rank-local centroid kd-tree moved the default Nitsche velocity answer by 6.6e-3, while the cell-geometry replacement is identical cell by cell from one through eight ranks.\n\nValidated the focused Nitsche regression at 2, 4 and 8 Open MPI ranks (10.99 s, 7.31 s and 9.60 s respectively). --- docs/developer/CHANGELOG.md | 5 +++++ .../test_1069_boundary_normal_parallel.py | 16 +++++++--------- 2 files changed, 12 insertions(+), 9 deletions(-) diff --git a/docs/developer/CHANGELOG.md b/docs/developer/CHANGELOG.md index d1af2c6f8..2e5d432c9 100644 --- a/docs/developer/CHANGELOG.md +++ b/docs/developer/CHANGELOG.md @@ -506,6 +506,11 @@ in `Stokes_Constrained` (#224), then made parallel-correct. minimum radius, restoring correct stiffness on graded and adapted meshes (#275). +- The local size now comes from each cell's own geometry instead of a kd-tree + over the centroids held by the current MPI rank. The old field changed at + partition boundaries and moved the default ``local_h=True`` Nitsche velocity + answer by 6.6e-3 between rank counts; the replacement is cell-by-cell + identical from one to eight ranks (#569, #687). - `mesh.boundary_slip` API with `BoundingSurface` objects for boundary tangent-slip (#225); `Surface.influence_function` respects finite edges (#241). diff --git a/tests/parallel/test_1069_boundary_normal_parallel.py b/tests/parallel/test_1069_boundary_normal_parallel.py index 2eb43c712..036f39a45 100644 --- a/tests/parallel/test_1069_boundary_normal_parallel.py +++ b/tests/parallel/test_1069_boundary_normal_parallel.py @@ -333,14 +333,12 @@ def _nitsche_annulus_diagnostics(): leakage. Both are stable from tolerance 1e-9 to 1e-12, so neither is the linear solve. - ``local_h=False`` is deliberate and it is not a workaround for this fix. The - default ``local_h=True`` scales the Nitsche penalty by ``mesh.cell_size()``, which - is built from ``Mesh._get_mesh_sizes`` — a kd-tree query against THIS RANK's cell - centroids, and so partition-dependent in its own right (on this mesh the field's - sum is 26.0822 at np=1, 26.1211 at np=2, 26.1386 at np=4, and its max moves at - np=4). That is a SEPARATE defect from the boundary normal, it is not what #564 is - about, and leaving it in would make this test measure the two together. See the - TODO(BUG) on ``Mesh._assemble_cell_size``. + This test now leaves ``local_h`` at its default ``True``. Before #569/#687, + doing so mixed the boundary-normal regression with a second partition-dependent + input from ``mesh.cell_size()``; this test therefore had to disable the public + default. The cell-local geometric size is now partition independent, so retaining + the default jointly guards the normal assembly and the Nitsche penalty path users + actually run. """ RI, RO = 0.5, 1.0 mesh = uw.meshing.Annulus(radiusInner=RI, radiusOuter=RO, cellSize=0.12, qdegree=3) @@ -357,7 +355,7 @@ def _nitsche_annulus_diagnostics(): y / r * sympy.cos(4 * theta) * (r - RI) * (RO - r) * 40.0]]) stokes.add_essential_bc((0.0, 0.0), "Lower") # default normal= is the assembled one — that is what is under test - stokes.add_nitsche_bc(0.0, "Upper", local_h=False) + stokes.add_nitsche_bc(0.0, "Upper") stokes.tolerance = 1.0e-9 stokes.petsc_options["snes_type"] = "ksponly" stokes.solve()