[BUG] - Navier-Stokes solver produces PETSc error when run in PARALLEL if using a regular UnstructuredSimplexBox or StructuredQuadBox mesh

[jcgraciosa]

Describe the bug
The Navier-Stokes solver produces the following "PETSc has generated inconsistent data" error when run in parallel using a REGULAR UnstructuredSimplexBox and StructuredQuadBox:

[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
[0]PETSC ERROR: PETSc has generated inconsistent data
[0]PETSC ERROR: Point 341: 2.0625 -6.32996e-25 0. not located in mesh
[0]PETSC ERROR: WARNING! There are unused option(s) set! Could be the program crashed before usage or a spelling mistake, etc!

• What version of the code?
  Latest development branch
• What system are you running the code on?
  Can be reproduced in a local Mac installation
• Steps to reproduce?
  This stripped version of Couette flow can be used to replicate the issue:

import sympy
import underworld3 as uw

Re         = 100
resolution = 16
U_wall     = 1.0
rho        = 1.0
nu         = 1.0 / Re
mu         = rho * nu

width       = 4.0
height      = 1.0

minX, maxX  = 0., width
minY, maxY  = 0., height

uw.pprint(minX, minY)
uw.pprint(maxX, maxY)

# interpolation bug occurs with regular set to True
mesh = uw.meshing.UnstructuredSimplexBox(
    minCoords=(minX, minY),
    maxCoords=(maxX, maxY),
    cellSize=1. / resolution,
    qdegree=3,
    regular=True,
)

# or with quadbox
#mesh = uw.meshing.StructuredQuadBox(
#   minCoords=(minX, minY),
#   maxCoords=(maxX, maxY),
#   elementRes=(resolution, resolution),
#   qdegree=3,
#)

mesh.dm.view()

v = uw.discretisation.MeshVariable("U", mesh, mesh.dim, degree=2)
p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1, continuous=False)

ns = uw.systems.NavierStokesSLCN(
    mesh, velocityField=v, pressureField=p, rho=rho, order=1
)
ns.constitutive_model = uw.constitutive_models.ViscousFlowModel
ns.constitutive_model.Parameters.viscosity = mu
ns.bodyforce = sympy.Matrix([0, 0])

ns.add_dirichlet_bc((U_wall, 0.0), "Top")
ns.add_dirichlet_bc((0.0,   0.0), "Bottom")

ns.tolerance = 1.0e-6

ns.petsc_options["snes_monitor"] = None
ns.petsc_options["snes_converged_reason"] = None
ns.petsc_options["ksp_monitor"] = None

dx = mesh.get_min_radius()
dt_adv = 1.0 * dx / U_wall
dt_diff = 1.0 * dx**2 / nu
dt = min(dt_adv, dt_diff)

uw.pprint(f"dt: {dt}")

ns.solve(timestep=dt, zero_init_guess=True)

The issue can be reproduced with 3 procs:

pixi run mpirun -n 3 python3 couette-stripped.py