Feature/netgen periodic meshes by UZerbinati · Pull Request #5217 · firedrakeproject/firedrake

UZerbinati · 2026-07-05T17:49:46Z

When a Netgen mesh carries periodic identifications, ngsPETSc now produces a vertex-merged (periodic) DMPlex. This finishes the job on the Firedrake side: the Mesh() netgen branch detects periodicity, builds the intermediate topology without reordering, and attaches a discontinuous DG1 coordinate field (MeshGeometry._periodic_coordinates) carrying the un-wrapped geometry via _postprocess_periodic_mesh, exactly as Firedrake's own periodic meshes. Each Firedrake cell is matched to its Netgen element by geometry (hashing the wrapped per-cell coordinates), which is robust to reordering. The DG coordinate space uses the equispaced variant required by _set_dg_coordinates. High-order curving of periodic meshes raises NotImplementedError for now. Adds regression tests for a periodic square (one and two directions), a periodic cylinder, and the too-coarse error path. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

Add a "Periodic Meshes" section to the Netgen integration demo showing how to declare OCC periodic identifications and import the resulting mesh into Firedrake. The worked example is the periodic-cylinder ("screw pinch") reduced model of a tokamak: a cylinder whose two circular ends are identified, making the axial (toroidal) direction periodic. It labels the lateral wall, checks the geometry survives the vertex merge (volume 2*pi^2), and solves a manufactured Helmholtz problem with Dirichlet data on the wall and periodicity on the ends, demonstrating seam continuity. Also documents the too-coarse ValueError and the degree==1 limitation. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

pbrubeck · 2026-07-09T19:09:08Z

