DG Labelling

Signed-off-by: Umberto Zerbinati <[email protected]>

Author: Umberto Zerbinati

ngsPETSc/plex.py

@@ -183,14 +183,24 @@ def createPETScDMPlex(self, mesh):
                 plex = PETSc.DMPlex().createFromCellList(dim, T, V, comm=comm)
                 plex.setName(self.name)
                 vStart, _ = plex.getDepthStratum(0)
+                nMat = len(self.ngMesh.GetRegionNames(dim=2)) # number of materials
+                nBnd = len(self.ngMesh.GetRegionNames(dim=1)) # number of boundaries
                 if surfMesh:
                     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.surfaces[1]))
+                        if nMat > 1:
+                            cone  = plex.getCone(join[0])
+                            for i in range(len(cone)):
+                                plex.setLabelValue(FACE_SETS_LABEL, cone[i], nBnd+int(e.index))
                 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 nMat > 1:
+                            cone  = plex.getCone(join[0])
+                            for i in range(len(cone)):
+                                plex.setLabelValue(FACE_SETS_LABEL, cone[i], nBnd+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))
@@ -209,9 +219,16 @@ def createPETScDMPlex(self, mesh):
                 plex = PETSc.DMPlex().createFromCellList(2, T, V, comm=comm)
                 plex.setName(self.name)
                 vStart, _ = plex.getDepthStratum(0)   # vertices
+                nMat = len(self.ngMesh.GetRegionNames(dim=2)) # number of materials
+                nBnd = len(self.ngMesh.GetRegionNames(dim=1)) # number of boundaries
                 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 we have more than one material, we set the facet with different labels
+                    if nMat > 1:
+                        cone  = plex.getCone(join[0])
+                        for i in range(len(cone)):
+                            plex.setLabelValue(FACE_SETS_LABEL, cone[i], nBnd+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])

ngsPETSc/utils/firedrake/meshes.py

@@ -255,3 +255,17 @@ def createFromTopology(self, topology, name, comm):
         #Adding refine_marked_elements and curve_field methods
         setattr(fd.MeshGeometry, "refine_marked_elements", refineMarkedElements)
         setattr(fd.MeshGeometry, "curve_field", curveField)
+        #Adding labels for boundary regions and regions
+        Regions = self.meshMap.ngMesh.GetRegionNames(dim=geometric_dim)
+        bndRegions = self.meshMap.ngMesh.GetRegionNames(dim=geometric_dim-1)
+        numberRegions = len(Regions)
+        numberBndRegions = len(bndRegions)
+        self.firedrakeMesh.labels = {}
+        for i in range(numberBndRegions):
+            if bndRegions[i] not in self.firedrakeMesh.labels:
+                self.firedrakeMesh.labels[bndRegions[i]] = []
+            self.firedrakeMesh.labels[bndRegions[i]].append(i+1)
+        for i in range(numberRegions):
+            if Regions[i] not in self.firedrakeMesh.labels:
+                self.firedrakeMesh.labels[Regions[i]] = []
+            self.firedrakeMesh.labels[Regions[i]].append(i+1+numberBndRegions)