Netgen: curved meshes from Mesh constructor and refine_marked_elements

firedrake/mesh.py

@@ -2994,9 +2994,10 @@ def refine_marked_elements(self, mark, netgen_flags=None):
         utils.check_netgen_installed()
 
         if netgen_flags is None:
-            netgen_flags = {}
+            netgen_flags = self.netgen_flags
         DistParams = self._distribution_parameters
-        els = {2: self.netgen_mesh.Elements2D, 3: self.netgen_mesh.Elements3D}
+        netgen_mesh = self.netgen_mesh.Copy()
+        els = {2: netgen_mesh.Elements2D, 3: netgen_mesh.Elements3D}
         dim = self.geometric_dimension
         refine_faces = netgen_flags.get("refine_faces", False)
         if dim in [2, 3]:
@@ -3019,12 +3020,14 @@ def refine_marked_elements(self, mark, netgen_flags=None):
                             else:
                                 el.refine = False
                         if not refine_faces and dim == 3:
-                            self.netgen_mesh.Elements2D().NumPy()["refine"] = 0
-                        self.netgen_mesh.Refine(adaptive=True)
+                            netgen_mesh.Elements2D().NumPy()["refine"] = 0
+                        netgen_mesh.Refine(adaptive=True)
                         mark = mark-np.ones(mark.shape)
-                    return fd.Mesh(self.netgen_mesh, distribution_parameters=DistParams, comm=self.comm)
+                    return fd.Mesh(netgen_mesh, distribution_parameters=DistParams,
+                                   netgen_flags=netgen_flags, comm=self.comm)
                 return fd.Mesh(netgen.libngpy._meshing.Mesh(dim),
-                               distribution_parameters=DistParams, comm=self.comm)
+                               distribution_parameters=DistParams,
+                               netgen_flags=netgen_flags, comm=self.comm)
         else:
             raise NotImplementedError("No implementation for dimension other than 2 and 3.")
 
@@ -3431,6 +3434,24 @@ def Mesh(meshfile, **kwargs):
 
     if from_netgen:
         mesh.netgen_mesh = netgen_firedrake_mesh.meshMap.ngMesh
+        mesh.netgen_flags = netgen_flags
+
+        # Curve the mesh, if requested
+        degree = netgen_flags.get("degree", 1)
+        if degree != 1:
+            permutation_tol = netgen_flags.get("permutation_tol", 1e-8)
+            cg = netgen_flags.get("cg", False)
+            ho_field = mesh.curve_field(
+                order=degree,
+                permutation_tol=permutation_tol,
+                cg_field=cg
+            )
+            temp = Mesh(ho_field, distribution_parameters=mesh._distribution_parameters,
+                        comm=mesh.comm)
+            temp.netgen_mesh = mesh.netgen_mesh
+            temp.netgen_flags = netgen_flags
+            temp._tolerance = mesh.tolerance
+            mesh = temp
 
     mesh.submesh_parent = submesh_parent
     mesh._tolerance = tolerance

firedrake/mg/netgen.py

@@ -246,8 +246,6 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None):
     logger.info(f"\tOrder of the hierarchy: {order}")
     if isinstance(order, int):
         order = [order]*(levs+1)
-    permutation_tol = flags.get("permutation_tol", 1e-8)
-    location_tol = flags.get("location_tol", 1e-8)
     refType = flags.get("refinement_type", "uniform")
     logger.info(f"\tRefinement type: {refType}")
     optMoves = flags.get("optimisation_moves", False)
@@ -266,15 +264,11 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None):
         parameters.update(mesh._distribution_parameters)
     parameters["partition"] = False
     # Curve the mesh
-    if order[0] > 1:
-        ho_field = mesh.curve_field(
-            order=order[0],
-            permutation_tol=permutation_tol,
-            cg_field=cg
-        )
-        temp = fd.Mesh(ho_field, distribution_parameters=parameters, comm=comm)
-        temp.netgen_mesh = mesh.netgen_mesh
-        temp._tolerance = mesh.tolerance
+    if order[0] != mesh.coordinates.function_space().ufl_element().degree():
+        temp_flags = dict(flags)
+        temp_flags['degree'] = order[0]
+        temp = fd.Mesh(mesh.netgen_mesh, distribution_parameters=parameters,
+                       netgen_flags=temp_flags, comm=comm)
         mesh = temp
     # Make a plex (cdm) without overlap.
     dm_cell_type, = mesh.dm_cell_types
@@ -313,17 +307,19 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None):
                 raise ValueError("Only 2D and 3D meshes can be optimised.")
         mesh.netgen_mesh = ngmesh
         # Curve the mesh
-        if order[l+1] > 1:
+        if order[l+1] != mesh.coordinates.function_space().ufl_element().degree():
             logger.info("\t\t\tCurving the mesh ...")
             tic = time.time()
             if snap == "geometry":
+                temp_flags = dict(flags)
+                temp_flags['degree'] = order[l+1]
                 mesh = fd.Mesh(
-                    mesh.curve_field(order=order[l+1],
-                                     location_tol=location_tol,
-                                     permutation_tol=permutation_tol),
+                    mesh.netgen_mesh,
                     distribution_parameters=parameters,
+                    netgen_flags=temp_flags,
                     comm=comm)
             elif snap == "coarse":
+                ho_field = meshes[0].coordinates
                 mesh = snapToCoarse(ho_field, mesh, order[l+1], snap_smoothing, cg)
             toc = time.time()
             logger.info(f"\t\t\tMeshed curved. Time taken: {toc-tic}")