Surface meshes

In the present section we will show how to generate a surface mesh using Open Cascade and how to solve the heat equation on it. We will start by generating a surface mesh for a torus.

from netgen.occ import *
from netgen.meshing import MeshingStep
from math import sin, cos, pi
import netgen.gui

def Curve(t): return Pnt(0, 3+1.5*cos(t), sin(t))
n = 100
pnts = [Curve(2*pi*t/n) for t in range(n+1)]

spline = SplineApproximation(pnts)
f = Face(Wire(spline))

torus = f.Revolve(Axis((0,0,0), Z), 360)

ngmesh = OCCGeometry(torus).GenerateMesh(maxh=0.3, perfstepsend=MeshingStep.MESHSURFACE)

Once we have generated a surface mesh using Netgen, we can simply pass it to Firedrake and using the curve_field method to generate a curved surface mesh.

from firedrake import *
mesh = Mesh(Mesh(ngmesh).curve_field(3))

We begin defining the initial condition for the problem, in particular, we will consider a Gaussian bump centered at a point on the torus.

x, y, z = SpatialCoordinate(mesh)
V = FunctionSpace(mesh, "CG", 3)
u0 = assemble(interpolate(exp(-10*((x-0.400126)**2 + (y+4.22785)**2 + (z-0.550597)**2)),V))
u = Function(V)
u.interpolate(u0)
out = VTKFile("output/heat.pvd")
out.write(u, time=0.0)

We now use Irksome, to describe the time-dependent problem we aim to solve, i.e.

\[\begin{split}\begin{align*} \frac{\partial u}{\partial t} - \nabla \cdot \nabla u &= 0 \quad \text{in} \quad \Omega \times (0, T], \\ u &= u_0 \quad \text{at} \quad t = 0. \end{align*}\end{split}\]
from irksome import Dt
v = TestFunction(V)
F = inner(Dt(u), v)*dx + inner(grad(u), grad(v))*dx

We now need to define the final time we want to reach, and the number of time steps we want to take and the solver parameters. We will then construct a TimeStepper that will take care of time-stepping the PDE.

from irksome import GaussLegendre, Dt, MeshConstant, TimeStepper
MC = MeshConstant(mesh)
T = 1.0
N = 100
dt = MC.Constant(T / N)
t = MC.Constant(0.0)
luparams = {"mat_type": "aij",
        "ksp_type": "preonly",
        "pc_type": "lu"}
butcher_tableau = GaussLegendre(1)
ns = butcher_tableau.num_stages
stepper = TimeStepper(F, butcher_tableau, t, dt, u,
                  solver_parameters=luparams)
while (float(t) < T):
    if (float(t) + float(dt) > T):
        dt.assign(T - float(t))
    stepper.advance()
    print(float(t))
    t.assign(float(t) + float(dt))

    out.write(u, time=float(t))