+    topology_reorder = False if netgen_periodic else reorder
+    topology = MeshTopology(plex, name=plex.getName(), reorder=topology_reorder,


Suggested change

topology_reorder = False if netgen_periodic else reorder

topology = MeshTopology(plex, name=plex.getName(), reorder=topology_reorder,

reorder = False if netgen_periodic else reorder

topology = MeshTopology(plex, name=plex.getName(), reorder=reorder,

pbrubeck · 2026-07-09T19:10:21Z

+    @PETSc.Log.EventDecorator()
+    def _periodic_coordinates(self, permutation_tol=1e-8):


Move this to netgen.py

pbrubeck · 2026-07-09T19:12:27Z

+        # 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


Do as the comment says and supress reordering here. Remove the comment, we don't like AI slop

pbrubeck · 2026-07-09T19:13:19Z

+            # discontinuous coordinate field carrying the un-wrapped geometry.
+            if degree != 1:
+                raise NotImplementedError(
+                    "High-order curving of periodic netgen meshes is not supported yet."


Support this, I support you

pbrubeck · 2026-07-09T19:14:14Z

+            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


create a helper function in netgen.py doing everything here

…ic-meshes

…ry tutorials. This takes some of the demo from netgen_meshes to make it less lengthy.

pbrubeck · 2026-08-03T12:25:04Z

+
+   from netgen.occ import Cylinder, OCCGeometry, Pnt, Z, gp_Trsf, gp_Vec
+   from netgen.meshing import IdentificationType
+   from math import pi as PI


Suggested change

from math import pi as PI

…updated figures.

pbrubeck · 2026-08-03T23:11:21Z


+    @PETSc.Log.EventDecorator()
+    def _periodic_coordinates(self):
+        '''Return a discontinuous coordinate field for a periodic netgen mesh.


We should not add a netgen-only function to this file.

pbrubeck · 2026-08-03T23:12:36Z

+        # 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]):


This is an antipattern (see AGENTS.md)

pbrubeck · 2026-08-03T23:13:34Z

This work should be properly re-implemented on top of #5277 and NGSolve/ngsPETSc#125

…lem. I have rewritten verbose sentences to make the demo clearer and easier to use

UZerbinati

@ThighamOxford any chance you would like to take a hit also at tackling the other comments Pablo left ? I'll be thinking about/working on high order :)

UZerbinati · 2026-08-04T21:57:50Z

+
+    - \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)`.


Suggested change

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)`.

on :math:`\Omega = [0, 2 \pi]`. This problem has a trivial nullspace of constants, thus if we fix the constant to be zero then this problem has an analytical solution of :math:`u(x) = \sin(x)`.

UZerbinati · 2026-08-04T21:58:53Z

+
+    - \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)`.


What do you mean when you say that the nullspace is fixed ? Does the problem have a nullspace ?

The 1D periodic poisson problem with periodic data at the boundary has a nullspace of constants. I then clarify that for this 2D periodic poisson problem the nullspace is fixed due to additional dirichlet boundary conditions on the top and bottom of the box.

UZerbinati · 2026-08-04T22:00:24Z

+   :alt: Finite element solution to 2D periodic Poisson problem.
+
+   Finite element solution to 2D periodic Poisson problem.
+


If I remember correctly there are a bunch of features that do not work on extruded meshes, it would be nice to point to this issue, no ? Probably more detail can be found in the documentation for extruded meshes ?

UZerbinati · 2026-08-04T22:02:08Z

+
+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.


Suggested change

automatically -- no extra flag on the ``Mesh`` constructor is required.

automatically: no extra flag on the ``Mesh`` constructor is required.

UZerbinati · 2026-08-04T22:02:54Z

+.. 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``


Suggested change

whole period, its two ends are identified and the cell collapses; Firedrake then raises a ``ValueError``

whole period, its two ends are identified and the cell collapses. Firedrake then raises a ``ValueError``

UZerbinati · 2026-08-04T22:03:16Z

+   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.


I need to work on fixing this :)

UZerbinati · 2026-08-04T22:03:38Z

+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`: ::


Suggested change

merge intact --- the volume of the cylinder is :math:`\pi r^2 h = 2\pi^2`: ::

merge intact: the volume of the cylinder is :math:`\pi r^2 h = 2\pi^2`: ::

UZerbinati · 2026-08-04T22:08:08Z

+
+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


I believe there is no need to use gp_Vec tuple should be accepted now (bad Claude...).
Also Translation should be importable without gp_Trsf.

UZerbinati · 2026-08-04T22:15:01Z

+   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",


I'm pretty sure there is also a way of identifying surfaces in OCC (at least it is there for CSG3D: see here)that does not require to specify IndentificationType and transformation. Maybe it would be useful to showcase that. Ahh ! Found it here

Co-authored-by: Umberto Zerbinati <umberto.zerbinati@maths.ox.ac.uk>

Co-authored-by: Pablo Brubeck <brubeck@protonmail.com>

Umberto Zerbinati and others added 2 commits June 25, 2026 14:16

UZerbinati marked this pull request as draft July 5, 2026 17:50

UZerbinati mentioned this pull request Jul 5, 2026

Support importing periodic Netgen meshes NGSolve/ngsPETSc#124

Draft

pbrubeck reviewed Jul 9, 2026

View reviewed changes

Umberto Zerbinati and others added 7 commits July 10, 2026 16:07

Destroy object the garbage collector would have issue with

1f5925c

Merge branch 'uz/netgen_garbage_collector' into feature/netgen-period…

fc9e9ea

…ic-meshes

merge conflict

4407328

Merge origin/main into feature/netgen-periodic-meshes

8bd5b23

Resolve merge conflict in Netgen demo

ab28e92

fix permutation_tol

7557a73

Written a periodic_meshes demo for Firedrake and placed in introducto…

6a46b63

…ry tutorials. This takes some of the demo from netgen_meshes to make it less lengthy.

pbrubeck reviewed Aug 3, 2026

View reviewed changes

ThighamOxford added 3 commits August 3, 2026 13:26

Added periodic_meshes to test_demos_run.py

3426bd6

Fixed code block that was not intedned to be run in test

68bf2ba

I have edited the demo to match the Firedrake code style, and I have …

0ffaaa3

…updated figures.

pbrubeck reviewed Aug 3, 2026

View reviewed changes

Comment thread demos/periodic_meshes/periodic_meshes.py.rst Outdated

ThighamOxford added 2 commits August 4, 2026 11:55

I have simplified exposition by making examples 3 and 4 the same prob…

5410681

…lem. I have rewritten verbose sentences to make the demo clearer and easier to use

Fixed a code typo

a077ec3

UZerbinati commented Aug 4, 2026

View reviewed changes

ThighamOxford and others added 2 commits August 6, 2026 16:02

Update demos/periodic_meshes/periodic_meshes.py.rst

da596ef

Co-authored-by: Umberto Zerbinati <umberto.zerbinati@maths.ox.ac.uk>

Update demos/periodic_meshes/periodic_meshes.py.rst

e565df2

Co-authored-by: Pablo Brubeck <brubeck@protonmail.com>

		topology_reorder = False if netgen_periodic else reorder
		topology = MeshTopology(plex, name=plex.getName(), reorder=topology_reorder,

		@PETSc.Log.EventDecorator()
		def _periodic_coordinates(self, permutation_tol=1e-8):


		- \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)`.


		- \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)`.

		:alt: Finite element solution to 2D periodic Poisson problem.

		Finite element solution to 2D periodic Poisson problem.

	automatically -- no extra flag on the ``Mesh`` constructor is required.
	automatically: no extra flag on the ``Mesh`` constructor is required.

	whole period, its two ends are identified and the cell collapses; Firedrake then raises a ``ValueError``
	whole period, its two ends are identified and the cell collapses. Firedrake then raises a ``ValueError``

	merge intact --- the volume of the cylinder is :math:`\pi r^2 h = 2\pi^2`: ::
	merge intact: the volume of the cylinder is :math:`\pi r^2 h = 2\pi^2`: ::


		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

Uh oh!

Conversation

UZerbinati commented Jul 5, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

pbrubeck commented Aug 3, 2026

Uh oh!

Uh oh!

UZerbinati left a comment • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

UZerbinati commented Jul 5, 2026 •

edited

Loading

UZerbinati left a comment •

edited

Loading