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

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)

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)