diff --git a/demos/netgen/netgen_mesh.py.rst b/demos/netgen/netgen_mesh.py.rst index 6d0f488c83..88b2c1fafa 100755 --- a/demos/netgen/netgen_mesh.py.rst +++ b/demos/netgen/netgen_mesh.py.rst @@ -411,3 +411,102 @@ 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. + diff --git a/demos/periodic_meshes/Example1.png b/demos/periodic_meshes/Example1.png new file mode 100644 index 0000000000..2f6001ea4c Binary files /dev/null and b/demos/periodic_meshes/Example1.png differ diff --git a/demos/periodic_meshes/Example2.png b/demos/periodic_meshes/Example2.png new file mode 100644 index 0000000000..56f7386584 Binary files /dev/null and b/demos/periodic_meshes/Example2.png differ diff --git a/demos/periodic_meshes/Example3.png b/demos/periodic_meshes/Example3.png new file mode 100644 index 0000000000..218de4213c Binary files /dev/null and b/demos/periodic_meshes/Example3.png differ diff --git a/demos/periodic_meshes/Example4.png b/demos/periodic_meshes/Example4.png new file mode 100644 index 0000000000..28f40a971e Binary files /dev/null and b/demos/periodic_meshes/Example4.png differ diff --git a/demos/periodic_meshes/Example5.png b/demos/periodic_meshes/Example5.png new file mode 100644 index 0000000000..53504dd44f Binary files /dev/null and b/demos/periodic_meshes/Example5.png differ diff --git a/demos/periodic_meshes/Example6.png b/demos/periodic_meshes/Example6.png new file mode 100644 index 0000000000..9280cac5be Binary files /dev/null and b/demos/periodic_meshes/Example6.png differ diff --git a/demos/periodic_meshes/periodic_meshes.py.rst b/demos/periodic_meshes/periodic_meshes.py.rst new file mode 100644 index 0000000000..7aa9e2a350 --- /dev/null +++ b/demos/periodic_meshes/periodic_meshes.py.rst @@ -0,0 +1,449 @@ +Periodic meshes in Firedrake +============================ + +This tutorial was contributed by `Thomas Higham `__ and `Umberto Zerbinati `__. + +The purpose of this demo is to summarise the support for periodic meshes in Firedrake. +Firedrake can build periodic meshes for simple one- and two-dimensional geometries and can also periodically extrude two-dimensional meshes into three-dimensional prism meshes. +Firedrake also has support for periodic tetrahedral meshes generated in Netgen. + +We begin by importing the necessary libraries: :: + + from firedrake import * + +1D Periodic Poisson Problem +--------------------------- +We solve the 1D periodic Poisson problem + +.. math:: + + - \Delta u = \sin(x), \quad u(0) = u(2 \pi), + +on :math:`\Omega = [0, 2 \pi]`. This problem has a trivial nullspace of constants; if we fix the constant to be zero then this problem has an analytical solution of :math:`u(x) = \sin(x)`. +We use the Firedrake function ``PeriodicIntervalMesh`` to build the mesh and implement the usual weak formulation of the Poisson equation: :: + + mesh = PeriodicIntervalMesh(32, 2*pi) + x, = SpatialCoordinate(mesh) + + V = FunctionSpace(mesh, "CG", 1) + u = TrialFunction(V) + v = TestFunction(V) + + uh = Function(V, name="Numerical") + u_exact = sin(x) + + a = dot(grad(u), grad(v)) * dx + L = inner(sin(x), v) * dx + + problem = LinearVariationalProblem(a, L, uh) + +We use the ``VectorSpaceBasis`` function to tell PETSc that the problem has a +constant nullspace:: + + nullspace = VectorSpaceBasis( + constant=True, + comm=mesh.comm, + ) + + solver = LinearVariationalSolver( + problem, + nullspace=nullspace, + transpose_nullspace=nullspace, + ) + + solver.solve() + + print(f"L2 error = {errornorm(u_exact, uh, norm_type='L2'):.3e}") + +We plot the solution below. + + +.. figure:: Example1.png + :align: center + :alt: Finite element solution to 1D periodic Poisson problem. + + Finite element solution to 1D periodic Poisson problem. + +2D Periodic Poisson Problem +--------------------------- + +In two dimensions, we can choose which boundaries are periodic. We consider +an :math:`x`-periodic Poisson problem, + +.. math:: + + - \Delta u = (1 + \pi^2) \sin(x) \sin(\pi y), \quad u(0,y) = u(2 \pi, y), \quad u(x, 0) = u(x, 1) = 0, + +on :math:`\Omega = [0, 2 \pi] \times [0,1]`. The nullspace is fixed for this problem. The exact solution is :math:`u(x,y) = \sin(x)\sin(\pi y)`. +We use the Firedrake function ``PeriodicRectangleMesh`` to build the mesh: if you call this function without calling the argument ``PeriodicRectangleMesh`` then the rectangle will be periodic in both :math:`x` and :math:`y`. +In our case we only want periodicity in :math:`x`. :: + + nx = 32 + ny = 16 + + mesh = PeriodicRectangleMesh(nx, ny, 2*pi, 1.0, direction="x") + x, y = SpatialCoordinate(mesh) + + V = FunctionSpace(mesh, "CG", 1) + u = TrialFunction(V) + v = TestFunction(V) + uh = Function(V, name="Numerical") + u_exact = sin(x) * sin(pi*y) + + f = (1 + pi**2)*sin(x)*sin(pi*y) + a = inner(grad(u), grad(v))*dx + L = f*v*dx + + # Homogeneous Dirichlet BCs on y = 0 and y = 1. + bcs = DirichletBC(V, 0.0, (3, 4)) + + solve(a == L, uh, bcs=bcs) + + print(f"L2 error = {errornorm(u_exact, uh, norm_type='L2'):.3e}") + VTKFile("output/rectangle_poisson.pvd").write(uh) + + +We plot the solution below. + + +.. figure:: Example2.png + :align: center + :alt: Finite element solution to 2D periodic Poisson problem. + + Finite element solution to 2D periodic Poisson problem. + +3D Periodic Helmholtz Problem +----------------------------- + +For the simplest 3D shapes (cube and cuboid) Firedrake is able to generate tetrahedral periodic meshes using the commands ``PeriodicUnitCubeMesh`` and ``PeriodicBoxMesh`` respectively. +For cylindrical objects we have to construct a 2D cross-section and extrude, forming a mesh with prismatic elements. To construct any periodic mesh with tetrahedral elements we need to use Netgen, which is discussed in the next section. + +For now we solve the 3D :math:`z`-periodic Helmholtz problem + +.. math:: + + (I - \Delta) u = f , \quad u(x,y,0) = u(x, y, 2 \pi), \quad u\vert_{\text{walls}} = 0, + +on the cylindrical domain + +.. math:: + + \Omega = + \left\{ + (x,y,z)\in\mathbb{R}^3 : + x^2+y^2\le 1,\; + 0\le z\le 2\pi + \right\}. + +We can manufacture a right-hand side by choosing the solution +:math:`u_{\text{exact}}(x,y,z) = (1-x^2-y^2)\cos(z)` which gives :math:`f = (I - \Delta) u_{\text{exact}}`. + + +We build a 2D mesh for the desired cross section of the cylinder and then we use the command ``ExtrudedMesh`` with flag ``periodic=True`` to generate a mesh of prisms. :: + + refinements = 2 + base = UnitDiskMesh(refinements) + + mesh = ExtrudedMesh( + base, + layers=32, + layer_height=2*pi/32, + periodic=True, + ) + +The two important arguments are ``layers`` and ``layer_height``, which tell Firedrake how far to extrude the mesh. By default ``layer_height`` is 1/``layers``. :: + + x, y, z = SpatialCoordinate(mesh) + + V = FunctionSpace(mesh, "CG", 1) + + u = TrialFunction(V) + v = TestFunction(V) + + uh = Function(V, name="Numerical") + u_exact = (1 - x*x - y*y) * cos(z) + + f = u_exact - div(grad(u_exact)) + + u = TrialFunction(V) + v = TestFunction(V) + a = (inner(u, v) + inner(grad(u), grad(v))) * dx + L = inner(f, v) * dx + + + bc = DirichletBC(V, 0.0, "on_boundary") + + solve(a == L, uh, bcs=bc) + + print(f"L2 error = {errornorm(u_exact, uh, norm_type='L2'):.3e}") + VTKFile("output/cylinder_fd_Helmholtz.pvd").write(uh) + +We plot the solution below. + + +.. figure:: Example3.png + :align: center + :alt: Cross-section of the finite element solution to a 3D periodic Helmholtz problem. + + Cross-section of the finite element solution to a 3D periodic Helmholtz problem using prismatic elements. + +Periodic Meshes From Netgen +--------------------------- + +Netgen can help us generate a periodic mesh of tetrahedral elements. + +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 ``PeriodicBoxMesh``, and it is also available for any Netgen geometry carrying +periodic identifications. + +Identifications are declared on the geometry, before meshing, with the OCC +``Identify`` method: + +.. code-block:: python + + 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. + +To construct a periodic cylinder of length :math:`2\pi` we identify the two end caps of the cylinder 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 + + height = 2 * pi + cyl = Cylinder(Pnt(0, 0, 0), Z, r=1.0, h=height) + + # 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 + cyl.faces.Min(Z).Identify(cyl.faces.Max(Z), "toroidal", + IdentificationType.PERIODIC, + gp_Trsf.Translation(gp_Vec(0, 0, height))) + ngmsh = OCCGeometry(cyl).GenerateMesh(maxh=0.4) + msh = Mesh(ngmsh) + +.. 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 ``z``-direction. We can verify the geometry survived the +merge intact --- the volume of the cylinder is :math:`\pi r^2 h = 2\pi^2`: :: + + volume = assemble(Constant(1.0) * dx(domain=msh)) + print( + f"cylinder volume: {volume:.4f} " + f"(exact 2*pi**2 = {2*pi**2:.4f})" + ) + +We solve again the same Helmholtz problem as in the previous section. :: + + V = FunctionSpace(msh, "CG", 2) + x, y, z = SpatialCoordinate(msh) + u_exact = cos(z) * (1 - x**2 - y**2) + f = u_exact - div(grad(u_exact)) + + u = TrialFunction(V) + v = TestFunction(V) + a = (inner(u, v) + inner(grad(u), grad(v))) * dx + L = inner(f, v) * dx + +We look up the id of the ``"wall"`` boundary with ``GetRegionNames`` :: + + labels = [i + 1 for i, name in enumerate(ngmsh.GetRegionNames(codim=1)) if name == "wall"] + bc = DirichletBC(V, 0, labels) + + uh = Function(V) + solve(a == L, uh, bcs=bc) + + print(f"L2 error = {errornorm(u_exact, uh, norm_type='L2'):.3e}") + VTKFile("output/cylinder_helmholtz.pvd").write(uh) + + +We plot the solution below with a "crinkle cut" cross-section to inspect the solution. + + +.. figure:: Example4.png + :align: center + :alt: Finite element solution to 3D periodic Helmholtz problem using tetrahedral elements. + + Finite element solution to 3D periodic Helmholtz problem. + +Netgen Tokamak Example +--------------------------- + +Helmholtz problems on periodic cylinders are used in simplified models of magnetic confinement fusion devices, where the straight cylinder acts as a large-aspect-ratio approximation of a torus. +To obtain a more realistic geometry, we construct the domain in cylindrical coordinates and choose a cross-section representative of a tokamak plasma. After solving we transform the solution to Cartesian coordiantes. + +We describe the tokamak in cylindrical coordinates :math:`(R,\phi,Z)` where :math:`\phi` is the toroidal angle. +The cross-section is constructed in the :math:`(R,Z)` plane. Rather than extruding by :math:`\phi` we extrude by the corresponding toroidal arc-length + +.. math:: + + s = R_0 \phi, + +where :math:`R_0` is called the major radius. This makes one complete revolution of the tokamak have length :math:`2 \pi R_0`. +We parametrize the tokamak cross-section using the formula + +.. math:: + + R(\theta) = R0 + a \cos(\theta + \sin^{-1}(\delta) \sin(\theta)), \quad Z(\theta) = a \kappa \sin(\theta). + +where :math:`\delta` is the triangularity, :math:`\kappa` is the elongation, :math:`a` is the tokamak minor radius, and :math:`\theta \in [0, 2\pi]`. :: + + from netgen.occ import ( + OCCGeometry, WorkPlane, Axes, Pnt, Z, X, + gp_Trsf, gp_Vec + ) + from netgen.meshing import IdentificationType + import math + + # Geometry parameters - large aspect ratio tokamak + R0 = 3.0 + a = 1.0 + kappa = 2.0 + delta = 0.3 + + n_boundary = 40 + maxh = 0.5 + + # Build a tokamak cross section and extrude periodically + + alpha = math.asin(delta) + + boundary_points = [] + for i in range(n_boundary): + theta = 2.0 * pi * i / n_boundary + R = R0 + a * math.cos(theta + alpha * math.sin(theta)) + Zc = kappa * a * math.sin(theta) + boundary_points.append((R, Zc)) + + wp = WorkPlane(Axes((0, 0, 0), n=Z, h=X)) + + # Start at the first boundary point, then draw a closed polyline + R0p, Z0p = boundary_points[0] + wp.MoveTo(R0p, Z0p) + for R, Zc in boundary_points[1:]: + wp.LineTo(R, Zc) + wp.Close() + + face = wp.Face() + + # Extrude in periodic phi direction and identify end caps. + + height = 2 * pi * R0 + solid = face.Extrude(height * Z) + + # Side wall(s) + for f in solid.faces: + f.name = "wall" + + bottom = solid.faces.Min(Z) + top = solid.faces.Max(Z) + bottom.name = "bottom" + top.name = "top" + + # Periodic identification of the end caps + bottom.Identify( + top, + "periodic_phi", + IdentificationType.PERIODIC, + gp_Trsf.Translation(gp_Vec(0, 0, height)), + ) + + # Mesh and convert to Firedrake + ngmsh = OCCGeometry(solid).GenerateMesh(maxh=maxh) + msh = Mesh(ngmsh, name="TokamakPeriodic") + +The physical cylindrical coordinate convention is :math:`(R,\phi,Z)`. +However, the Netgen mesh stores its coordinate components in the order +:math:`(R,Z,s)`: the first two components are the coordinates of the +cross-section and the third is the toroidal arc-length. We therefore unpack +``SpatialCoordinate`` as ``R, Zc, s``. + +We solve the Helmholtz problem + +.. math:: + + (I-\Delta)u=f, \quad u(r, 0, z) = u(r, 2 \pi R_0, z), \quad u|_{\text{walls}} = 0, + +with source term +:math:`f(R,s,Z)=RZ\cos(s)`. :: + + V = FunctionSpace(msh, "CG", 2) + R, Zc, s = SpatialCoordinate(msh) + + f = R * Zc * cos(s) + + 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) + + uh = Function(V, name="Solution") + solve(a == L, uh, bcs=bc) + VTKFile("output/tokamak_helmholtz.pvd").write(uh) + +We plot the output in :math:`(R, s, Z)`, a cylindrical representation: + + +.. figure:: Example5.png + :align: center + :width: 85% + :alt: Finite element solution to 3D periodic Helmholtz problem in cylindrical tokamak geometry with tetrahedral elements. + + Finite element solution to 3D periodic Helmholtz problem in cylindrical tokamak geometry. + +Since the solution is represented as a function on the mesh, we can obtain a +Cartesian visualisation by interpolating the cylindrical-to-Cartesian map + +.. math:: + + (R,s,Z) \longmapsto + \bigl( R\cos(s / R_0), R\sin(s/ R_0), Z \bigr) + +into the mesh coordinate field. The mesh components must again be unpacked in +their stored order, :math:`(R,Z,s)`:: + + R, Zc, s = SpatialCoordinate(msh) + msh.coordinates.interpolate(as_vector(( + R*cos(s/R0), + R*sin(s/R0), + Zc, + ))) + VTKFile("TokamakCartesianSolution.pvd").write(uh) + +We plot a cross-section "crinkle cut" of the solution in Cartesian coordinates: + +.. figure:: Example6.png + :align: center + :width: 105% + :alt: Finite element solution to 3D periodic Helmholtz problem in Cartesian tokamak geometry. + + Cross-section of a finite element solution to a 3D periodic Helmholtz problem in Cartesian tokamak geometry with tetrahedral elements. + diff --git a/docs/source/intro_tut.rst b/docs/source/intro_tut.rst index 93b3d3afeb..fe9b46e3b8 100644 --- a/docs/source/intro_tut.rst +++ b/docs/source/intro_tut.rst @@ -17,3 +17,4 @@ Introductory Tutorials A linear shallow water equations example using a Strang timestepping scheme. A linear wave equation with optional mass lumping. Creating Firedrake-compatible meshes in Gmsh. + Creating periodic meshes. diff --git a/firedrake/mesh.py b/firedrake/mesh.py index 881cdfe523..a3bf059fdc 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -3076,6 +3076,78 @@ def curve_field(self, order, permutation_tol=None, 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): + '''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. + ''' + 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] + )[0] + data[cell_nodes[i]] = unwrapped[e][permutation] + return new_coordinates + @PETSc.Log.EventDecorator() def make_mesh_from_coordinates(coordinates, name, tolerance=0.5): @@ -3324,6 +3396,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): @@ -3381,6 +3460,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}: @@ -3393,6 +3473,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) @@ -3416,7 +3501,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"), @@ -3428,14 +3517,31 @@ 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: - permutation_tol = netgen_flags.get("permutation_tol", None) + 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 + coordinates = mesh._periodic_coordinates() + 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 + elif degree != 1: cg = netgen_flags.get("cg", None) coordinates = mesh.curve_field( order=degree, - permutation_tol=permutation_tol, cg_field=cg, ) # Do not redistribute the mesh diff --git a/firedrake/netgen.py b/firedrake/netgen.py index c0c0369779..ff56a0ef59 100644 --- a/firedrake/netgen.py +++ b/firedrake/netgen.py @@ -70,12 +70,18 @@ def netgen_distribute(V: firedrake.functionspaceimpl.WithGeometryBase, for i in np.ndindex(V.shape): di = netgen_data[(..., *i)].flatten() vec0[:len(di)] = di - _, vec = plex.distributeField(sf, section0, vec0) - arr = vec.getArray() + section_i, vec_i = plex.distributeField(sf, section0, vec0) + arr = vec_i.getArray() if plex_data is None: plex_data = np.empty(arr.shape + V.shape, dtype=dtype) plex_data[(..., *i)] = arr + section_i.destroy() + vec_i.destroy() plex_data = plex_data.reshape(-1, *nshape[1:]) + # Destroy the remaining transients created in this call. + sfBCInv.destroy() + section0.destroy() + vec0.destroy() return plex_data diff --git a/tests/firedrake/demos/test_demos_run.py b/tests/firedrake/demos/test_demos_run.py index f53a611c2b..0b1a9321b2 100644 --- a/tests/firedrake/demos/test_demos_run.py +++ b/tests/firedrake/demos/test_demos_run.py @@ -43,6 +43,7 @@ Demo(("netgen", "netgen_mesh"), ["mumps", "netgen", "slepc", "vtk"]), Demo(("nonlinear_QG_winddrivengyre", "qg_winddrivengyre"), ["vtk"]), Demo(("parallel-printing", "parprint"), []), + Demo(("periodic_meshes", "periodic_meshes"), ["netgen", "vtk"]), Demo(("poisson", "poisson_mixed"), ["vtk"]), Demo(("patch", "poisson_mg_patches"), []), Demo(("patch", "stokes_vanka_patches"), []), diff --git a/tests/firedrake/regression/test_netgen.py b/tests/firedrake/regression/test_netgen.py index 64abf21c31..d0a6e305d5 100644 --- a/tests/firedrake/regression/test_netgen.py +++ b/tests/firedrake/regression/test_netgen.py @@ -362,3 +362,90 @@ def adapt(mesh, eta): 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)