diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 62c411630..8b8735956 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/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): """ diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index 4b11dbd6e..b0985e3b8 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/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)`. diff --git a/src/underworld3/meshing/annulus.py b/src/underworld3/meshing/annulus.py index 1a3631d67..cf8d73e91 100644 --- a/src/underworld3/meshing/annulus.py +++ b/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 diff --git a/src/underworld3/meshing/segmented.py b/src/underworld3/meshing/segmented.py index 4c555191f..b5d264415 100644 --- a/src/underworld3/meshing/segmented.py +++ b/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 diff --git a/src/underworld3/meshing/spherical.py b/src/underworld3/meshing/spherical.py index 732705d86..7d74116db 100644 --- a/src/underworld3/meshing/spherical.py +++ b/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