Skip to content

Fix XDMF output to explicitly write cell-to-vertex topology and ensure strict ParaView compatibility#205

Open
tiannh7 wants to merge 5 commits into
underworldcode:developmentfrom
tiannh7:fix-xdmf-topology-validation
Open

Fix XDMF output to explicitly write cell-to-vertex topology and ensure strict ParaView compatibility#205
tiannh7 wants to merge 5 commits into
underworldcode:developmentfrom
tiannh7:fix-xdmf-topology-validation

Conversation

@tiannh7
Copy link
Copy Markdown

@tiannh7 tiannh7 commented May 24, 2026

This PR resolves a critical issue where generating XDMF files could silently output PETSc DMPlex's internal point numbering instead of zero-indexed vertex connectivity. This malformed connectivity causes abrupt crashes in ParaView 5.13+ and renders corrupted meshes in older versions.

Key changes:

  • Standard-compliant /viz generation: Ported the _write_mesh_viz_groups generation logic into write_timestep. It now correctly gathers and constructs true /viz/geometry/vertices and /viz/topology/cells across MPI ranks. This restores strict XDMF schema compliance and is 100% backward-compatible with all older versions of ParaView.
  • Proactive validation: Added safety checks in checkpoint_xdmf to explicitly warn users if the connectivity goes out of bounds (cells.max() >= numVertices), if numCorners mismatches the topology type, or if Field attribute lengths (Node/Cell center) are malformed.
  • Unit tests: Added a new test case to enforce that /viz groups are consistently generated and XDMF reference paths correctly point to them without triggering validation warnings.

While this PR hardens our existing XDMF/HDF5 pipeline and ensures stability, it might be worth exploring alternative or complementary formats (such as native .vtu / .vtk) down the road. Integrating tools like pyvista or meshio could potentially offer a more integrated, single-file visualization workflow. This is just a thought for the future.

Copilot AI review requested due to automatic review settings May 24, 2026 15:18
@tiannh7 tiannh7 requested a review from lmoresi as a code owner May 24, 2026 15:18
Copy link
Copy Markdown
Contributor

Copilot AI left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Pull request overview

Note

Copilot was unable to run its full agentic suite in this review.

Adds ParaView-compatible “/viz” HDF5 geometry/topology outputs for meshes written via write_timestep, and validates that the generated XDMF references these groups.

Changes:

  • Write /viz/geometry/vertices and /viz/topology/cells into the mesh HDF5 when create_xdmf is enabled.
  • Add MPI gather-based assembly of global vertex lists and dense connectivity for visualization.
  • Add a regression test ensuring /viz datasets exist and that the XDMF points to them.

Reviewed changes

Copilot reviewed 2 out of 2 changed files in this pull request and generated 9 comments.

File Description
tests/test_0005_xdmf_compat.py Adds a regression test validating /viz datasets and XDMF references.
src/underworld3/discretisation/discretisation_mesh.py Implements /viz mesh group writer, hooks it into write_timestep, and adds extra XDMF warnings/validation.

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

f"to DMPlex vertex count ({n_local_vertices}) for {mesh_h5_path}."
)

local_vertices = numpy.column_stack((vertex_gids, coords_local))
Comment on lines +4120 to +4131
pStart, _ = dm.getDepthStratum(0)
cStart, cEnd = dm.getHeightStratum(0)
vertex_numbering = dm.getVertexNumbering().getIndices()
vertex_gids = _petsc_numbering_to_global_ids(vertex_numbering)
cell_num_points = mesh.element.entities[mesh.dim]

cell_points_list = []
for cell_id in range(cStart, cEnd):
points = numpy.asarray(
dm.getTransitiveClosure(cell_id)[0][-cell_num_points:],
dtype=numpy.int64,
)
Comment on lines +4169 to +4176
gathered_vertices = uw.mpi.comm.gather(local_vertices, root=0)
gathered_cells = uw.mpi.comm.gather(local_cells, root=0)
uw.mpi.barrier()

if uw.mpi.rank == 0:
import h5py

vertex_blocks = [block for block in gathered_vertices if block is not None and block.size > 0]
cell_blocks = [block for block in gathered_cells if block is not None and block.size > 0]

if vertex_blocks:
all_vertices = numpy.vstack(vertex_blocks)
Comment on lines +4189 to +4206
vertex_lookup = {}
for row in all_vertices:
gid = int(row[0])
if gid not in vertex_lookup:
vertex_lookup[gid] = numpy.asarray(row[1:], dtype=numpy.float64)

ordered_gids = numpy.array(sorted(vertex_lookup.keys()), dtype=numpy.int64)
if ordered_gids.size:
ordered_vertices = numpy.vstack([vertex_lookup[int(gid)] for gid in ordered_gids])
else:
ordered_vertices = numpy.empty((0, mesh.dim), dtype=numpy.float64)

dense_lookup = {int(gid): idx for idx, gid in enumerate(ordered_gids.tolist())}
if all_cells.size:
dense_cells = numpy.asarray(
[[dense_lookup[int(gid)] for gid in row] for row in all_cells],
dtype=numpy.int64,
)
Comment on lines +4189 to +4206
vertex_lookup = {}
for row in all_vertices:
gid = int(row[0])
if gid not in vertex_lookup:
vertex_lookup[gid] = numpy.asarray(row[1:], dtype=numpy.float64)

ordered_gids = numpy.array(sorted(vertex_lookup.keys()), dtype=numpy.int64)
if ordered_gids.size:
ordered_vertices = numpy.vstack([vertex_lookup[int(gid)] for gid in ordered_gids])
else:
ordered_vertices = numpy.empty((0, mesh.dim), dtype=numpy.float64)

dense_lookup = {int(gid): idx for idx, gid in enumerate(ordered_gids.tolist())}
if all_cells.size:
dense_cells = numpy.asarray(
[[dense_lookup[int(gid)] for gid in row] for row in all_cells],
dtype=numpy.int64,
)
Comment on lines +4326 to +4334
cells_data = cells[...]
c_min, c_max = cells_data.min(), cells_data.max()
if c_min < 0 or c_max >= numVertices:
warnings.warn(
f"XDMF connectivity is invalid! cells max {c_max} >= "
f"numVertices {numVertices} or min {c_min} < 0. ParaView will likely crash. "
f"Ensure cell-to-vertex connectivity is written."
)

Comment on lines +4321 to +4324
warnings.warn(
"Using raw '/topology/cells' for XDMF. This may not be Paraview-compatible. "
"Expected '/viz/topology/cells'."
)
Comment on lines +4137 to +4142
if mesh.dim == 3:
if 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]
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants