Support for 3D hierarchy, periodic hierarchy and added digits option.

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

ngsPETSc/utils/firedrake.py

@@ -77,7 +77,7 @@ def refineMarkedElements(self, mark):
     else:
         raise NotImplementedError("No implementation for dimension other than 2 and 3.")
 
-def curveField(self, order):
+def curveField(self, order, digits=8):
     '''
     This method returns a curved mesh as a Firedrake function.
 
@@ -104,7 +104,7 @@ def curveField(self, order):
     V = newFunctionCoordinates.dat.data
     getIdx = self._cell_numbering.getOffset
     refPts = np.array(refPts)
-    rnd = lambda x: round(x, 8)
+    rnd = lambda x: round(x, digits)
     if self.geometric_dimension() == 2:
         #Mapping to the physical domain
         physPts = np.ndarray((len(self.netgen_mesh.Elements2D()),
@@ -219,7 +219,7 @@ def createFromTopology(self, topology, name):
         setattr(fd.MeshGeometry, "refine_marked_elements", refineMarkedElements)
         setattr(fd.MeshGeometry, "curve_field", curveField)
 
-def NetgenHierarchy(ngmesh, levs, order=1, comm=fd.COMM_WORLD):
+def NetgenHierarchy(ngmesh, levs, order=1, digits=8, adaptive=False, comm=fd.COMM_WORLD):
     '''
     This function creates a Firedrake mesh hierarchy from Netgen/NGSolve meshes.
 
@@ -235,43 +235,37 @@ def NetgenHierarchy(ngmesh, levs, order=1, comm=fd.COMM_WORLD):
     if mesh.comm.size > 1 and mesh._grown_halos:
         raise RuntimeError("Cannot refine parallel overlapped meshes ")
     #Computing useful quantities
-    fd.FunctionSpace(mesh, "DG", 0)
-    coarse_index = mesh._cell_numbering.getOffset
-    if ngmesh.dim == 2:
-        number_coarse_cells = len(ngmesh.Elements2D())
-    elif ngmesh.dim == 3:
-        pass
+    ngElements = {2: ngmesh.Elements2D, 3: ngmesh.Elements3D}
+    number_coarse_cells = len(ngElements[ngmesh.dim]())
     #Curving linear mesh
-    mesh = fd.Mesh(mesh.curve_field(order=order))
+    mesh = fd.Mesh(mesh.curve_field(order=order, digits=digits))
     meshes += [mesh]
     for lv in range(levs):
+        #Streightening the mesh
+        ngmesh.Curve(1)
         #We refine the netgen mesh uniformly
-        ngmesh.Refine(adaptive=False)
+        ngmesh.Refine(adaptive=adaptive)
         #We construct the refined linear mesh and curve it
         mesh = fd.Mesh(ngmesh, comm=comm)
-        if ngmesh.dim == 2:
-            number_fine_cells = len(ngmesh.Elements2D())
-        elif ngmesh.dim == 3:
-            pass
-        mesh = fd.Mesh(mesh.curve_field(order=order))
+        number_fine_cells = len(ngElements[ngmesh.dim]())
+        mesh = fd.Mesh(mesh.curve_field(order=order, digits=digits))
         meshes += [mesh]
         #We populate the coarse to fine map
         fine_cells_per_coarse_cell = number_fine_cells // number_coarse_cells
         c2f = np.zeros((number_coarse_cells,fine_cells_per_coarse_cell),dtype=np.int32)
         f2c = np.zeros((number_fine_cells,1),dtype=np.int32)
         coarse_counter = [0]*number_coarse_cells
-        if ngmesh.dim == 2:
-            V = ngmesh.Coordinates()
-            T = ngmesh.Elements2D().NumPy()["nodes"]
-            T = np.array([list(np.trim_zeros(a, 'b')) for a in list(T)])-1
-            for i in range(number_fine_cells):
-                coarse_index = meshes[-2].locate_cell(sum(V[T[i]])/3)
-                fine_index = meshes[-1].locate_cell(sum(V[T[i]])/3)
-                c2f[coarse_index][coarse_counter[coarse_index]] = fine_index
-                coarse_counter[coarse_index] += 1
-                f2c[fine_index][0] = coarse_index
-            coarse_to_fine_cells.append(c2f)
-            fine_to_coarse_cells.append(f2c)
+        V = ngmesh.Coordinates()
+        T = ngElements[ngmesh.dim]().NumPy()["nodes"]
+        T = np.array([list(np.trim_zeros(a, 'b')) for a in list(T)])-1
+        for i in range(number_fine_cells):
+            coarse_index = meshes[-2].locate_cell(sum(V[T[i]])/len(T[i]))
+            fine_index = meshes[-1].locate_cell(sum(V[T[i]])/len(T[i]))
+            c2f[coarse_index][coarse_counter[coarse_index]] = fine_index
+            coarse_counter[coarse_index] += 1
+            f2c[fine_index][0] = coarse_index
+        coarse_to_fine_cells.append(c2f)
+        fine_to_coarse_cells.append(f2c)
         number_coarse_cells = number_fine_cells
 
     coarse_to_fine_cells = dict((Fraction(i, refinements_per_level), c2f)