Rotated and block constraints are exclusive; remove four dead label lookups

src/underworld3/cython/petsc_generic_snes_solvers.pyx

@@ -3545,8 +3545,6 @@ class SNES_Scalar(SolverBaseClass):
             value = mesh.boundaries[bc.boundary].value
             ind = value
 
-            bc_label = self.dm.getLabel(boundary)
-            bc_is = bc_label.getStratumIS(value)
             self.natural_bcs[index] = self.natural_bcs[index]._replace(boundary_label_val=value)
 
             # use type 5 bc for `DM_BC_ESSENTIAL_FIELD` enum
@@ -3685,9 +3683,6 @@ class SNES_Scalar(SolverBaseClass):
             boundary = bc.boundary
             value = mesh.boundaries[bc.boundary].value
 
-            bc_label = mesh.dm.getLabel(boundary)
-            bc_is = bc_label.getStratumIS(value)
-
             if bc.fn_f is not None:
 
                 bd_F0  = sympy.Array(bc.fn_f)
@@ -4490,8 +4485,6 @@ class SNES_Vector(SolverBaseClass):
             value = mesh.boundaries[bc.boundary].value
             ind = value
 
-            bc_label = self.dm.getLabel(boundary)
-            bc_is = bc_label.getStratumIS(value)
             self.natural_bcs[index] = self.natural_bcs[index]._replace(boundary_label_val=value)
 
             # use type 5 bc for `DM_BC_ESSENTIAL_FIELD` enum
@@ -5297,8 +5290,6 @@ class SNES_MultiComponent(SolverBaseClass):
             value = mesh.boundaries[bc.boundary].value
             ind = value
 
-            bc_label = self.dm.getLabel(boundary)
-            bc_is = bc_label.getStratumIS(value)
             self.natural_bcs[index] = self.natural_bcs[index]._replace(boundary_label_val=value)
 
             bc_type = 6
@@ -6108,6 +6099,69 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
         # this attrib records if we need to re-setup
         self.is_setup = False
 
+    def _reject_mixed_constraint_mechanisms(self, adding):
+        """Refuse a rotated constraint and a block constraint on one solver (#464).
+
+        The rotated driver builds its own index-set fieldsplit over exactly two
+        fields (``rotated_bc._solve_rotated_iterative``), and ``build_rotation``
+        addresses velocity and pressure by field number. A block constraint
+        registers a multiplier field of its own, and those DOFs are in neither
+        index set — the preconditioner would then be built over a strict subset
+        of the operator's rows, with nothing said about it.
+
+        The two are alternative ways to impose the same wall-normal condition,
+        so asking for both is a configuration error rather than a case to
+        support. Supporting it would need a third split for the multipliers and
+        a ``build_rotation`` that knows about them.
+
+        Parameters
+        ----------
+        adding : str
+            Name of the method being called, so the message can say which
+            mechanism is already in place and which one was refused.
+        """
+
+        rotated = list(getattr(self, "_rotated_freeslip_bcs", None) or []) + list(
+            getattr(self, "_fault_contact_faults", None) or []
+        )
+        multipliers = list(getattr(self, "_multipliers", None) or [])
+
+        if adding == "solve":
+            # The dispatch reads both lists, so it can only report the pair.
+            if not (rotated and multipliers):
+                return
+            raise RuntimeError(
+                f"solve(): this solver carries {len(rotated)} rotated "
+                f"(free-slip or fault contact) and {len(multipliers)} "
+                f"block-constraint boundary condition(s). The rotated solve "
+                f"splits velocity and pressure by field number and the "
+                f"multiplier fields lie outside that split, so the "
+                f"preconditioner would cover only part of the operator. Both "
+                f"impose the same wall-normal condition — use one of them "
+                f"(issue #464)."
+            )
+
+        if adding == "add_constraint_bc":
+            if not rotated:
+                return
+            present, refused = len(rotated), "a block constraint"
+            present_kind = "rotated (free-slip or fault contact)"
+        else:
+            if not multipliers:
+                return
+            present, refused = len(multipliers), "a rotated constraint"
+            present_kind = "block-constraint multiplier"
+
+        raise RuntimeError(
+            f"{adding}(): this solver already carries {present} "
+            f"{present_kind} boundary condition(s), so it cannot also take "
+            f"{refused}. The rotated solve splits velocity and pressure by "
+            f"field number and a block constraint adds a multiplier field "
+            f"outside that split, so the preconditioner would cover only part "
+            f"of the operator. Both impose the same wall-normal condition — "
+            f"use one of them (issue #464)."
+        )
+
     def add_rotated_freeslip_bc(self, conds=None, boundary=None, normal=None):
         r"""Add STRONG free-slip (:math:`\mathbf{u}\cdot\hat{\mathbf n}=0`) by rotating
         the boundary velocity DOFs into a per-node (normal, tangential) frame and
@@ -6163,6 +6217,8 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
         DeprecationWarning: the string becomes ``boundary`` and a second
         positional argument, if present, becomes ``normal``.
         """
