Feature/netgen periodic meshes

demos/netgen/netgen_mesh.py.rst

@@ -411,3 +411,99 @@ It is also possible to construct high-order meshes using the ``SplineGeometry``,
 .. figure:: Example7.png
    :align: center
    :alt: Example of a curved mesh of order 2 generated from a geometry described using Netgen CSG2d.
+
+Periodic Meshes
+---------------
+Netgen can identify pairs of vertices lying on opposite boundaries of a geometry as being *the same* point.
+When such a mesh is imported into Firedrake, the identified vertices are merged in the mesh topology, so that
+a continuous (CG) function space automatically shares its degrees of freedom across the seam: the mesh is
+genuinely **periodic**. This is exactly the representation Firedrake uses for its built-in
+``PeriodicRectangleMesh``/``PeriodicBoxMesh``, and it is now available for any Netgen geometry carrying
+periodic identifications.
+
+Identifications are declared on the geometry, before meshing, with the OCC ``Identify`` method::
+
+   shape_a.Identify(shape_b, name, IdentificationType.PERIODIC, transformation)
+
+where ``transformation`` is the rigid motion (typically a translation) that maps ``shape_a`` onto ``shape_b``.
+Netgen then meshes the two boundaries compatibly and records the vertex pairs; Firedrake consumes them
+automatically -- no extra flag on the ``Mesh`` constructor is required.
+
+As a physically motivated example we build the *periodic cylinder*, the classic reduced ("screw pinch") model
+of a tokamak plasma column. A tokamak is a torus, so the plasma is periodic in the toroidal direction; in the
+large-aspect-ratio limit one straightens a toroidal section into a cylinder and identifies its two circular
+ends, recovering periodicity along the axis. We take the axial (toroidal) coordinate to run over :math:`[0, 2\pi)`
+and identify the two end caps by a translation of :math:`2\pi` along ``Z``::
+
+   from netgen.occ import Cylinder, OCCGeometry, Pnt, Z, gp_Trsf, gp_Vec
+   from netgen.meshing import IdentificationType
+   from math import pi as PI
+
+   cyl = Cylinder(Pnt(0, 0, 0), Z, r=1.0, h=2 * PI)
+   # Label the lateral wall, then the two end caps that we will identify.
+   for face in cyl.faces:
+       face.name = "wall"
+   cyl.faces.Min(Z).name = "bottom"
+   cyl.faces.Max(Z).name = "top"
+   # Identify the bottom cap with the top cap: a translation of 2*pi along Z
+   # maps one onto the other, making the axial direction periodic.
+   cyl.faces.Min(Z).Identify(cyl.faces.Max(Z), "toroidal",
+                             IdentificationType.PERIODIC,
+                             gp_Trsf.Translation(gp_Vec(0, 0, 2 * PI)))
+   ngmsh = OCCGeometry(cyl).GenerateMesh(maxh=0.4)
+   msh = Mesh(ngmsh)
+   VTKFile("output/Tokamak.pvd").write(msh)
+
+.. warning::
+
+   The mesh must contain at least a handful of cells along each periodic direction. If a single cell spans a
+   whole period, its two ends are identified and the cell collapses; Firedrake then raises a ``ValueError``
+   asking you to refine along the periodic direction. Here the axis has length :math:`2\pi` and ``maxh=0.4``
+   gives roughly sixteen cells along it, which is ample. Only ``degree == 1`` periodic meshes are supported
+   for now.
+
+Because the two end caps have been identified, no boundary markers survive on them: the seam has become an
+*interior* set of facets, and the only labelled boundary that remains is the lateral wall. This is what makes
+a continuous field wrap around continuously in the axial direction. We can verify the geometry survived the
+merge intact -- the volume of the cylinder is :math:`\pi r^2 h = 2\pi^2` -- while the ends carry no exterior
+facets::
+
+   volume = assemble(Constant(1.0) * dx(domain=msh))
+   PETSc.Sys.Print(f"cylinder volume: {volume:.4f}  (exact 2*pi**2 = {2 * PI**2:.4f})")
+
+To show that the periodicity is doing real work, we solve a Helmholtz problem whose exact solution is periodic
+in the axial coordinate and vanishes on the lateral wall,
+
+.. math::
+
+   u_{\text{ex}}(x, y, z) = \cos(z)\,\bigl(1 - x^2 - y^2\bigr),
+
+so that we can impose a homogeneous Dirichlet condition on the wall while relying on the identified ends for
+continuity along the axis. We look up the id of the ``"wall"`` boundary with ``GetRegionNames`` (as in the
+Poisson example above) and manufacture the right-hand side :math:`f = u_{\text{ex}} - \Delta u_{\text{ex}}` for
+:math:`(I - \Delta)u = f`::
+
+   V = FunctionSpace(msh, "CG", 2)
+   x, y, z = SpatialCoordinate(msh)
+   uex = cos(z) * (1 - x**2 - y**2)
+   f = uex - div(grad(uex))
+
+   u = TrialFunction(V)
+   v = TestFunction(V)
+   a = (inner(u, v) + inner(grad(u), grad(v))) * dx
+   L = inner(f, v) * dx
+
+   labels = [i + 1 for i, name in enumerate(ngmsh.GetRegionNames(codim=1)) if name == "wall"]
+   bc = DirichletBC(V, 0, labels)
+
+   sol = Function(V)
+   solve(a == L, sol, bcs=bc)
+   VTKFile("output/TokamakSolution.pvd").write(sol)
+
+   error = sqrt(assemble(inner(sol - uex, sol - uex) * dx))
+   PETSc.Sys.Print(f"L2 error: {error:.2e}")
+
+The recovered solution is continuous across the identified ends: opening ``output/TokamakSolution.pvd`` in
+ParaView, the field wraps seamlessly from the top cap back to the bottom, exactly as a toroidal mode should.
+Had the ends *not* been identified, the same computation would leave an artificial jump at the seam and the
+manufactured solution would not be recovered.

