Mesh Hierarchy support when working with Netgen mesh 2D (#3314)

* Mesh Hierarchy support when working with Netgen mesh

firedrake/mg/mesh.py

@@ -9,7 +9,6 @@
 from firedrake.cython import mgimpl as impl
 from .utils import set_level
 
-
 __all__ = ("HierarchyBase", "MeshHierarchy", "ExtrudedMeshHierarchy", "NonNestedHierarchy",
            "SemiCoarsenedExtrudedHierarchy")
 
@@ -80,28 +79,58 @@ def __getitem__(self, idx):
 
 def MeshHierarchy(mesh, refinement_levels,
                   refinements_per_level=1,
+                  netgen_flags=False,
                   reorder=None,
                   distribution_parameters=None, callbacks=None,
                   mesh_builder=firedrake.Mesh):
     """Build a hierarchy of meshes by uniformly refining a coarse mesh.
 
-    :arg mesh: the coarse :func:`~.Mesh` to refine
-    :arg refinement_levels: the number of levels of refinement
-    :arg refinements_per_level: the number of refinements for each
-        level in the hierarchy.
-    :arg distribution_parameters: options controlling mesh
-        distribution, see :func:`~.Mesh` for details.  If ``None``,
-        use the same distribution parameters as were used to
-        distribute the coarse mesh, otherwise, these options override
-        the default.
-    :arg reorder: optional flag indicating whether to reorder the
-         refined meshes.
-    :arg callbacks: A 2-tuple of callbacks to call before and
+    Parameters
+    ----------
+    mesh : MeshGeometry
+        the coarse mesh to refine
+    refinement_levels : int
+        the number of levels of refinement
+    refinements_per_level : int
+        the number of refinements for each level in the hierarchy.
+    netgen_flags : bool, dict
+        either a bool or a dictionary containing options for Netgen.
+        If not False the hierachy is constructed using ngsPETSc, if
+        None hierarchy constructed in a standard manner.
+    distribution_parameters : dict
+        options controlling mesh distribution, see :py:func:`.Mesh`
+        for details.  If ``None``, use the same distribution
+        parameters as were used to distribute the coarse mesh,
+        otherwise, these options override the default.
+    reorder : bool
+        optional flag indicating whether to reorder the
+        refined meshes.
+    callbacks : tuple
+        A 2-tuple of callbacks to call before and
         after refinement of the DM.  The before callback receives
         the DM to be refined (and the current level), the after
         callback receives the refined DM (and the current level).
-    :arg mesh_builder: Function to turn a DM into a ``Mesh``. Used by pyadjoint.
+    mesh_builder
+        Function to turn a DM into a ``Mesh``. Used by pyadjoint.
+    Returns
+    -------
+    A :py:class:`HierarchyBase` object representing the
+    mesh hierarchy.
     """
+
+    if (isinstance(netgen_flags, bool) and netgen_flags) or isinstance(netgen_flags, dict):
+        try:
+            from ngsPETSc import NetgenHierarchy
+        except ImportError:
+            raise ImportError("Unable to import netgen and ngsPETSc. Please ensure that netgen and ngsPETSc\
+                            are installed and available to Firedrake. You can do this via \
+                            firedrake-update --netgen.")
+        if hasattr(mesh, "netgen_mesh"):
+            return NetgenHierarchy(mesh, refinement_levels, flags=netgen_flags)
+        else:
+            raise RuntimeError("Cannot create a NetgenHierarchy from a mesh that has not been generated by\
+                                Netgen.")
+
     cdm = mesh.topology_dm
     cdm.setRefinementUniform(True)
     dms = []

tests/conftest.py

@@ -20,6 +20,9 @@ def pytest_configure(config):
     config.addinivalue_line(
         "markers",
         "skipplot: mark as skipped if matplotlib is not installed")
+    config.addinivalue_line(
+        "markers",
+        "skipnetgen: mark as skipped if netgen and ngsPETSc is not installed")
 
 
 def pytest_collection_modifyitems(session, config, items):
@@ -39,6 +42,15 @@ def pytest_collection_modifyitems(session, config, items):
     except ImportError:
         ml_backend = False
 
+    try:
+        import netgen
+        del netgen
+        import ngsPETSc
+        del ngsPETSc
+        netgen_installed = True
+    except ImportError:
+        netgen_installed = False
+
     for item in items:
         if complex_mode:
             if item.get_closest_marker("skipcomplex") is not None:
@@ -57,6 +69,10 @@ def pytest_collection_modifyitems(session, config, items):
             if item.get_closest_marker("skipplot") is not None:
                 item.add_marker(pytest.mark.skip(reason="Test cannot be run unless Matplotlib is installed"))
 
+        if not netgen_installed:
+            if item.get_closest_marker("skipnetgen") is not None:
+                item.add_marker(pytest.mark.skip(reason="Test cannot be run unless Netgen and ngsPETSc are installed"))
+
 
 @pytest.fixture(scope="module", autouse=True)
 def check_empty_tape(request):

tests/multigrid/test_netgen_gmg.py

@@ -0,0 +1,101 @@
+from firedrake import *
+import pytest
+
+
+def create_netgen_mesh_circle():
+    from netgen.geom2d import Circle, CSG2d
+    geo = CSG2d()
+
+    circle = Circle(center=(0, 0), radius=1.0, mat="mat1", bc="circle")
+    geo.Add(circle)
+
+    ngmesh = geo.GenerateMesh(maxh=0.75)
+    return ngmesh
+
+
+@pytest.mark.skipcomplex
+@pytest.mark.skipnetgen
+def test_netgen_mg_circle():
+    ngmesh = create_netgen_mesh_circle()
+    mesh = Mesh(ngmesh)
+    nh = MeshHierarchy(mesh, 2, netgen_flags={"degree": 3})
+    mesh = nh[-1]
+
+    V = FunctionSpace(mesh, "CG", 3)
+
+    u = TrialFunction(V)
+    v = TestFunction(V)
+
+    a = inner(grad(u), grad(v))*dx
+    labels = [i+1 for i, name in enumerate(ngmesh.GetRegionNames(codim=1)) if name in ["circle"]]
+    bcs = DirichletBC(V, zero(), labels)
+    x, y = SpatialCoordinate(mesh)
+
+    f = 4+0*x
+    L = f*v*dx
+    exact = (1-x**2-y**2)
+
+    u = Function(V)
+    solve(a == L, u, bcs=bcs, solver_parameters={"ksp_type": "cg",
+                                                 "pc_type": "mg"})
+    expect = Function(V).interpolate(exact)
+    assert (norm(assemble(u - expect)) <= 1e-6)
+
+
+@pytest.mark.skipcomplex
+@pytest.mark.skipnetgeun
+def test_netgen_mg_circle_non_uniform_degree():
+    ngmesh = create_netgen_mesh_circle()
+    mesh = Mesh(ngmesh)
+    nh = MeshHierarchy(mesh, 2, netgen_flags={"degree": [1, 2, 3]})
+    mesh = nh[-1]
+
+    V = FunctionSpace(mesh, "CG", 3)
+
+    u = TrialFunction(V)
+    v = TestFunction(V)
+
+    a = inner(grad(u), grad(v))*dx
+    labels = [i+1 for i, name in enumerate(ngmesh.GetRegionNames(codim=1)) if name in ["circle"]]
+    bcs = DirichletBC(V, zero(), labels)
+    x, y = SpatialCoordinate(mesh)
+
+    f = 4+0*x
+    L = f*v*dx
+    exact = (1-x**2-y**2)
+
+    u = Function(V)
+    solve(a == L, u, bcs=bcs, solver_parameters={"ksp_type": "cg",
+                                                 "pc_type": "mg"})
+    expect = Function(V).interpolate(exact)
+    assert (norm(assemble(u - expect)) <= 1e-6)
+
+
+@pytest.mark.skipcomplex
+@pytest.mark.parallel
+@pytest.mark.skipnetgen
+def test_netgen_mg_circle_parallel():
+    ngmesh = create_netgen_mesh_circle()
+    mesh = Mesh(ngmesh)
+    nh = MeshHierarchy(mesh, 2, netgen_flags={"degree": 3})
+    mesh = nh[-1]
+
+    V = FunctionSpace(mesh, "CG", 3)
+
+    u = TrialFunction(V)
+    v = TestFunction(V)
+
+    a = inner(grad(u), grad(v))*dx
+    labels = [i+1 for i, name in enumerate(ngmesh.GetRegionNames(codim=1)) if name in ["circle"]]
+    bcs = DirichletBC(V, zero(), labels)
+    x, y = SpatialCoordinate(mesh)
+
+    f = 4+0*x
+    L = f*v*dx
+    exact = (1-x**2-y**2)
+
+    u = Function(V)
+    solve(a == L, u, bcs=bcs, solver_parameters={"ksp_type": "cg",
+                                                 "pc_type": "mg"})
+    expect = Function(V).interpolate(exact)
+    assert norm(assemble(u - expect)) <= 1e-6