From 8bd9464c33eca1dae7ea4a6b123ce50c82e8003a Mon Sep 17 00:00:00 2001 From: ThighamOxford Date: Fri, 31 Jul 2026 12:54:58 +0100 Subject: [PATCH 1/3] Resolve remaining Netgen demo conflict --- demos/netgen/netgen_mesh.py.rst | 100 ++++++++++++++++++++++++++++++-- 1 file changed, 94 insertions(+), 6 deletions(-) diff --git a/demos/netgen/netgen_mesh.py.rst b/demos/netgen/netgen_mesh.py.rst index 88b2c1fafa..48915da475 100755 --- a/demos/netgen/netgen_mesh.py.rst +++ b/demos/netgen/netgen_mesh.py.rst @@ -1,7 +1,5 @@ Netgen integration in Firedrake -=============================== - -This tutorial was contributed by `Umberto Zerbinati `__ and `Patrick E. Farrell `__. This tutorial is based on the documentation for Netgen in the `NGSolve Docs `__. +=================This tutorial was contributed by `Umberto Zerbinati `__ and `Patrick E. Farrell `__. This tutorial is based on the documentation for Netgen in the `NGSolve Docs `__. The purpose of this demo is to summarise how to construct and use a Netgen mesh in Firedrake. Netgen is an automated 2D/3D open-source mesher developed by Joachim Schöberl, which is capable of dealing both with simplicial and quadrilateral dominated mesh. @@ -452,7 +450,7 @@ and identify the two end caps by a translation of :math:`2\pi` along ``Z``:: 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) + VTKFile("output/PeriodicCylinder.pvd").write(msh) .. warning:: @@ -499,14 +497,104 @@ Poisson example above) and manufacture the right-hand side :math:`f = u_{\text{e sol = Function(V) solve(a == L, sol, bcs=bc) - VTKFile("output/TokamakSolution.pvd").write(sol) + VTKFile("output/PeriodicCylinderSolution.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 +The recovered solution is continuous across the identified ends: opening ``output/PeriodicCylinderSolution.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. +We can also solve a Helmholtz problem on a realistic tokamak geometry, although we no longer have an analytical solution to test against. +We now work in cylindrical coordinate space (R, phi, Z) where we construct a cross-section in the (R,Z) plane with the appropriate tokamak shape and then extrude in phi by :math:`2\pi`. + + from firedrake import * + from netgen.occ import ( + OCCGeometry, WorkPlane, Axes, Pnt, Z, X, + gp_Trsf, gp_Vec + ) + from netgen.meshing import IdentificationType + from math import pi as PI + import math + + # Geometry parameters - large aspect ratio tokamak + R0 = 3.0 + a = 1.0 + kappa = 2.0 + delta = 0.3 + + n_boundary = 40 + mesh_size = 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 * math.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. + + phi_length = 2 * PI + solid = face.Extrude(phi_length * 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, phi_length)), + ) + + # Mesh and convert to Firedrake + ngmsh = OCCGeometry(solid).GenerateMesh(maxh=mesh_size) + msh = Mesh(ngmsh, name="TokamakPeriodic") + VTKFile("output/Tokamak_3d_Mesh.pvd").write(msh) + +Now we solve a Helmholtz problem again with source :math:`f = \cos(z)xy`. + + V = FunctionSpace(msh, "CG", 2) + x, y, z = SpatialCoordinate(msh) + + f = cos(z)* x*y + + 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) From 1f7b23db235d8381dd55d28cbaeefecf6aae2d0e Mon Sep 17 00:00:00 2001 From: Thomas Higham Date: Fri, 31 Jul 2026 13:00:07 +0100 Subject: [PATCH 2/3] Apply suggestion from @ThighamOxford --- demos/netgen/netgen_mesh.py.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/demos/netgen/netgen_mesh.py.rst b/demos/netgen/netgen_mesh.py.rst index 48915da475..a7acf9c746 100755 --- a/demos/netgen/netgen_mesh.py.rst +++ b/demos/netgen/netgen_mesh.py.rst @@ -1,5 +1,5 @@ Netgen integration in Firedrake -=================This tutorial was contributed by `Umberto Zerbinati `__ and `Patrick E. Farrell `__. This tutorial is based on the documentation for Netgen in the `NGSolve Docs `__. +This tutorial was contributed by `Umberto Zerbinati `__ and `Patrick E. Farrell `__. This tutorial is based on the documentation for Netgen in the `NGSolve Docs `__. The purpose of this demo is to summarise how to construct and use a Netgen mesh in Firedrake. Netgen is an automated 2D/3D open-source mesher developed by Joachim Schöberl, which is capable of dealing both with simplicial and quadrilateral dominated mesh. From ece6c14ddb2bec1be74a3dfa784fcd9c5e4599b8 Mon Sep 17 00:00:00 2001 From: ThighamOxford Date: Fri, 31 Jul 2026 13:02:52 +0100 Subject: [PATCH 3/3] Fix title underline formatting in Netgen demo --- demos/netgen/netgen_mesh.py.rst | 2 ++ 1 file changed, 2 insertions(+) diff --git a/demos/netgen/netgen_mesh.py.rst b/demos/netgen/netgen_mesh.py.rst index a7acf9c746..008703c5fa 100755 --- a/demos/netgen/netgen_mesh.py.rst +++ b/demos/netgen/netgen_mesh.py.rst @@ -1,4 +1,6 @@ Netgen integration in Firedrake +=============================== + This tutorial was contributed by `Umberto Zerbinati `__ and `Patrick E. Farrell `__. This tutorial is based on the documentation for Netgen in the `NGSolve Docs `__. The purpose of this demo is to summarise how to construct and use a Netgen mesh in Firedrake.