firedrake/mesh.py

@@ -3055,6 +3055,80 @@ def curve_field(self, order, permutation_tol=1e-8, cg_field=None):
         new_coordinates.dat.data_wo_with_halos[broken_indices] = own_curved_points
         return new_coordinates
 
+    @PETSc.Log.EventDecorator()
+    def _periodic_coordinates(self, permutation_tol=1e-8):
+        '''Return a discontinuous coordinate field for a periodic netgen mesh.
+
+        A periodic netgen mesh is converted by ngsPETSc into a DMPlex whose
+        topology is periodic (the identified vertices are merged) but whose
+        continuous coordinates are "wrapped" at the periodic seam. This method
+        builds the discontinuous (DG1) coordinate field carrying each cell's
+        true, un-wrapped corner coordinates, which is then attached to the mesh
+        (see :func:`~.utility_meshes._postprocess_periodic_mesh`).
+
+        This method requires that the mesh has been constructed from a netgen
+        mesh that carries periodic identifications.
+
+        :arg permutation_tol: tolerance used to match the reference element nodes.
+        '''
+        utils.check_netgen_installed()
+        from firedrake.netgen import find_permutation
+        from firedrake.function import Function
+        from firedrake.functionspace import VectorFunctionSpace
+        from ngsPETSc.plex import buildPeriodicVertexMap
+        from ngsPETSc.utils.utils import trim_util
+
+        if not hasattr(self, "netgen_mesh"):
+            raise ValueError("Cannot build periodic coordinates for a mesh that "
+                             "has not been generated by netgen.")
+
+        # Netgen element -> vertex connectivity (0-based).
+        if self.topological_dimension == 2:
+            ng_element = self.netgen_mesh.Elements2D()
+        else:
+            ng_element = self.netgen_mesh.Elements3D()
+        conn = trim_util(ng_element.NumPy()["nodes"])
+
+        # The same vertex merging ngsPETSc applied when it built the periodic plex.
+        old_to_new, survivors, _ = buildPeriodicVertexMap(self.netgen_mesh)
+        coords = self.netgen_mesh.Coordinates()
+        # `unwrapped` is the true geometry (what we want to store); `wrapped` is the
+        # merged/representative geometry, which coincides with Firedrake's continuous
+        # coordinates on the periodic plex.
+        unwrapped = coords[conn]
+        wrapped = coords[survivors][old_to_new[conn]]
+
+        # Index netgen cells by their (rounded) set of wrapped vertex coordinates,
+        # so each Firedrake cell can be matched to its netgen element by geometry
+        # alone. This is robust to mesh reordering and parallel redistribution.
+        def cell_key(pts):
+            return tuple(sorted(tuple(np.round(p, 8)) for p in pts))
+        lookup = {cell_key(wrapped[e]): e for e in range(wrapped.shape[0])}
+
+        # Build the DG1 coordinate field, initialised from the (wrapped) continuous
+        # coordinates, then overwrite every cell with its un-wrapped geometry. The
+        # equispaced DG element matches the layout expected by _set_dg_coordinates
+        # (as used by Firedrake's own periodic meshes).
+        broken_space = VectorFunctionSpace(
+            self, finat.ufl.FiniteElement("DG", self.ufl_cell(), 1, variant="equispaced")
+        )
+        new_coordinates = Function(broken_space).interpolate(self.coordinates)
+        cell_nodes = new_coordinates.cell_node_map().values
+        data = new_coordinates.dat.data
+        # Snapshot the wrapped coordinates (in Firedrake node order) before
+        # overwriting, so the per-cell matching always sees the interpolated values.
+        wrapped_fd = data[cell_nodes].real.copy()
+
+        for i in range(cell_nodes.shape[0]):
+            fd_nodes = wrapped_fd[i]
+            e = lookup[cell_key(fd_nodes)]
+            # permutation taking the netgen node order to this cell's node order
+            permutation = find_permutation(
+                wrapped[e][np.newaxis], fd_nodes[np.newaxis], tol=permutation_tol
+            )[0]
+            data[cell_nodes[i]] = unwrapped[e][permutation]
+        return new_coordinates
+
 
 @PETSc.Log.EventDecorator()
 def make_mesh_from_coordinates(coordinates, name, tolerance=0.5):