+        self._reject_mixed_constraint_mechanisms("add_rotated_freeslip_bc")
+
         if isinstance(conds, str):
             # legacy boundary-first call: (boundary[, normal])
             if boundary is not None:
@@ -6241,6 +6297,8 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
         """
         from underworld3.utilities import fault_contact
 
+        self._reject_mixed_constraint_mechanisms("add_fault_bc")
+
         if not isinstance(boundary, str):
             raise TypeError(
                 f"add_fault_bc() requires the fault's boundary name string; "
@@ -9344,6 +9402,14 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
         """
 
 
+        # Checked here as well as at registration (#464). The registration check
+        # gives the better message — it knows which call was refused — but it
+        # only covers what goes through the solver's own methods, and
+        # `fault_contact` writes `_fault_contact_faults` directly. This runs
+        # before any setup reads either list, so an unsupported pair costs
+        # nothing before it is refused.
+        self._reject_mixed_constraint_mechanisms("solve")
+
         if homotopy:
             # The march runs a SEQUENCE of ordinary solves at successively sharper
             # yield surfaces; each one re-enters this method with homotopy=False.

src/underworld3/systems/solvers.py

@@ -2659,6 +2659,8 @@ def add_constraint_bc(self, conds=None, boundary=None, normal=None, screening=No
         DeprecationWarning; see
         ``SolverBaseClass._value_first_bc_args``.
         """
+        self._reject_mixed_constraint_mechanisms("add_constraint_bc")
+
         conds, boundary = self._value_first_bc_args(
             "add_constraint_bc", conds, boundary, alias=g)
         g = conds if conds is not None else 0.0

tests/test_0062_rotated_constraint_exclusivity.py

@@ -0,0 +1,82 @@
+"""A solver takes rotated constraints or block constraints, not both (#464).
+
+The rotated driver builds its own index-set fieldsplit over velocity and
+pressure and addresses them by field number. A block constraint registers a
+multiplier field, whose DOFs are then in neither index set, so the
+preconditioner covers a strict subset of the operator's rows and says nothing
+about it. Both mechanisms impose the same wall-normal condition.
+
+The two controls matter as much as the two refusals: a guard that also refused
+each mechanism on its own would pass a test that only checked that mixing
+raises. Neutering `_reject_mixed_constraint_mechanisms` and repeating the first
+case shows what the guard is for — the solver accepts one rotated boundary
+condition and one multiplier together, and says nothing.
+"""
+
+import pytest
+
+import underworld3 as uw
+
+pytestmark = [pytest.mark.level_1, pytest.mark.tier_a]
+
+
+@pytest.fixture
+def stokes():
+    """A constrained Stokes solver, which is the class that has both APIs."""
+
+    mesh = uw.meshing.UnstructuredSimplexBox(
+        minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.5, qdegree=2
+    )
+    v = uw.discretisation.MeshVariable("Uc", mesh, mesh.dim, degree=2)
+    p = uw.discretisation.MeshVariable("Pc", mesh, 1, degree=1)
+
+    return uw.systems.Stokes_Constrained(mesh, velocityField=v, pressureField=p)
+
+
+def test_block_constraint_after_rotated_freeslip_is_refused(stokes):
+    stokes.add_rotated_freeslip_bc(0.0, "Left")
+
+    with pytest.raises(RuntimeError, match=r"add_constraint_bc.*#464"):
+        stokes.add_constraint_bc(0.0, "Bottom")
+
+
+def test_rotated_freeslip_after_block_constraint_is_refused(stokes):
+    stokes.add_constraint_bc(0.0, "Bottom")
+
+    with pytest.raises(RuntimeError, match=r"add_rotated_freeslip_bc.*#464"):
+        stokes.add_rotated_freeslip_bc(0.0, "Left")
+
+
+def test_several_rotated_boundaries_are_still_allowed(stokes):
+    """The control: the guard is about mixing mechanisms, not about counting."""
+
+    stokes.add_rotated_freeslip_bc(0.0, "Left")
+    stokes.add_rotated_freeslip_bc(0.0, "Right")
+
+    assert len(stokes._rotated_freeslip_bcs) == 2
+
+
+def test_several_block_constraints_are_still_allowed(stokes):
+    """The other control."""
+
+    stokes.add_constraint_bc(0.0, "Bottom")
+    stokes.add_constraint_bc(0.0, "Top")
+
+    assert len(stokes._multipliers) == 2
+
+
+def test_the_dispatch_refuses_the_pair_even_if_registration_was_bypassed(stokes):
+    """The guarantee, as opposed to the message.
+
+    The registration checks only cover what goes through the solver's own
+    methods, and `fault_contact` writes `_fault_contact_faults` directly. The
+    solve dispatch is where both lists are read together, so it carries the
+    check as well. Reaching it here means bypassing registration, which is
+    exactly the case the dispatch exists to catch.
+    """
+
+    stokes.add_rotated_freeslip_bc(0.0, "Left")
+    stokes._multipliers.append(object())
+
+    with pytest.raises(RuntimeError, match=r"solve\(\).*#464"):
+        stokes.solve()