Reorganise how netgen meshes are curved

Author: pefarrell

firedrake/mesh.py

@@ -2994,7 +2994,7 @@ 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
         netgen_mesh = self.netgen_mesh.Copy()
         els = {2: netgen_mesh.Elements2D, 3: netgen_mesh.Elements3D}
@@ -3023,9 +3023,11 @@ def refine_marked_elements(self, mark, netgen_flags=None):
                             netgen_mesh.Elements2D().NumPy()["refine"] = 0
                         netgen_mesh.Refine(adaptive=True)
                         mark = mark-np.ones(mark.shape)
-                    return fd.Mesh(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.")
 
@@ -3432,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

@@ -266,15 +266,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,15 +309,16 @@ 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":
                 mesh = snapToCoarse(ho_field, mesh, order[l+1], snap_smoothing, cg)