Fix DG1 XDMF export with element-local simplex topology

docs/developer/subsystems/checkpointing-system.md

@@ -30,13 +30,31 @@ Optional payloads are controlled by explicit flags:
 
 ### Visualisation and Coordinate Remap
 
+Continuous fields use the standard mesh vertices, and DG0 fields use cell data.
+DG1 fields on full-dimensional triangular and tetrahedral meshes use a second
+grid named `DG1` in the same XDMF file. Each simplex has its own three or four
+physical vertices: the saved linear polynomial is evaluated within that cell,
+without averaging traces across shared edges or faces. Interior DG interpolation
+nodes are not mistaken for the physical mesh vertices.
+
+The DG1 visualization arrays (`vertices`, `cells`, `values`) live under `/dg1`
+in each variable HDF5 file. Tensor visualization uses a nine-component 3-by-3
+layout (zero-padded in 2D); `/fields` and PETSc reload data retain their native
+layout and precision. Open the one `.xdmf` file in ParaView and select the
+`domain` or `DG1` block for the corresponding fields. Do not merge coincident
+points or apply point-averaging filters if discontinuities must be preserved.
+
+Higher-degree discontinuous fields, tensor-product DG cells, embedded manifolds,
+and integration-point fields are not supported by this DG1 exporter. They raise
+an explicit error when visualization is requested; `create_xdmf=False` still
+allows native checkpoint output. Parallel export uses owned cells only.
+
 ```python
 mesh.write_timestep(
     "output",
     index=100,
     outputPath="output",
     meshVars=[velocity, pressure, temperature],
-    time=100.0,
     create_xdmf=True,
 )
 ```

src/underworld3/discretisation/discretisation_mesh.py

@@ -4877,11 +4877,31 @@ def write_timestep(
             exists. If ``True``, write an indexed mesh file for this timestep.
         create_xdmf
             Write ParaView/XDMF-compatible datasets and companion XDMF file.
+            DG1 on full-dimensional triangles/tetrahedra uses a separate grid
+            with independent vertices per cell, preserving jumps without
+            smoothing. Visualization-only arrays live under ``/dg1`` in the
+            variable files; native checkpoint and reload data are unchanged.
+            Higher-order discontinuous and non-simplex DG visualization are
+            not supported (use ``create_xdmf=False`` for native-only output).
         petsc_reload
             Write PETSc DMPlex section/vector metadata for reload with
             ``MeshVariable.read_checkpoint()``.
 
         """
+        if create_xdmf:
+            for var in meshVars or []:
+                integration_point = getattr(var, "is_integration_point", False)
+                if integration_point or (not var.continuous and var.degree > 0):
+                    if (
+                        var.degree != 1 or not self.isSimplex
+                        or self.dim not in (2, 3) or self.cdim != self.dim
+                        or integration_point
+                    ):
+                        raise NotImplementedError(
+                            "DG XDMF supports degree-one fields on full-dimensional "
+                            "triangle/tetrahedron meshes only; use create_xdmf=False "
+                            "for native-only checkpoints."
+                        )
         options = PETSc.Options()
         options.setValue("viewer_hdf5_sp_output", True)
         options.setValue("viewer_hdf5_collective", False)
@@ -9632,6 +9652,8 @@ def _write_compat_groups(mesh, var, var_h5_path):
     Uses ``uw.function.write_vertices_to_viewer`` (PETSc interpolation +
     ViewerHDF5) for continuous variables, and
     ``uw.function.write_cell_field_to_viewer`` for cell/DG-0 variables.
+    DG1 uses ``/dg1`` with disconnected simplex vertices and nodal traces,
+    never the one-value-per-cell compatibility path.
     PETSc handles all parallel I/O natively.
 
     Vertex coordinates are also written to ``/vertex_fields/coordinates``
@@ -9649,8 +9671,9 @@ def _write_compat_groups(mesh, var, var_h5_path):
     """
     import underworld3 as uw
 
-    is_cell = (not var.continuous) or (var.degree == 0)
-    group = "cell_fields" if is_cell else "vertex_fields"
+    is_dg1 = not var.continuous and var.degree == 1
+    is_cell = var.degree == 0
+    group = "dg1" if is_dg1 else ("cell_fields" if is_cell else "vertex_fields")
 
     # Some PETSc versions (3.21+) write /vertex_fields/ or /cell_fields/
     # automatically during var.write().  Remove any pre-existing group so
@@ -9668,14 +9691,29 @@ def _write_compat_groups(mesh, var, var_h5_path):
         var_h5_path, "a", comm=PETSc.COMM_WORLD,
     )
 
