# 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. 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. We will first have a look at how different geometries can be implemented using Netgen, then we will consider as an example the Poisson problem and discuss boundary conditions. Finally we will show how to use mesh refinement features included in Netgen to construct an adaptive scheme.

## Installing Netgen¶

This demo requires the NGSolve/Netgen suite to be installed. This is most easily achieved by providing the optional –netgen flag to either firedrake-install (for a new installation), or firedrake-update (to add the NGSolve/Netgen suite to an existing installation).

from firedrake import *
try:
import netgen
except ImportError:
import sys
warning("Unable to import NetGen.")
sys.exit(0)


## Constructive Solid Geometry¶

Using the Constructive Solid Geometry (CSG) features implemented in Netgen, we can construct a geometry starting from basic sets such as rectangles and disks. Let’s have a look at an example. First we need to import the Netgen library required in order to construct a SplineGeometry:

from netgen.geom2d import SplineGeometry
geo = SplineGeometry()


Now using the two of the predefined CSG geometries included in Netgen, a rectangle and a circle, we construct a square with an inscribed disk that we will later flag for refinement, using the SetMaterial and SetDomainMaxH methods:

geo = SplineGeometry()
p2=(1, 1),
bc="rectangle",
leftdomain=1,
rightdomain=0)
r=0.5,
bc="circle",
leftdomain=2,
rightdomain=1)

# Flagging for the inside of the disk with a different material IDs
geo.SetMaterial(1, "outer")
geo.SetMaterial(2, "inner")
geo.SetDomainMaxH(2, 0.02)


Notice that the leftdomain and rightdomain attribute in the AddRectangle and AddCircle methods are used to set a domain index for the domain on the left and right side respectively of the rectangle and circle perimeter. It is worth mentioning that the perimeters are parametrised in a counterclockwise direction. We can now construct a mesh for the geometry we have defined and save it to a PVD file for visualisation. We will do so using the GenerateMesh method inside of the SplineGeometry class:

ngmsh = geo.GenerateMesh(maxh=0.1)
# Generating a Firedrake mesh from the NetGen mesh
msh = Mesh(ngmsh)
File("output/MeshExample1.pvd").write(msh)


We can also “hardwire” a geometry in Netgen, specifying the points and the splines making the boundary of the domain. In particular we add points to an instance of the SplineGeometry class using the method AppendPoint, while a spline can be added using the Append method. When using the Append method we can specify the type of Spline and its boundary label (in this particular case we use the “bnd” tag for all the curves making up the boundary). Here is showed how to construct a “Pacman”-like domain:

geo = SplineGeometry()
pnts = [(0, 0), (1, 0), (1, 1),
(0, 1), (-1, 1), (-1, 0),
(-1, -1), (0, -1)]
p1, p2, p3, p4, p5, p6, p7, p8 = [geo.AppendPoint(*pnt) for pnt in pnts]
curves = [[["line", p1, p2], "line"],
[["spline3", p2, p3, p4], "curve"],
[["spline3", p4, p5, p6], "curve"],
[["spline3", p6, p7, p8], "curve"],
[["line", p8, p1], "line"]]
[geo.Append(c, bc=bc) for c, bc in curves]
ngmsh = geo.GenerateMesh(maxh=0.4)
msh = Mesh(ngmsh)
File("output/MeshExample2.pvd").write(msh)


## Example: Poisson Problem¶

Let’s now have a look at some features that can be useful when solving a variational problem using a Netgen mesh. We will consider as example the Poisson problem with constant source data and homogeneous boundary conditions on the boundary of the PacMan. The only method new to the reader should be the GetRegionNames which allows to find the IDs of the boundary we have labeled in Netgen. As usual we begin defining the FunctionSpace that will be used in our discretisation, defining trial and test functions and the source data:

V = FunctionSpace(msh, "CG", 1)
u = TrialFunction(V)
v = TestFunction(V)

f = Function(V)
x, y = SpatialCoordinate(msh)
f.assign(1)


Now we can define the bilinear form and linear function that characterize the weak formulation of the Poisson problem in abstract form, i.e.

$\text{find } u\in H^1_0(\Omega) \text{ s.t. } a(u,v) := \int_{\Omega} \nabla u\cdot \nabla v \; d\vec{x} = L(v) := \int_{\Omega} fv\; d\vec{x}\qquad v\in H^1_0(\Omega).$

In code this becomes:

a = inner(grad(u), grad(v))*dx
L = inner(f, v) * dx


Now we are ready to assemble the stiffness matrix for the problem. Since we want to enforce Dirichlet boundary conditions we construct a DirichletBC object and we use the GetRegionNames method from the Netgen mesh in order to map the label we have given when describing the geometry to the PETSc DMPLEX IDs. In particular if we look for the IDs of boundary element labeled either “line” or “curve” we would get:

