From cce29a462284441fa1a6be02664642b5881e12fa Mon Sep 17 00:00:00 2001 From: Tyagi Date: Sat, 12 Sep 2026 19:33:30 +1000 Subject: [PATCH] Fix DG1 XDMF export with element-local simplex topology Write independent physical vertices and cell-local linear field traces for DG1 triangles and tetrahedra. Use native cell maps and owned-cell PETSc output to preserve discontinuities without coordinate matching or inter-element averaging. Expose continuous and DG1 fields as separate grids in one XDMF file, support scalar/vector/tensor layouts, and retain native checkpoint data unchanged. Reject unsupported DG visualization explicitly while preserving native-only output. Add serial and MPI regression coverage for affine fields, jumps, orientation, ownership, tensor packing, and checkpoint reload. Document supported layouts and ParaView usage. Validated 24 focused serial tests, four DG tests on eight ranks, and actual ParaView 2D/3D reads. --- .../subsystems/checkpointing-system.md | 20 ++- .../discretisation/discretisation_mesh.py | 100 +++++++++++- src/underworld3/function/field_projection.py | 40 +++++ tests/test_0005_xdmf_dg1.py | 142 ++++++++++++++++++ 4 files changed, 298 insertions(+), 4 deletions(-) create mode 100644 tests/test_0005_xdmf_dg1.py diff --git a/docs/developer/subsystems/checkpointing-system.md b/docs/developer/subsystems/checkpointing-system.md index c0502b894..c1445ba41 100644 --- a/docs/developer/subsystems/checkpointing-system.md +++ b/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, ) ``` diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index 7f045997e..679f025b4 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/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,7 +9691,10 @@ 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) @@ -9676,6 +9702,18 @@ def _write_compat_groups(mesh, var, var_h5_path): 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 = ( + '' + f'" if dg_vars else "" xdmf_start = f""" @@ -9818,6 +9863,7 @@ def checkpoint_xdmf( &MeshData;:/{geomPath}/vertices + {collection_start} + + + &{first.clean_name}_Data;:/dg1/cells + + + + + &{first.clean_name}_Data;:/dg1/vertices + + +""" + 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""" + + + &{var.clean_name}_Data;:/dg1/values + + +""" + dg_grid += " " xdmf_end = f""" + {dg_grid} + {collection_end} """ diff --git a/src/underworld3/function/field_projection.py b/src/underworld3/function/field_projection.py index 28fc4809d..91a01b2cc 100644 --- a/src/underworld3/function/field_projection.py +++ b/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) diff --git a/tests/test_0005_xdmf_dg1.py b/tests/test_0005_xdmf_dg1.py new file mode 100644 index 000000000..163aea3ca --- /dev/null +++ b/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 + )