@@ -3303,6 +3377,13 @@ def Mesh(meshfile, **kwargs):
 
     :param netgen_flags: The dictionary of flags to be passed to ngsPETSc.
 
+    If the Netgen mesh carries periodic identifications (e.g. created with
+    ``shape.Identify(..., IdentificationType.PERIODIC, ...)``) the resulting
+    Firedrake mesh is periodic: the identified vertices are merged and a
+    discontinuous coordinate field carries the un-wrapped geometry. The mesh
+    must be fine enough that no cell spans a full period, and high-order curving
+    (``degree != 1``) of periodic meshes is not currently supported.
+
     When the mesh is read from a file the following mesh formats
     are supported (determined, case insensitively, from the
     filename extension):
@@ -3360,6 +3441,7 @@ def Mesh(meshfile, **kwargs):
     # they all immediately call a petsc4py which in turn uses a PETSc
     # internal comm
     geometric_dim = kwargs.get("dim", None)
+    netgen_periodic = False
     if isinstance(meshfile, PETSc.DMPlex):
         plex = meshfile
         if MPI.Comm.Compare(user_comm, plex.comm.tompi4py()) not in {MPI.CONGRUENT, MPI.IDENT}:
@@ -3372,6 +3454,11 @@ def Mesh(meshfile, **kwargs):
         netgen_firedrake_mesh = FiredrakeMesh(meshfile, netgen_flags, user_comm)
         plex = netgen_firedrake_mesh.meshMap.petscPlex
         plex.setName(_generate_default_mesh_topology_name(name))