labels = [i+1 for i, name in enumerate(ngmsh.GetRegionNames(codim=1)) if name in ["line","curve"]]
bc = DirichletBC(V, 0, labels)
print(labels)


We then proceed to solve the problem:

sol = Function(V)
solve(a == L, sol, bcs=bc)
File("output/Poisson.pvd").write(sol)


## Mesh Refinement¶

In this section we will discuss how to use the mesh refinement methods wrapped from Netgen C++ interface. In particular we will be considering a Laplace eigenvalue problem on the same PacMan domain presented above, i.e.:

$\text{Find } u \in H^1_0(\Omega) \text{ and } \lambda \in \mathbb{R} \text{ s.t. } \int_{\Omega} \nabla u\cdot\nabla v\;d\vec{x} = \lambda \int_{\Omega}uv\;d\vec{x}\qquad \forall v\in H^1_0(\Omega).$

This script is based on a code developed by Professor Daniele Boffi and based on a code from Professor Douglas Arnold for the source problem. We begin by defining some quantities of interest such as the desired tolerance, the maximum number of iterations and the exact eigenvalue:

from firedrake.petsc import PETSc
from slepc4py import SLEPc
import numpy as np

tolerance = 1e-16
max_iterations = 10
exact = 3.375610652693620492628**2


We create a function to solve the eigenvalue problem using SLEPc. We begin initialising the FunctionSpace, the bilinear forms and linear functionals needed in the variational problem. Then a SLEPc Eigenvalue Problem Solver (EPS) is initialised and set up to use a shift and invert (SINVERT) spectral transformation where the preconditioner factorisation is computed using MUMPS:

def Solve(msh, labels):
V = FunctionSpace(msh, "CG", 2)
u = TrialFunction(V)
v = TestFunction(V)
m = (u*v)*dx
uh = Function(V)
bc = DirichletBC(V, 0, labels)
A = assemble(a, bcs=bc)
M = assemble(m, bcs=bc, weight=0.)
Asc, Msc = A.M.handle, M.M.handle
E = SLEPc.EPS().create()
E.setType(SLEPc.EPS.Type.ARNOLDI)
E.setProblemType(SLEPc.EPS.ProblemType.GHEP)
E.setDimensions(1, SLEPc.DECIDE)
E.setOperators(Asc, Msc)
ST = E.getST()
ST.setType(SLEPc.ST.Type.SINVERT)
PC = ST.getKSP().getPC()
PC.setType("lu")
PC.setFactorSolverType("mumps")
E.setST(ST)
E.solve()
vr, vi = Asc.getVecs()
with uh.dat.vec_wo as vr:
lam = E.getEigenpair(0, vr, vi)
return (lam, uh, V)


We will also need a function that mark the elements that need to be marked according to an error indicator, i.e.

$\eta = \sum_{K\in \mathcal{T}_h(\Omega)} h^2\int_K|\lambda u_h + \Delta u_h|^2\;d\vec{x}+\frac{h}{2}\int_{E\subset \partial K} | [\![ \nabla u\cdot n_E]\!] | ^2\; ds$

In order to do so we begin by computing the value of the indicator using a piecewise constant function space:

def Mark(msh, uh, lam):
W = FunctionSpace(msh, "DG", 0)
# Both the error indicator and the marked element vector will be DG0 field.
w = TestFunction(W)
n = FacetNormal(V.mesh())
h = CellDiameter(msh)
# Assembling the error indicator.
eta = assemble(h**2*R_T**2*w*dx +
(h("+")+h("-"))*(R_dT("+")-R_dT("-"))**2*(w("+")+w("-"))*dS)
frac = .95
delfrac = .05
part = .2
mark = Function(W)
# Filling in the marked element vector using eta.
with mark.dat.vec as markedVec:
with eta.dat.vec as etaVec:
sum_eta = etaVec.sum()
if sum_eta < tolerance:
return markedVec
eta_max = etaVec.max()
sct, etaVec0 = PETSc.Scatter.toZero(etaVec)
markedVec0 = etaVec0.duplicate()
sct(etaVec, etaVec0)
if etaVec.getComm().getRank() == 0:
eta = etaVec0.getArray()
marked = np.zeros(eta.size, dtype='bool')
sum_marked_eta = 0.
#Marking strategy
while sum_marked_eta < part*sum_eta:
new_marked = (~marked) & (eta > frac*eta_max)
sum_marked_eta += sum(eta[new_marked])
marked += new_marked
frac -= delfrac
markedVec0.getArray()[:] = 1.0*marked[:]
sct(markedVec0, markedVec, mode=PETSc.Scatter.Mode.REVERSE)
return mark


