splits option in Netgen falgs

Signed-off-by: Umberto Zerbinati <zerbinati@maths.ox.ac.uk>

ngsPETSc/utils/firedrake.py

@@ -29,6 +29,15 @@ class comp:
 
 from ngsPETSc import MeshMapping
 
+def flagsUtils(flags, option, default):
+    '''
+    utility fuction used to parse Netgen flag options
+    '''
+    try:
+        return flags[option]
+    except KeyError:
+        return default
+
 def refineMarkedElements(self, mark):
     '''
     This method is used to refine a mesh based on a marking function
@@ -165,6 +174,9 @@ def curveField(self, order, tol=1e-8):
                             newFunctionCoordinates.sub(2).dat.data[datIdx] = curvedPhysPts[i][j][2]
     return newFunctionCoordinates
 
+splitTypes = {"Alfeld": lambda x: x.SplitAlfeld(),
+              "Powell-Sabin": lambda x: x.SplitPowellSabin()}
+
 class FiredrakeMesh:
     '''
     This class creates a Firedrake mesh from Netgen/NGSolve meshes.
@@ -175,15 +187,21 @@ class FiredrakeMesh:
     '''
     def __init__(self, mesh, netgen_flags, user_comm=fd.COMM_WORLD):
         self.comm = user_comm
+        #Parsing netgen flags
+        if not isinstance(netgen_flags, dict):
+            netgen_flags = {}
+        split2tets = flagsUtils(netgen_flags, "split_to_tets", False)
+        split = flagsUtils(netgen_flags, "split", False)
+        #Checking the mesh format
         if isinstance(mesh,(ngs.comp.Mesh,ngm.Mesh)):
-            try:
-                if netgen_flags["purify_to_tets"]:
-                    mesh.Split2Tets()
-            except KeyError:
-                warnings.warn("No purify_to_tets flag found, mesh will not be purified to tets.")
+            if split2tets:
+                mesh = mesh.Split2Tets()
+            if split:
+                splitTypes[split](mesh)
             self.meshMap = MeshMapping(mesh, comm=self.comm)
         else:
             raise ValueError("Mesh format not recognised.")
+        #We quadrilateralise the mesh or apply a transform if requested
         try:
             if netgen_flags["quad"]:
                 transform = PETSc.DMPlexTransform().create(comm=self.comm)
@@ -248,15 +266,6 @@ def snapToNetgenDMPlex(ngmesh, petscPlex):
     else:
         raise NotImplementedError("Snapping to Netgen meshes is only implemented for 2D meshes.")
 
-def flagsUtils(flags, option, default):
-    '''
-    utility fuction used to parse Netgen flag options
-    '''
-    try:
-        return flags[option]
-    except KeyError:
-        return default
-
 def uniformRefinementRoutine(ngmesh, cdm):
     '''
     Routing called inside of NetgenHierarchy to compute refined ngmesh and plex.
@@ -271,6 +280,33 @@ def uniformRefinementRoutine(ngmesh, cdm):
     rdm.removeLabel("pyop2_ghost")
     return (rdm, ngmesh)
 
