Add mesh.nullspace_rotations and stokes.petsc_use_nullspace

src/underworld3/cython/petsc_generic_snes_solvers.pyx

@@ -3553,6 +3553,35 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
         self.is_setup = False
         self._saddle_preconditioner = function
 
+    @property
+    def petsc_use_nullspace(self):
+        """Enable full nullspace handling: constant pressure + mesh rotation modes.
+
+        Convenience property that enables the pressure nullspace and
+        auto-populates velocity nullspace modes from ``mesh.nullspace_rotations``.
+
+        For finer control, use ``petsc_use_pressure_nullspace`` and
+        ``petsc_velocity_nullspace_basis`` separately.
+        """
+        return (self._petsc_use_pressure_nullspace
+                or len(self._petsc_velocity_nullspace_basis) > 0)
+
+    @petsc_use_nullspace.setter
+    def petsc_use_nullspace(self, value):
+        value = bool(value)
+        self._petsc_use_pressure_nullspace = value
+        if value:
+            # Sync velocity nullspace basis to mesh rotations (even if empty)
+            # to avoid stale modes carrying over between meshes/runs
+            if hasattr(self.mesh, 'nullspace_rotations'):
+                self.petsc_velocity_nullspace_basis = self.mesh.nullspace_rotations or ()
+            else:
+                self._petsc_velocity_nullspace_basis = ()
+        else:
+            self._petsc_velocity_nullspace_basis = ()
+        self._reset_stokes_nullspace()
+        self.is_setup = False
+
     @property
     def petsc_use_pressure_nullspace(self):
         """

src/underworld3/discretisation/discretisation_mesh.py

@@ -664,6 +664,11 @@ def mesh_update_callback(array, change_context):
         self._lvec = None
         self.petsc_fe = None
 
+        # Rigid-body rotation null modes for this geometry.
+        # Mesh factories override this for closed surfaces (annulus, sphere).
+        # Each entry is a SymPy Matrix velocity field in mesh coordinates.
+        self._nullspace_rotations = []
+
         self.degree = degree
         self.qdegree = qdegree
 
@@ -1531,6 +1536,31 @@ def t(self):
         """
         return self._t
 
+    @property
+    def nullspace_rotations(self):
+        """Symbolic velocity fields for rigid-body rotation null modes.
+
+        Returns a list of SymPy Matrix expressions in mesh Cartesian
+        coordinates. Empty for meshes with no rotation nullspace (boxes,
+        wedge segments with walls). Set by mesh factory functions for
+        closed surfaces (annulus, spherical shell, etc.).
+
+        Each entry represents a rigid rotation: v = omega x r.
+
+        Returns
+        -------
+        list of sympy.Matrix
+            Velocity fields for each independent rotation mode.
+
+        Examples
+        --------
+        >>> annulus = uw.meshing.Annulus(...)
+        >>> annulus.nullspace_rotations  # [Matrix([-y, x])]
+        >>> shell = uw.meshing.SphericalShell(...)
+        >>> shell.nullspace_rotations   # 3 rotation matrices
+        """
+        return self._nullspace_rotations
+
     @property
     def r(self) -> Tuple[sympy.vector.BaseScalar]:
         r"""Tuple of coordinate scalars :math:`(x, y)` or :math:`(x, y, z)`.

src/underworld3/meshing/annulus.py

@@ -538,6 +538,10 @@ class boundary_normals(Enum):
 
     new_mesh.boundary_normals = boundary_normals
 
+    # Full annulus: rigid rotation about z-axis
+    x, y = new_mesh.X
+    new_mesh._nullspace_rotations = [sympy.Matrix([-y, x])]
+
     return new_mesh
 
 
@@ -1110,6 +1114,10 @@ class boundary_normals(Enum):
 
     new_mesh.boundary_normals = boundary_normals
 
+    # Full annulus with spokes: rigid rotation about z-axis
+    x, y = new_mesh.X
+    new_mesh._nullspace_rotations = [sympy.Matrix([-y, x])]
+
     return new_mesh
 
 
@@ -1392,6 +1400,10 @@ class boundary_normals(Enum):
 
     new_mesh.boundary_normals = boundary_normals
 
+    # Full annulus with internal boundary: rigid rotation about z-axis
+    x, y = new_mesh.X
+    new_mesh._nullspace_rotations = [sympy.Matrix([-y, x])]
+
     return new_mesh
 
 
@@ -1671,4 +1683,8 @@ class boundary_normals(Enum):
 
     new_mesh.boundary_normals = boundary_normals
 
+    # Full disc with internal boundaries: rigid rotation about z-axis
+    x, y = new_mesh.X
+    new_mesh._nullspace_rotations = [sympy.Matrix([-y, x])]
+
     return new_mesh

src/underworld3/meshing/segmented.py

@@ -676,6 +676,14 @@ class boundary_normals(Enum):
 
     new_mesh.boundary_normals = boundary_normals
 
+    # Full segmented spherical shell: 3 rigid rotation modes
+    x, y, z = new_mesh.X
+    new_mesh._nullspace_rotations = [
+        sympy.Matrix([0, -z, y]),
+        sympy.Matrix([z, 0, -x]),
+        sympy.Matrix([-y, x, 0]),
+    ]
+
     return new_mesh
 
 
@@ -1088,4 +1096,12 @@ class boundary_normals(Enum):
 
     new_mesh.boundary_normals = boundary_normals
 
+    # Solid sphere: 3 rigid rotation modes
+    x, y, z = new_mesh.X
+    new_mesh._nullspace_rotations = [
+        sympy.Matrix([0, -z, y]),
+        sympy.Matrix([z, 0, -x]),
+        sympy.Matrix([-y, x, 0]),
+    ]
+
     return new_mesh

src/underworld3/meshing/spherical.py

@@ -269,6 +269,14 @@ class boundary_normals(Enum):
         Upper = 12
         Centre = 1
 
+    # Full spherical shell: 3 rigid rotation modes
+    x, y, z = new_mesh.X
+    new_mesh._nullspace_rotations = [
+        sympy.Matrix([0, -z, y]),   # rotation about x
+        sympy.Matrix([z, 0, -x]),   # rotation about y
+        sympy.Matrix([-y, x, 0]),   # rotation about z
+    ]
+
     return new_mesh
 
 
@@ -475,6 +483,14 @@ class boundary_normals(Enum):
         Upper = 13
         Centre = 1
 
+    # Full spherical shell with internal boundary: 3 rigid rotation modes
+    x, y, z = new_mesh.X
+    new_mesh._nullspace_rotations = [
+        sympy.Matrix([0, -z, y]),
+        sympy.Matrix([z, 0, -x]),
+        sympy.Matrix([-y, x, 0]),
+    ]
+
     return new_mesh
 
 
@@ -1011,4 +1027,12 @@ class boundary_normals(Enum):
 
     new_mesh.boundary_normals = boundary_normals
 
+    # Full cubed sphere: 3 rigid rotation modes
+    x, y, z = new_mesh.X
+    new_mesh._nullspace_rotations = [
+        sympy.Matrix([0, -z, y]),
+        sympy.Matrix([z, 0, -x]),
+        sympy.Matrix([-y, x, 0]),
+    ]
+
     return new_mesh