+        # A periodic netgen mesh produces a vertex-merged (periodic) topology that
+        # is finished off with a discontinuous coordinate field below. That field
+        # is built on the un-reordered topology (as for Firedrake's own periodic
+        # meshes), so suppress reordering here and reapply it in postprocessing.
+        netgen_periodic = len(netgen_firedrake_mesh.meshMap.ngMesh.GetIdentifications()) > 0
 
     else:
         basename, ext = os.path.splitext(meshfile)
@@ -3395,7 +3482,11 @@ def Mesh(meshfile, **kwargs):
         plex.setName(_generate_default_mesh_topology_name(name))
     # Create mesh topology
     submesh_parent = kwargs.get("submesh_parent", None)
-    topology = MeshTopology(plex, name=plex.getName(), reorder=reorder,
+    # A periodic netgen mesh is finished off with a discontinuous coordinate field
+    # built on the un-reordered topology; the requested reordering is reapplied in
+    # _postprocess_periodic_mesh.
+    topology_reorder = False if netgen_periodic else reorder
+    topology = MeshTopology(plex, name=plex.getName(), reorder=topology_reorder,
                             distribution_parameters=distribution_parameters,
                             distribution_name=kwargs.get("distribution_name"),
                             permutation_name=kwargs.get("permutation_name"),
@@ -3407,9 +3498,30 @@ def Mesh(meshfile, **kwargs):
         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:
+        periodic = len(mesh.netgen_mesh.GetIdentifications()) > 0
+        if periodic:
+            # ngsPETSc produced a periodic (vertex-merged) topology; attach the
+            # discontinuous coordinate field carrying the un-wrapped geometry.
+            if degree != 1:
+                raise NotImplementedError(
+                    "High-order curving of periodic netgen meshes is not supported yet."
+                )
+            from firedrake.utility_meshes import _postprocess_periodic_mesh
+            permutation_tol = netgen_flags.get("permutation_tol", 1e-8)
+            coordinates = mesh._periodic_coordinates(permutation_tol=permutation_tol)
+            temp = _postprocess_periodic_mesh(coordinates,
+                                              mesh.comm,
+                                              distribution_parameters,
+                                              reorder,
+                                              name,
+                                              kwargs.get("distribution_name"),
+                                              kwargs.get("permutation_name"))
+            temp.netgen_mesh = mesh.netgen_mesh
+            temp.netgen_flags = mesh.netgen_flags
+            mesh = temp
+        # Curve the mesh, if requested
+        elif degree != 1:
             permutation_tol = netgen_flags.get("permutation_tol", 1e-8)
             cg = netgen_flags.get("cg", None)
             coordinates = mesh.curve_field(

tests/firedrake/regression/test_netgen.py

@@ -342,3 +342,92 @@
             break
         mesh = adapt(mesh, eta)
     assert error_estimators[-1] < 0.06
+
+
+def _occ_periodic_square(maxh, directions="x"):
+    from netgen.occ import Rectangle, OCCGeometry, X, Y, gp_Trsf, gp_Vec
+    from netgen.meshing import IdentificationType
+    shape = Rectangle(1, 1).Face()
+    shape.edges.Min(X).name, shape.edges.Max(X).name = "left", "right"
+    shape.edges.Min(Y).name, shape.edges.Max(Y).name = "bottom", "top"
+    if "x" in directions:
+        shape.edges.Min(X).Identify(shape.edges.Max(X), "px",
+                                    IdentificationType.PERIODIC,
+                                    gp_Trsf.Translation(gp_Vec(1, 0, 0)))
+    if "y" in directions:
+        shape.edges.Min(Y).Identify(shape.edges.Max(Y), "py",
+                                    IdentificationType.PERIODIC,
+                                    gp_Trsf.Translation(gp_Vec(0, 1, 0)))
+    return OCCGeometry(shape, dim=2).GenerateMesh(maxh=maxh)
+
+
+def _solve_periodic_helmholtz(mesh, uex):
+    # Solve (I - div grad) u = f, where f is manufactured so the exact solution
+    # is the periodic field uex; return the L2 error.
+    V = FunctionSpace(mesh, "CG", 1)
+    u, v = TrialFunction(V), TestFunction(V)
+    f = uex - div(grad(uex))
+    a = (inner(u, v) + inner(grad(u), grad(v))) * dx
+    L = inner(f, v) * dx
+    uh = Function(V)
+    solve(a == L, uh)
+    return sqrt(assemble(inner(uh - uex, uh - uex) * dx))
+
+
+@pytest.mark.skipnetgen
+def test_netgen_periodic_square():
+    # A periodic netgen mesh identifies the seam DOFs, so the periodic boundary
+    # markers are absent and the geometry is recovered exactly.
+    mesh = Mesh(_occ_periodic_square(0.1, directions="x"))
+    assert abs(assemble(Constant(1.0) * dx(domain=mesh)) - 1.0) < 1e-10
+
+    x, y = SpatialCoordinate(mesh)
+    # exact solution periodic in x; the seam continuity is essential to recover it.
+    err = _solve_periodic_helmholtz(mesh, sin(2 * pi * x))
+    assert err < 5e-2
+
+
+@pytest.mark.skipnetgen
+def test_netgen_periodic_square_both_directions():
+    mesh = Mesh(_occ_periodic_square(0.1, directions="xy"))
+    assert abs(assemble(Constant(1.0) * dx(domain=mesh)) - 1.0) < 1e-10
+    # both pairs identified: no exterior facets remain (a torus).
+    assert mesh.exterior_facets.set.total_size == 0
+    x, y = SpatialCoordinate(mesh)
+    err = _solve_periodic_helmholtz(mesh, sin(2 * pi * x) * cos(2 * pi * y))
+    assert err < 5e-2
+
+
+@pytest.mark.skipnetgen
+def test_netgen_periodic_cylinder():
+    # Periodic along the axis of a cylinder: a curved-boundary periodic mesh.
+    from netgen.occ import Cylinder, OCCGeometry, Pnt, gp_Trsf, gp_Vec, Z
+    from netgen.meshing import IdentificationType
+    cyl = Cylinder(Pnt(0, 0, 0), Z, r=1.0, h=1.0)
+    cyl.faces.Min(Z).Identify(cyl.faces.Max(Z), "pz",
+                              IdentificationType.PERIODIC,
+                              gp_Trsf.Translation(gp_Vec(0, 0, 1)))
+    mesh = Mesh(OCCGeometry(cyl).GenerateMesh(maxh=0.15))
+    # Volume approaches pi as the polygonal boundary is refined.
+    assert abs(assemble(Constant(1.0) * dx(domain=mesh)) - pi) < 2e-2
+    x, y, z = SpatialCoordinate(mesh)
+    err = _solve_periodic_helmholtz(mesh, sin(2 * pi * z))
+    assert err < 2e-1
+
+
+@pytest.mark.skipnetgen
+def test_netgen_periodic_too_coarse():
+    # A mesh too coarse along the periodic direction produces seam-spanning cells
+    # that collapse on merging; this must raise a clear error rather than build a
+    # broken mesh.
+    from netgen.occ import Box, OCCGeometry, X, gp_Trsf, gp_Vec, Pnt
+    from netgen.meshing import IdentificationType
+    box = Box(Pnt(0, 0, 0), Pnt(1, 1, 1))
+    box.faces.Min(X).Identify(box.faces.Max(X), "px",
+                              IdentificationType.PERIODIC,
+                              gp_Trsf.Translation(gp_Vec(1, 0, 0)))
+    ngmesh = OCCGeometry(box).GenerateMesh(maxh=0.4)
+    with pytest.raises(ValueError, match="degenerate"):
+        Mesh(ngmesh)
+
+