High-order coordinates via interpolated DMPlex

ngsPETSc/plex.py

@@ -16,7 +16,7 @@ class ngs:
         "dummy class"
         class comp:
             "dummy class"
-            Mesh = type(None)
+            Mesh = type("_MissingNGSolveMesh", (), {})
 
 FACE_SETS_LABEL = "Face Sets"
 CELL_SETS_LABEL = "Cell Sets"
@@ -38,21 +38,31 @@ def __init__(self, mesh, comm=None, geo=None, name="Default"):
         elif isinstance(comm, PETSc.Comm):
             comm = comm.tompi4py()
 
+        source_type = None
         if isinstance(mesh, ngs.comp.Mesh):
             mesh = mesh.ngmesh
-
-        if isinstance(mesh, ngm.Mesh):
+        if comm.rank == 0:
+            if isinstance(mesh, ngm.Mesh):
+                source_type = "netgen"
+            elif isinstance(mesh, PETSc.DMPlex):
+                source_type = "plex"
+        source_type = comm.bcast(source_type, root=0)
+
+        if source_type == "netgen":
             ngmesh = mesh
             plex = createPETScDMPlex(ngmesh, comm, name)
-        elif isinstance(mesh, PETSc.DMPlex):
+        elif source_type == "plex":
             plex = mesh
             ngmesh = createNetgenMesh(plex, geo)
         else:
             raise TypeError("Mesh format not recognised.")
         self.petscPlex = plex
         self.ngMesh = ngmesh
         self.comm = comm
-        self.geo = self.ngMesh.GetGeometry()
+        self.geo = (
+            self.ngMesh.GetGeometry()
+            if source_type == "plex" or comm.rank == 0 else None
+        )
         self.geoInfo = bool(self.geo)
 
 
@@ -180,23 +190,28 @@ def createPETScDMPlex(ngMesh, comm, name):
     :arg ngMesh: the serial Netgen mesh object to be converted
     :arg comm: the MPI.Comm object
 
-    :returns: a tuple of Netgen mesh and DMPlex
+    :returns: the interpolated PETSc DMPlex
     """
-    if len(ngMesh.GetIdentifications()) > 0:
-        warnings.warn("Periodic meshes are not supported by ngsPETSc" , RuntimeWarning)
-    els = {
-        0: ngMesh.Elements0D,
-        1: ngMesh.Elements1D,
-        2: ngMesh.Elements2D,
-        3: ngMesh.Elements3D,
-    }
-    gdim = ngMesh.dim
-    tdim = gdim
-    cells = els[tdim]()
-    while len(cells) == 0 and tdim > 0:
-        tdim -= 1
+    if comm.rank == 0:
+        els = {
+            0: ngMesh.Elements0D,
+            1: ngMesh.Elements1D,
+            2: ngMesh.Elements2D,
+            3: ngMesh.Elements3D,
+        }
+        if len(ngMesh.GetIdentifications()) > 0:
+            warnings.warn("Periodic meshes are not supported by ngsPETSc", RuntimeWarning)
+        gdim = ngMesh.dim
+        tdim = gdim
         cells = els[tdim]()
-    tdim = comm.bcast(tdim, root=0)
+        while len(cells) == 0 and tdim > 0:
+            tdim -= 1
+            cells = els[tdim]()
+    else:
+        gdim = None
+        tdim = None
+        cells = None
+    gdim, tdim = comm.bcast((gdim, tdim), root=0)
     if comm.rank == 0:
         cells_np = cells.NumPy()
         # Netgen always stores coordinates as float64. createFromCellList performs

tests/test_plex.py

@@ -7,6 +7,7 @@
 from netgen.csg import unit_cube
 
 from petsc4py import PETSc
+from mpi4py import MPI
 import pytest
 
 from ngsPETSc import MeshMapping
@@ -18,6 +19,23 @@ def _plex_number_of_points(plex, h=0, local=False):
         np = plex.getComm().tompi4py().allreduce(np)
     return np
 
+
+def test_ngs_plex_root_only():
+    """Construct and refine an interpolated DMPlex from rank-zero Netgen data."""
+    comm = MPI.COMM_WORLD
+    mesh = unit_square.GenerateMesh(maxh=1.) if comm.rank == 0 else None
+    plex = MeshMapping(mesh, comm=comm).petscPlex
+
+    assert plex.getDimension() == 2
+    assert _plex_number_of_points(plex, h=0) > 0
+    assert _plex_number_of_points(plex, h=1) > 0
+
+    coarse_cells = _plex_number_of_points(plex, h=0)
+    plex.setRefinementUniform(True)
+    refined = plex.refine()
+    assert _plex_number_of_points(refined, h=0) == 4 * coarse_cells
+
+
 @pytest.mark.mpi_skip
 def test_ngs_plex_2d():
     '''