New features before Firedrake user meeting

ngsPETSc/plex.py

@@ -33,9 +33,10 @@ class MeshMapping:
 
     '''
 
-    def __init__(self, mesh=None, comm=MPI.COMM_WORLD, name="Default"):
+    def __init__(self, mesh=None, comm=MPI.COMM_WORLD, geo=None, name="Default"):
         self.name = name
         self.comm = comm
+        self.geo = geo
         if isinstance(mesh,(ngs.comp.Mesh,ngm.Mesh)):
             self.createPETScDMPlex(mesh)
         elif isinstance(mesh,PETSc.DMPlex):
@@ -53,6 +54,8 @@ def createNGSMesh(self, plex):
         self.petscPlex = plex
         ngMesh = ngm.Mesh(dim=plex.getCoordinateDim())
         self.ngMesh = ngMesh
+        if self.geo:
+            self.ngMesh.SetGeometry(self.geo)
 
         if plex.getDimension() == 2:
             coordinates = plex.getCoordinates().getArray().reshape([-1,2])
@@ -94,8 +97,12 @@ def createNGSMesh(self, plex):
                 for j in bcIndices:
                     bcIndex = plex.getCone(j)-vStart
                     if len(bcIndex) == 2:
-                        edge = ngm.Element1D([v+1 for v in bcIndex],index=bcLabel)
-                        self.ngMesh.Add(edge)
+                        edge = ngm.Element1D([v+1 for v in bcIndex],
+                                             index=bcLabel,edgenr=bcLabel-1)
+                        if self.geo:
+                            self.ngMesh.Add(edge, project_geominfo=True)
+                        else:
+                            self.ngMesh.Add(edge)
         elif plex.getDimension() == 3:
             coordinates = plex.getCoordinates().getArray().reshape([-1,3])
             self.ngMesh.AddPoints(coordinates)
@@ -139,9 +146,12 @@ def createNGSMesh(self, plex):
                             faces = faces + [face]
                         else:
                             faces = faces + [[face[0],face[2],face[1]]]
-                #fd = self.ngmesh.Add(ngm.FaceDescriptor(bc=bcLabel))
-                self.ngMesh.Add(ngm.FaceDescriptor(bc=bcLabel))
-                self.ngMesh.AddElements(dim=2, index=bcLabel,
+                self.ngMesh.Add(ngm.FaceDescriptor(bc=bcLabel, surfnr=bcLabel))
+                if self.geo:
+                    self.ngMesh.AddElements(dim=2, index=bcLabel,
+                        data=np.asarray(faces,dtype=np.int32), base=0, project_geometry=True)
+                else:
+                    self.ngMesh.AddElements(dim=2, index=bcLabel, edgenr=bcLabel-1,
                                         data=np.asarray(faces,dtype=np.int32), base=0)
         else:
             raise NotImplementedError("No implementation for dimension greater than 3.")
@@ -159,6 +169,7 @@ def createPETScDMPlex(self, mesh):
         else:
             self.ngMesh = mesh
         comm = self.comm
+        self.geo = self.ngMesh.GetGeometry()
         if self.ngMesh.dim == 3:
             if comm.rank == 0:
                 V = self.ngMesh.Coordinates()

ngsPETSc/utils/firedrake/hierarchies.py

@@ -13,37 +13,34 @@
 
 from netgen.meshing import MeshingParameters
 
+from ngsPETSc.plex import MeshMapping
 from ngsPETSc.utils.firedrake.meshes import flagsUtils
 
 def snapToNetgenDMPlex(ngmesh, petscPlex):
     '''
     This function snaps the coordinates of a DMPlex mesh to the coordinates of a Netgen mesh.
     '''
-    if petscPlex.getDimension() == 2:
-        ngCoordinates = ngmesh.Coordinates()
-        petscCoordinates = petscPlex.getCoordinatesLocal().getArray().reshape(-1, ngmesh.dim)
-        for i, pt in enumerate(petscCoordinates):
-            j = np.argmin(np.sum((ngCoordinates - pt)**2, axis=1))
-            petscCoordinates[i] = ngCoordinates[j]
-        petscPlexCoordinates = petscPlex.getCoordinatesLocal()
-        petscPlexCoordinates.setArray(petscPlexCoordinates)
-        petscPlex.setCoordinatesLocal(petscPlexCoordinates)
-    else:
-        raise NotImplementedError("Snapping to Netgen meshes is only implemented for 2D meshes.")
+    ngCoordinates = ngmesh.Coordinates()
+    petscCoordinates = petscPlex.getCoordinatesLocal().getArray().reshape(-1, ngmesh.dim)
+    for i, pt in enumerate(petscCoordinates):
+        j = np.argmin(np.sum((ngCoordinates - pt)**2, axis=1))
+        petscCoordinates[i] = ngCoordinates[j]
+    petscPlexCoordinates = petscPlex.getCoordinatesLocal()
+    petscPlexCoordinates.setArray(petscPlexCoordinates)
+    petscPlex.setCoordinatesLocal(petscPlexCoordinates)
 
 def uniformRefinementRoutine(ngmesh, cdm):
     '''
     Routing called inside of NetgenHierarchy to compute refined ngmesh and plex.
     '''
-    #We refine the netgen mesh uniformly
-    ngmesh.Refine(adaptive=False)
     #We refine the DMPlex mesh uniformly
     cdm.setRefinementUniform(True)
     rdm = cdm.refine()
     rdm.removeLabel("pyop2_core")
     rdm.removeLabel("pyop2_owned")
     rdm.removeLabel("pyop2_ghost")
-    return (rdm, ngmesh)
+    mapping = MeshMapping(rdm, geo=ngmesh.GetGeometry())
+    return (rdm, mapping.ngMesh)
 
 def uniformMapRoutine(meshes):
     '''
@@ -111,8 +108,9 @@ 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:
-        raise NotImplementedError("Netgen hierachies are only implemented for 2D meshes.")
+    if mesh.geometric_dimension() > 3:
+        raise NotImplementedError("Netgen hierachies are only implemented \
+                                  for meshes with dimension greater than 3.")
     ngmesh = mesh.netgen_mesh
     comm = mesh.comm
     #Parsing netgen flags