High-order meshes

This tutorial is an improved version of the Netgen tutorial available in the Firedrake documentation, which can be found here. It is possible to construct high-order meshes for a geometry constructed in Netgen. In order to do so we need to use the curve_field method of a Firedrake Mesh object generated from a Netgen mesh. In particular, we need to pass the degree of the polynomial field we want to use to parametrise the coordinates of the domain to the curve_field method, which will return a Function constructed on a DG space for this purpose.

from firedrake import *
from netgen.occ import WorkPlane, OCCGeometry
import netgen
from mpi4py import MPI

wp = WorkPlane()
if COMM_WORLD.rank == 0:
    for i in range(6):
        wp.Line(0.6).Arc(0.4, 60)
    shape = wp.Face()
    ngmesh = OCCGeometry(shape,dim=2).GenerateMesh(maxh=1.)
else:
    ngmesh = netgen.libngpy._meshing.Mesh(2)
mesh = Mesh(Mesh(ngmesh, comm=COMM_WORLD).curve_field(4))
VTKFile("output/MeshExample5.pvd").write(mesh)
Example of a curved mesh of order 4 generated from a geometry described using Open Cascade WorkPlane.

High-order meshes are also supported in three dimensions; we just need to specify the correct dimension when constructing the OCCGeometry object. We will now show how to solve the Poisson problem on a high-order mesh, of order 3, for the unit sphere.

from netgen.occ import Sphere, Pnt
import netgen
from mpi4py import MPI

if COMM_WORLD.rank == 0:
    shape = Sphere(Pnt(0,0,0), 1)
    ngmesh = OCCGeometry(shape,dim=3).GenerateMesh(maxh=1.)
else:
    ngmesh = netgen.libngpy._meshing.Mesh(3)
mesh = Mesh(Mesh(ngmesh).curve_field(3))
# Solving the Poisson problem
VTKFile("output/MeshExample6.pvd").write(mesh)
x, y, z = SpatialCoordinate(mesh)
V = FunctionSpace(mesh, "CG", 3)
f = Function(V).interpolate(1+0*x)
u = TrialFunction(V)
v = TestFunction(V)
a = inner(grad(u), grad(v)) * dx
l = inner(f, v) * dx

sol = Function(V)

bc = DirichletBC(V, 0.0, [1])
A = assemble(a, bcs=bc)
b = assemble(l)
bc.apply(b)
solve(A, sol, b, solver_parameters={"ksp_type": "cg", "pc_type": "lu"})

VTKFile("output/Sphere.pvd").write(sol)
The solution of the Poisson problem solved on a mesh of order 3 for the unit sphere.

It is also possible to construct high-order meshes using the SplineGeometry, CSG2d and CSG classes.

from netgen.geom2d import CSG2d, Circle, Rectangle
import netgen
from mpi4py import MPI

if COMM_WORLD.rank == 0:
   geo = CSG2d()
   circle = Circle(center=(1,1), radius=0.1, bc="curve").Maxh(0.01)
   rect = Rectangle(pmin=(0,1), pmax=(1,2),
                    bottom="b", left="l", top="t", right="r")
   geo.Add(rect-circle)
   ngmesh = geo.GenerateMesh(maxh=0.2)
else:
    ngmesh = netgen.libngpy._meshing.Mesh(2)
mesh = Mesh(Mesh(ngmesh,comm=COMM_WORLD).curve_field(2))
VTKFile("output/MeshExample7.pvd").write(mesh)
Example of a curved mesh of order 2 generated from a geometry described using Netgen CSG2d.