Fix inverted hex cell faces in 3D mesh visualization

src/underworld3/visualisation/visualisation.py

@@ -109,6 +109,22 @@ def mesh_to_pv_mesh(mesh, jupyter_backend=None):
         cell_points = mesh.dm.getTransitiveClosure(cell_id)[0][-cell_num_points:]
         cell_points_list.append(cell_points - pStart)
 
+    # PETSc DMPlex uses opposite winding to VTK for 3D cells.
+    # Reorder vertices to match VTK convention (right-handed orientation)
+    # so that face normals, shading, clipping, and back-face culling work correctly.
+    #
+    # Hexahedra: PETSc winds bottom face CCW from above, VTK expects CW from above.
+    #   Fix: swap vertices 1 and 3 on the bottom face → [0,3,2,1,4,5,6,7]
+    # Tetrahedra: PETSc uses left-handed orientation, VTK expects right-handed.
+    #   Fix: swap vertices 1 and 2 → [0,2,1,3]
+    if mesh.dim == 3:
+        if mesh.dm.isSimplex():
+            tet_reorder = [0, 2, 1, 3]
+            cell_points_list = [pts[tet_reorder] for pts in cell_points_list]
+        else:
+            hex_reorder = [0, 3, 2, 1, 4, 5, 6, 7]
+            cell_points_list = [pts[hex_reorder] for pts in cell_points_list]
+
     try:
         import meshio
 

tests/test_0055_mesh_vtk_orientation.py

@@ -0,0 +1,97 @@
+"""Regression test for 3D mesh cell orientation in VTK/meshio output.
+
+PETSc DMPlex uses opposite vertex winding to VTK for 3D cells.
+mesh_to_pv_mesh() must reorder vertices so that cell Jacobians are
+positive (right-handed orientation), otherwise face normals, shading,
+clipping, and back-face culling break.
+
+See: https://github.com/underworldcode/underworld3/issues/40
+"""
+
+import pytest
+import numpy as np
+import underworld3 as uw
+
+pytestmark = [pytest.mark.level_1, pytest.mark.tier_a]
+
+
+def _tet_jacobian_sign(verts):
+    """Signed volume of tetrahedron. Positive = VTK right-handed."""
+    e01 = verts[1] - verts[0]
+    e02 = verts[2] - verts[0]
+    e03 = verts[3] - verts[0]
+    return np.dot(np.cross(e01, e02), e03)
+
+
+def _hex_jacobian_sign(verts):
+    """Trilinear Jacobian sign at vertex 0. Positive = VTK convention."""
+    e01 = verts[1] - verts[0]
+    e03 = verts[3] - verts[0]
+    e04 = verts[4] - verts[0]
+    return np.dot(np.cross(e01, e03), e04)
+
+
+def _get_viz_cell_points(mesh):
+    """Replicate the cell extraction and reordering from mesh_to_pv_mesh."""
+    cStart, cEnd = mesh.dm.getHeightStratum(0)
+    pStart, _ = mesh.dm.getDepthStratum(0)
+    n = mesh.element.entities[mesh.dim]
+
+    cell_points_list = []
+    for cell_id in range(cStart, cEnd):
+        pts = mesh.dm.getTransitiveClosure(cell_id)[0][-n:] - pStart
+        cell_points_list.append(pts)
+
+    if mesh.dim == 3:
+        if mesh.dm.isSimplex():
+            reorder = [0, 2, 1, 3]
+        else:
+            reorder = [0, 3, 2, 1, 4, 5, 6, 7]
+        cell_points_list = [pts[reorder] for pts in cell_points_list]
+
+    return cell_points_list
+
+
+def test_hex_orientation_3d():
+    mesh = uw.meshing.StructuredQuadBox(
+        minCoords=(0.0, 0.0, 0.0),
+        maxCoords=(1.0, 1.0, 1.0),
+        elementRes=(2, 2, 2),
+    )
+    coords = np.asarray(mesh.X.coords, dtype=np.double)
+    cell_points_list = _get_viz_cell_points(mesh)
+
+    for i, pts in enumerate(cell_points_list):
+        det = _hex_jacobian_sign(coords[pts])
+        assert det > 0, f"Hex cell {i} has negative Jacobian ({det})"
+
+
+def test_tet_orientation_3d():
+    mesh = uw.meshing.UnstructuredSimplexBox(
+        minCoords=(0.0, 0.0, 0.0),
+        maxCoords=(1.0, 1.0, 1.0),
+        cellSize=0.5,
+    )
+    coords = np.asarray(mesh.X.coords, dtype=np.double)
+    cell_points_list = _get_viz_cell_points(mesh)
+
+    for i, pts in enumerate(cell_points_list):
+        det = _tet_jacobian_sign(coords[pts])
+        assert det > 0, f"Tet cell {i} has negative Jacobian ({det})"
+
+
+def test_2d_unaffected():
+    """Verify 2D meshes are not modified by the 3D reordering."""
+    for mesh in [
+        uw.meshing.StructuredQuadBox(elementRes=(3, 3)),
+        uw.meshing.UnstructuredSimplexBox(),
+    ]:
+        coords = np.asarray(mesh.X.coords, dtype=np.double)
+        cell_points_list = _get_viz_cell_points(mesh)
+
+        for i, pts in enumerate(cell_points_list):
+            verts = coords[pts]
+            e01 = verts[1] - verts[0]
+            e02 = verts[2] - verts[0]
+            cross = e01[0] * e02[1] - e01[1] * e02[0]
+            assert cross > 0, f"2D cell {i} has wrong winding ({cross})"