It is now time to define the solve, mark and refine loop that is at the heart of the adaptive method described here:

for i in range(max_iterations):
print("level {}".format(i))
lam, uh, V = Solve(msh, labels)
mark = Mark(msh, uh, lam)
msh = msh.refine_marked_elements(mark)


Note that the mesh conforms to the CAD geometry as it is adaptively refined.

## Constructive Solid Geometry in 3D¶

In this section we will focus our attention on three dimensional constructive solid geometry. In particular we will look at the operators +,-,*~, which have been overridden to have a special meaning when applied to two instances of the class CSGeometry. It is important to notice that the same operators can be used also when working with a SplineGeometry and their action will have the same meaning that is presented here. The +,-,* operators have respectively the meaning of union, set difference, and intersection. We will build a cube using the planes intersection and remove from it a portion of sphere:

from netgen.csg import *
left = Plane(Pnt(0, 0, 0), Vec(-1, 0, 0))
right = Plane(Pnt(1, 1, 1), Vec(1, 0, 0))
front = Plane(Pnt(0, 0, 0), Vec(0, -1, 0))
back = Plane(Pnt(1, 1, 1), Vec(0, 1, 0))
bot = Plane(Pnt(0, 0, 0), Vec(0, 0, -1))
top = Plane(Pnt(1, 1, 1), Vec(0, 0, 1))
cube = left * right * front * back * bot * top
cube.bc("cube")
sphere = Sphere(Pnt(0.6, 0.6, 0.6), 0.5)
geo = CSGeometry()
ngmsh = geo.GenerateMesh(maxh=0.1)
msh = Mesh(ngmsh)
File("output/MeshExample3.pvd").write(msh)


Last we will have a look at the Netgen Open Cascade Technology interface, which has been recently included. We will follow the tutorial presented in the NetGen docs, which itself comes from the OCCT tutorial here. The idea is to draw a “flask” using the OCCT interface and solve the linear elasticity equations to compute the stress tensor on the flask subject to gravity. We begin importing the Netgen Open Cascade interface and constructing the bottom of the flask using many different method such as Axes, Face, Pnt, Segment, … (all the details this methods can be found in NetGen docs

from netgen.occ import *
myHeight = 70
myWidth = 50
myThickness = 30
pnt1 = Pnt(-myWidth / 2., 0, 0)
pnt2 = Pnt(-myWidth / 2., -myThickness / 4., 0)
pnt3 = Pnt(0, -myThickness / 2., 0)
pnt4 = Pnt(myWidth / 2., -myThickness / 4., 0)
pnt5 = Pnt(myWidth / 2., 0, 0)
seg1 = Segment(pnt1, pnt2)
arc = ArcOfCircle(pnt2, pnt3, pnt4)
seg2 = Segment(pnt4, pnt5)
wire = Wire([seg1, arc, seg2])
mirrored_wire = wire.Mirror(Axis((0, 0, 0), X))
w = Wire([wire, mirrored_wire])
f = Face(w)
f.bc("bottom")


Once the bottom part of the flask has been constructed we then extrude it to create the main body. We now construct the neck of the flask and fuse it with the main body:

body = f.Extrude(myHeight*Z)
body = body.MakeFillet(body.edges, myThickness / 12.0)
neckax = Axes(body.faces.Max(Z).center, Z)
myNeckHeight = myHeight / 10
body = body + neck
fmax = body.faces.Max(Z)
thickbody = body.MakeThickSolid([fmax], -myThickness / 50, 1.e-3)


Last we are left to construct the threading of the flask neck and fuse it to the rest of the flask body. In order to do this we are going to need the value of pi, which we grab from the Python math package:

import math
cyl1 = Cylinder(neckax, myNeckRadius * 0.99, 1).faces
cyl2 = Cylinder(neckax, myNeckRadius * 1.05, 1).faces
aPnt = Pnt(2 * math.pi, myNeckHeight / 2.0)
aDir = Dir(2 * math.pi, myNeckHeight / 4.0)
aMajor = 2 * math.pi
aMinor = myNeckHeight / 10
arc1 = Ellipse(anAx2d, aMajor, aMinor).Trim(0, math.pi)
arc2 = Ellipse(anAx2d, aMajor, aMinor/4).Trim(0, math.pi)
seg = Segment(arc1.start, arc1.end)
wire1 = Wire([Edge(arc1, cyl1), Edge(seg, cyl1)])
wire2 = Wire([Edge(arc2, cyl2), Edge(seg, cyl2)])

ngmsh = geo.GenerateMesh(maxh=5)