+def uniformMapRoutine(meshes):
+    '''
+    This function computes the coarse to fine and fine to coarse maps
+    for a uniform mesh hierarchy.
+    '''
+    refinements_per_level = 1
+    lgmaps = []
+    for i, m in enumerate(meshes):
+        no = impl.create_lgmap(m.topology_dm)
+        m.init()
+        o = impl.create_lgmap(m.topology_dm)
+        m.topology_dm.setRefineLevel(i)
+        lgmaps.append((no, o))
+    coarse_to_fine_cells = []
+    fine_to_coarse_cells = [None]
+    for (coarse, fine), (clgmaps, flgmaps) in zip(zip(meshes[:-1], meshes[1:]),
+                                                zip(lgmaps[:-1], lgmaps[1:])):
+        c2f, f2c = impl.coarse_to_fine_cells(coarse, fine, clgmaps, flgmaps)
+        coarse_to_fine_cells.append(c2f)
+        fine_to_coarse_cells.append(f2c)
+
+    coarse_to_fine_cells = dict((Fraction(i, refinements_per_level), c2f)
+                                for i, c2f in enumerate(coarse_to_fine_cells))
+    fine_to_coarse_cells = dict((Fraction(i, refinements_per_level), f2c)
+                                for i, f2c in enumerate(fine_to_coarse_cells))
+    return (coarse_to_fine_cells, fine_to_coarse_cells)
+
 def alfeldRefinementRoutine(ngmesh, cdm):
     '''
     Routing called inside of NetgenHierarchy to compute refined ngmesh and plex.
@@ -285,8 +321,15 @@ def alfeldRefinementRoutine(ngmesh, cdm):
     rdm = tr.apply(cdm)
     return (rdm, ngmesh)
 
-refinementTypes = {"uniform": uniformRefinementRoutine,
-                   "Alfeld": alfeldRefinementRoutine}
+def alfeldMapRoutine(meshes):
+    '''
+    This function computes the coarse to fine and fine to coarse maps
+    for a alfeld mesh hierarchy.
+    '''
+    raise NotImplementedError("Alfeld refinement is not implemented yet.")
+
+refinementTypes = {"uniform": (uniformRefinementRoutine, uniformMapRoutine),
+                   "Alfeld": (alfeldRefinementRoutine, alfeldMapRoutine)}
 
 def NetgenHierarchy(mesh, levs, flags):
     '''
@@ -317,7 +360,6 @@ def NetgenHierarchy(mesh, levs, flags):
     refType = flagsUtils(flags, "refinement_type", "uniform")
     #Firedrake quoantities
     meshes = []
-    refinements_per_level = 1
     coarse_to_fine_cells = []
     fine_to_coarse_cells = [None]
     params = {"partition": False}
@@ -333,7 +375,7 @@ def NetgenHierarchy(mesh, levs, flags):
     for l in range(levs):
         #Streightening the mesh
         ngmesh.Curve(1)
-        rdm, ngmesh = refinementTypes[refType](ngmesh, cdm)
+        rdm, ngmesh = refinementTypes[refType][0](ngmesh, cdm)
         cdm = rdm
         #We snap the mesh to the Netgen mesh
         snapToNetgenDMPlex(ngmesh, rdm)
@@ -347,24 +389,6 @@ def NetgenHierarchy(mesh, levs, flags):
                            distribution_parameters=params, comm=comm)
         meshes += [mesh]
     #We populate the coarse to fine map
-    lgmaps = []
-    for i, m in enumerate(meshes):
-        no = impl.create_lgmap(m.topology_dm)
-        m.init()
-        o = impl.create_lgmap(m.topology_dm)
-        m.topology_dm.setRefineLevel(i)
-        lgmaps.append((no, o))
-    coarse_to_fine_cells = []
-    fine_to_coarse_cells = [None]
-    for (coarse, fine), (clgmaps, flgmaps) in zip(zip(meshes[:-1], meshes[1:]),
-                                                zip(lgmaps[:-1], lgmaps[1:])):
-        c2f, f2c = impl.coarse_to_fine_cells(coarse, fine, clgmaps, flgmaps)
-        coarse_to_fine_cells.append(c2f)
-        fine_to_coarse_cells.append(f2c)
-
-    coarse_to_fine_cells = dict((Fraction(i, refinements_per_level), c2f)
-                                for i, c2f in enumerate(coarse_to_fine_cells))
-    fine_to_coarse_cells = dict((Fraction(i, refinements_per_level), f2c)
-                                for i, f2c in enumerate(fine_to_coarse_cells))
+    coarse_to_fine_cells, fine_to_coarse_cells = refinementTypes[refType][1](meshes)
     return fd.HierarchyBase(meshes, coarse_to_fine_cells, fine_to_coarse_cells,
-                         refinements_per_level, nested=False)
+                            1, nested=False)