Adding a netgen flag to specify the index value associated with each label

ngsPETSc/plex.py

@@ -33,11 +33,11 @@ class MeshMapping:
 
     '''
 
-    def __init__(self, mesh=None, comm=MPI.COMM_WORLD, name="Default"):
+    def __init__(self, mesh=None, comm=MPI.COMM_WORLD, name="default", labels={}): #pylint: disable=W0102
         self.name = name
         self.comm = comm
         if isinstance(mesh,(ngs.comp.Mesh,ngm.Mesh)):
-            self.createPETScDMPlex(mesh)
+            self.createPETScDMPlex(mesh, labels=labels)
         elif isinstance(mesh,PETSc.DMPlex):
             self.createNGSMesh(mesh)
         else:
@@ -146,19 +146,27 @@ def createNGSMesh(self, plex):
         else:
             raise NotImplementedError("No implementation for dimension greater than 3.")
 
-    def createPETScDMPlex(self, mesh):
+    def createPETScDMPlex(self, mesh, labels):
         '''
         This function generate an PETSc DMPlex from a Netgen/NGSolve mesh object
 
-        :arg plex: the Netgen/NGSolve mesh object to be converted into a PETSc
+        :arg mesh: the Netgen/NGSolve mesh object to be converted into a PETSc
                    DMPlex.
-
+        :arg labels: Dictionary with the new labels for the sets.
         '''
         if isinstance(mesh,ngs.comp.Mesh):
             self.ngMesh = mesh.ngmesh
         else:
             self.ngMesh = mesh
         comm = self.comm
+        # Map Netgen string to ids for the labels
+        newLabels = []
+        if len(labels) > 0:
+            for key in labels.keys():
+                for i, s in enumerate(self.ngMesh.GetRegionNames(dim=self.ngMesh.dim-1)):
+                    if key == s:
+                        newLabels = newLabels + [(i+1, labels[key])]
+        newLabels = dict(newLabels)
         if self.ngMesh.dim == 3:
             if comm.rank == 0:
                 V = self.ngMesh.Coordinates()
@@ -179,7 +187,10 @@ def createPETScDMPlex(self, mesh):
                 else:
                     for e in self.ngMesh.Elements2D():
                         join = plex.getFullJoin([vStart+v.nr-1 for v in e.vertices])
-                        plex.setLabelValue(FACE_SETS_LABEL, join[0], int(e.index))
+                        if e.index in newLabels:
+                            plex.setLabelValue(FACE_SETS_LABEL, join[0], int(newLabels[e.index]))
+                        else:
+                            plex.setLabelValue(FACE_SETS_LABEL, join[0], int(e.index))
                     for e in self.ngMesh.Elements1D():
                         join = plex.getJoin([vStart+v.nr-1 for v in e.vertices])
                         plex.setLabelValue(EDGE_SETS_LABEL, join[0], int(e.index))
@@ -200,12 +211,14 @@ def createPETScDMPlex(self, mesh):
                 vStart, _ = plex.getDepthStratum(0)   # vertices
                 for e in self.ngMesh.Elements1D():
                     join = plex.getJoin([vStart+v.nr-1 for v in e.vertices])
-                    plex.setLabelValue(FACE_SETS_LABEL, join[0], int(e.index))
+                    if e.index in newLabels:
+                        plex.setLabelValue(FACE_SETS_LABEL, join[0], int(newLabels[e.index]))
+                    else:
+                        plex.setLabelValue(FACE_SETS_LABEL, join[0], int(e.index))
                 if not (1 == self.ngMesh.Elements2D().NumPy()["index"]).all():
                     for e in self.ngMesh.Elements2D():
                         join = plex.getFullJoin([vStart+v.nr-1 for v in e.vertices])
                         plex.setLabelValue(CELL_SETS_LABEL, join[0], int(e.index))
-
                 self.petscPlex = plex
             else:
                 plex = PETSc.DMPlex().createFromCellList(2,

ngsPETSc/utils/firedrake.py

@@ -216,6 +216,7 @@ def __init__(self, mesh, netgen_flags, user_comm=fd.COMM_WORLD):
         split = flagsUtils(netgen_flags, "split", False)
         quad = flagsUtils(netgen_flags, "quad", False)
         optMoves = flagsUtils(netgen_flags, "optimisation_moves", False)
+        labels = flagsUtils(netgen_flags, "labels", {})
         #Checking the mesh format
         if isinstance(mesh,(ngs.comp.Mesh,ngm.Mesh)):
             if split2tets:
@@ -232,7 +233,7 @@ def __init__(self, mesh, netgen_flags, user_comm=fd.COMM_WORLD):
                 else:
                     raise ValueError("Only 2D and 3D meshes can be optimised.")
             #We create the plex from the netgen mesh
-            self.meshMap = MeshMapping(mesh, comm=self.comm)
+            self.meshMap = MeshMapping(mesh, comm=self.comm, labels=labels)
             #We apply the DMPLEX transform
             if quad:
                 newplex = splitToQuads(self.meshMap.petscPlex, mesh.dim, comm=self.comm)
@@ -364,7 +365,7 @@ def NetgenHierarchy(mesh, levs, flags):
         -tol, geometric tollerance adopted in snapToNetgenDMPlex.
         -refinement_type, the refinment type to be used: uniform (default), Alfeld
     '''
-    if mesh.geometric_dimension() == 3:
+    if mesh.dim == 3:
         raise NotImplementedError("Netgen hierachies are only implemented for 2D meshes.")
     ngmesh = mesh.netgen_mesh
     comm = mesh.comm