-    if is_cell:
+    if is_dg1:
+        from underworld3.function.field_projection import _write_dg1_to_viewer
+        _write_dg1_to_viewer(var, viewer)
+    elif is_cell:
         uw.function.write_cell_field_to_viewer(var, viewer)
     else:
         uw.function.write_vertices_to_viewer(var, viewer)
         uw.function.write_coordinates_to_viewer(mesh, viewer)
 
     viewer.destroy()
 
+    if is_dg1:
+        # Only topology is generated on rank zero; field values and vertices
+        # were written collectively by PETSc in the same owned-cell order.
+        if uw.mpi.rank == 0:
+            with h5py.File(var_h5_path, "a") as f:
+                nvertices = f["dg1/vertices"].shape[0]
+                ncorners = mesh.dim + 1
+                f["dg1"].create_dataset(
+                    "cells", data=numpy.arange(nvertices, dtype=numpy.int64).reshape(-1, ncorners)
+                )
+        uw.mpi.barrier()
+
 
 def checkpoint_xdmf(
     filename: str,
@@ -9802,6 +9840,13 @@ def checkpoint_xdmf(
     header += """
 ]>"""
 
+    dg_vars = [var for var in meshVars if not var.continuous and var.degree == 1]
+    collection_start = (
+        '<Grid Name="fields" GridType="Collection" CollectionType="Spatial">'
+        f'<Time Value="{index}"/>'
+        if dg_vars else ""
+    )
+    collection_end = "</Grid>" if dg_vars else ""
     xdmf_start = f"""
 <Xdmf>
   <Domain Name="domain">
@@ -9818,6 +9863,7 @@ def checkpoint_xdmf(
       &MeshData;:/{geomPath}/vertices
     </DataItem>
     <!-- ============================================================ -->
+      {collection_start}
       <Grid Name="domain" GridType="Uniform">
         <Topology
            TopologyType="{topology_type}"
@@ -9862,6 +9908,8 @@ def get_field_info(h5_filename, mesh_var, center):
 
     attributes = ""
     for var in meshVars:
+        if not var.continuous and var.degree == 1:
+            continue
         var_filename = filename + f".mesh.{var.clean_name}.{index:05}.h5"
 
         # Determine if data is stored on nodes (vertex_fields) or cells (cell_fields)
@@ -9946,8 +9994,54 @@ def get_field_info(h5_filename, mesh_var, center):
     """
         attributes += var_attribute
 
+    dg_grid = ""
+    if dg_vars:
+        first = dg_vars[0]
+        first_filename = filename + f".mesh.{first.clean_name}.{index:05}.h5"
+        with h5py.File(first_filename, "r") as f:
+            dg_cells = f["dg1/cells"].shape
+            dg_points = f["dg1/vertices"].shape
+        if dg_cells != (numCells, numCorners) or dg_points != (numCells * numCorners, spaceDim):
+            raise ValueError("DG1 visualization topology does not match the checkpoint mesh")
+        dg_grid = f"""
+    <Grid Name="DG1" GridType="Uniform">
+      <Topology TopologyType="{topology_type}" NumberOfElements="{numCells}">
+        <DataItem Format="HDF" NumberType="Int" Precision="8" Dimensions="{numCells} {numCorners}">
+          &{first.clean_name}_Data;:/dg1/cells
+        </DataItem>
+      </Topology>
+      <Geometry GeometryType="{geomType}">
+        <DataItem Format="HDF" NumberType="Float" Precision="8" Dimensions="{dg_points[0]} {spaceDim}">
+          &{first.clean_name}_Data;:/dg1/vertices
+        </DataItem>
+      </Geometry>
+"""
+        for var in dg_vars:
+            var_filename = filename + f".mesh.{var.clean_name}.{index:05}.h5"
+            with h5py.File(var_filename, "r") as f:
+                shape = f["dg1/values"].shape
+            if shape[0] != dg_points[0]:
+                raise ValueError(f"DG1 visualization size mismatch for {var.clean_name}")
+            components = shape[1] if len(shape) == 2 else 1
+            if var.vtype in (uw.VarType.TENSOR, uw.VarType.SYM_TENSOR):
+                kind = "Tensor"
+            elif var.vtype == uw.VarType.MATRIX:
+                kind = "Matrix"
+            else:
+                kind = "Scalar" if components == 1 else "Vector"
+            dimensions = " ".join(str(value) for value in shape)
+            dg_grid += f"""
+      <Attribute Name="{var.clean_name}" AttributeType="{kind}" Center="Node">
+        <DataItem Format="HDF" NumberType="Float" Precision="8" Dimensions="{dimensions}">
+          &{var.clean_name}_Data;:/dg1/values
+        </DataItem>
+      </Attribute>
+"""
+        dg_grid += "    </Grid>"
     xdmf_end = f"""
     </Grid>
+    {dg_grid}
+    {collection_end}
   </Domain>
 </Xdmf>
     """

src/underworld3/function/field_projection.py

@@ -350,3 +350,43 @@ def write_cell_field_to_viewer(
     mesh_var._sync_lvec_to_gvec()
     data = mesh_var._gvec.array.reshape(-1, nc).copy()
     _write_vec_to_group(viewer, data, name, group, PETSc.COMM_WORLD)
+
+
+def _write_dg1_to_viewer(mesh_var, viewer):
+    """Write owned simplex cells with independent physical vertices and DG1 traces.
+
+    Coordinate-section cell maps preserve element ownership and node ordering;
+    no point location, coordinate matching, or inter-element averaging is used.
+    Interior DG interpolation nodes define an affine polynomial, evaluated at
+    that same cell's physical vertices. Native checkpoint vectors are untouched.
+    """
+    mesh = mesh_var.mesh
+    if (
+        mesh_var.continuous or mesh_var.degree != 1 or not mesh.isSimplex
+        or mesh.dim not in (2, 3) or mesh.cdim != mesh.dim
+    ):
+        raise NotImplementedError("DG1 XDMF requires a full-dimensional triangle/tetrahedron mesh")
+    cstart, cend = mesh.dm.getHeightStratum(0)
+    owned = np.ones(cend - cstart, dtype=bool)
+    # Serial DMPlex may have an unset point SF; no cells are ghosts there.
+    if mesh.dm.comm.getSize() > 1:
+        _, leaves, remote = mesh.dm.getPointSF().getGraph()
+        if leaves is None:
+            leaves = np.arange(len(remote))
+        leaves = np.asarray(leaves)
+        owned[leaves[(leaves >= cstart) & (leaves < cend)] - cstart] = False
+    rows = mesh._cell_node_indices(1, False).reshape(-1, mesh.dim + 1)[owned]
+    vertex_rows = mesh._cell_node_indices(1, True).reshape(-1, mesh.dim + 1)[owned]
+    corners = mesh._get_coords_for_basis(1, True)[vertex_rows]
+    # Closure order is arbitrary; give VTK positively oriented simplices.
+    negative = np.linalg.det((corners[:, 1:] - corners[:, :1]).transpose(0, 2, 1)) < 0
+    corners[negative] = corners[negative][:, [0, 2, 1] if mesh.dim == 2 else [0, 2, 1, 3]]
+    nodes = mesh_var.coords[rows]
+    coefficients = mesh_var._lvec.array.reshape(-1, mesh_var.num_components)[rows]
+    matrix = (nodes[:, 1:] - nodes[:, :1]).transpose(0, 2, 1)
+    local = np.linalg.solve(matrix, (corners - nodes[:, :1]).transpose(0, 2, 1))
+    weights = np.concatenate((1 - local.sum(axis=1, keepdims=True), local), axis=1)
+    values = np.einsum("cij,cik->cjk", weights, coefficients).reshape(-1, mesh_var.num_components)
+    values = _repack_tensor_to_paraview(values, mesh_var.vtype, mesh.dim)
+    _write_vec_to_group(viewer, corners.reshape(-1, mesh.cdim), "vertices", "/dg1", PETSc.COMM_WORLD)
+    _write_vec_to_group(viewer, values, "values", "/dg1", PETSc.COMM_WORLD)

tests/test_0005_xdmf_dg1.py

@@ -0,0 +1,142 @@
+"""DG1 visualization preserves affine fields and jumps, including MPI ownership."""
+
+from pathlib import Path
+import xml.etree.ElementTree as ET
+
+import h5py
+import numpy as np
+import pytest
+import underworld3 as uw
+
+
+@pytest.mark.level_1
+@pytest.mark.tier_b
+@pytest.mark.parametrize("dim", [2, 3])
+def test_dg1_simplex_output(tmp_path, dim):
+    directory = Path(uw.mpi.comm.bcast(str(tmp_path), root=0))
+    mesh = uw.meshing.UnstructuredSimplexBox(
+        minCoords=(0.0,) * dim,
+        maxCoords=(1.0,) * dim,
+        cellSize=0.5,
+        regular=True,
+        qdegree=3,
+    )
+    scalar = uw.discretisation.MeshVariable("dg_scalar", mesh, 1, degree=1, continuous=False)
+    tensor = uw.discretisation.MeshVariable(
+        "dg_tensor",
+        mesh,
+        degree=1,
+        continuous=False,
+        vtype=uw.VarType.TENSOR,
+    )
+    pressure = uw.discretisation.MeshVariable("pressure", mesh, 1, degree=1)
+    vector = uw.discretisation.MeshVariable("dg_vector", mesh, dim, degree=1, continuous=False)
+    symmetric = uw.discretisation.MeshVariable(
+        "dg_symmetric", mesh, degree=1, continuous=False, vtype=uw.VarType.SYM_TENSOR
+    )
+    rows = mesh._cell_node_indices(1, False).reshape(-1, dim + 1)
+    coords = scalar.coords
+    offset = np.floor(coords[rows].mean(axis=1)[:, 0] * 7 + 1e-8)
+    scalar.array[rows, 0, 0] = 1 + coords[rows, 0] + 2 * coords[rows, 1] + offset[:, None]
+    tensor.array[:] = 0
+    tensor.array[:, 0, 0] = scalar.array[:, 0, 0]
+    tensor.array[:, 0, 1] = 3 + coords[:, 0]
+    tensor.array[:, 1, 0] = -2 + coords[:, 1]
+    tensor.array[:, 1, 1] = 5
+    pressure.array[:, 0, 0] = pressure.coords[:, 0]
+    vector.array[:, 0, :] = coords
+    symmetric.array[:] = 0
+    symmetric.array[:, 0, 0] = 2
+    symmetric.array[:, 1, 1] = 3
+    symmetric.array[:, 0, 1] = coords[:, 0]
+    original = np.array(scalar.array)
+    mesh.write_timestep(
+        "fields",
+        index=0,
+        outputPath=str(directory),
+        meshVars=[pressure, scalar, tensor, vector, symmetric],
+        petsc_reload=True,
+    )
+    restored = uw.discretisation.MeshVariable("restored", mesh, 1, degree=1, continuous=False)
+    restored.read_checkpoint(
+        str(directory / "fields.mesh.dg_scalar.00000.h5"), data_name="dg_scalar", same_layout=True
+    )
+    np.testing.assert_allclose(restored.array, original, rtol=1e-12, atol=1e-12)
+    if uw.mpi.rank != 0:
+        return
+    with h5py.File(directory / "fields.mesh.dg_scalar.00000.h5", "r") as handle:
+        points = handle["dg1/vertices"][:]
+        cells = handle["dg1/cells"][:]
+        values = handle["dg1/values"][:].reshape(-1)
+        native = handle["fields/dg_scalar"][:]
+        assert len(points) == len(cells) * (dim + 1)
+        assert len(np.unique(cells)) == len(points)
+        assert native.size == len(points)
+    with h5py.File(directory / "fields.mesh.00000.h5", "r") as handle:
+        assert len(cells) == len(handle["viz/topology/cells"])
+    corners = points[cells]
+    assert np.all(np.linalg.det((corners[:, 1:] - corners[:, :1]).transpose(0, 2, 1)) > 0)
+    centers = corners.mean(axis=1)
+    assert len(np.unique(np.round(centers, 10), axis=0)) == len(cells)
+    offset = np.repeat(np.floor(centers[:, 0] * 7 + 1e-8), dim + 1)
+    expected = 1 + points[:, 0] + 2 * points[:, 1] + offset
+    np.testing.assert_allclose(values, expected, rtol=1e-12, atol=1e-12)
+    # Repeated positions can have different values: these traces must not merge.
+    _, inverse = np.unique(np.round(points, 10), axis=0, return_inverse=True)
+    low = np.full(inverse.max() + 1, np.inf)
+    high = np.full(inverse.max() + 1, -np.inf)
+    np.minimum.at(low, inverse, values)
+    np.maximum.at(high, inverse, values)
+    assert np.max(high - low) >= 1
+    with h5py.File(directory / "fields.mesh.dg_tensor.00000.h5", "r") as handle:
+        tensor_values = handle["dg1/values"][:]
+        np.testing.assert_allclose(handle["dg1/vertices"][:], points)
+        assert tensor_values.shape == (len(points), 9)
+        np.testing.assert_allclose(tensor_values[:, 0], expected)
+        np.testing.assert_allclose(tensor_values[:, 1], 3 + points[:, 0])
+        np.testing.assert_allclose(tensor_values[:, 3], -2 + points[:, 1], atol=1e-12)
+        np.testing.assert_allclose(tensor_values[:, 4], 5)
+    tree = ET.parse(directory / "fields.mesh.00000.xdmf")
+    with h5py.File(directory / "fields.mesh.dg_vector.00000.h5", "r") as handle:
+        np.testing.assert_allclose(handle["dg1/values"][:], points, atol=1e-12)
+    with h5py.File(directory / "fields.mesh.dg_symmetric.00000.h5", "r") as handle:
+        sym_values = handle["dg1/values"][:]
+        assert sym_values.shape == (len(points), 9)
+        np.testing.assert_allclose(sym_values[:, 0], 2)
+        np.testing.assert_allclose(sym_values[:, 4], 3)
+        np.testing.assert_allclose(sym_values[:, 1], points[:, 0], atol=1e-12)
+        np.testing.assert_allclose(sym_values[:, 3], points[:, 0], atol=1e-12)
+    grids = tree.findall(".//Grid[@GridType='Uniform']")
+    assert len(grids) == 2
+    dg = next(grid for grid in grids if grid.get("Name") == "DG1")
+    assert {a.get("Name") for a in dg.findall("Attribute")} == {
+        "dg_scalar",
+        "dg_tensor",
+        "dg_vector",
+        "dg_symmetric",
+    }
+    assert all(a.get("Center") == "Node" for a in dg.findall("Attribute"))
+    for item in tree.findall(".//DataItem[@Format='HDF']"):
+        filename, dataset = item.text.strip().split(":", 1)
+        with h5py.File(directory / filename, "r") as handle:
+            assert tuple(map(int, item.get("Dimensions").split())) == handle[dataset].shape
+
+
+@pytest.mark.level_1
+@pytest.mark.tier_b
+@pytest.mark.parametrize("degree", [1, 2])
+def test_unsupported_dg_layout_fails_before_writing(tmp_path, degree):
+    directory = Path(uw.mpi.comm.bcast(str(tmp_path), root=0))
+    mesh = (
+        uw.meshing.StructuredQuadBox(elementRes=(2, 2))
+        if degree == 1
+        else uw.meshing.UnstructuredSimplexBox(cellSize=0.5, regular=True)
+    )
+    dg = uw.discretisation.MeshVariable("dg", mesh, 1, degree=degree, continuous=False)
+    with pytest.raises(NotImplementedError, match="DG.*XDMF"):
+        mesh.write_timestep("unsupported", index=0, outputPath=str(directory), meshVars=[dg])
+    assert not (directory / "unsupported.mesh.00000.h5").exists()
+    # Unsupported visualization must not prevent native-only checkpoints.
+    mesh.write_timestep(
+        "native", index=0, outputPath=str(directory), meshVars=[dg], create_xdmf=False
+    )