This page was generated from unit-6.1.2-surfacepde/surface_pdes.ipynb.

6.1.2 Surface PDEs

Surface Poisson equation

Homogeneous Dirichlet data

from ngsolve import *
from netgen.csg import *
from netgen.meshing import MeshingStep
from ngsolve.webgui import Draw

Given a half-sphere we want to solve the surface Poisson equation

\[\int_S\nabla_{\Gamma}u\cdot\nabla_{\Gamma}v\,ds = \int_S fv\,ds,\qquad u=0 \text{ on } \Gamma_D\]

geo          = CSGeometry()
sphere       = Sphere(Pnt(0,0,0), 1)
bot          = Plane(Pnt(0,0,0), Vec(0,0,-1))
finitesphere = sphere * bot

geo.AddSurface(sphere, finitesphere.bc("surface"))
geo.NameEdge(sphere,bot, "bottom")

mesh = Mesh(geo.GenerateMesh(maxh=0.3))
mesh.Curve(2)
Draw(mesh)

BaseWebGuiScene

Therefore, we define the usual \(H^1\) finite element space and use the dirichlet_bbnd flag to indicate the BBoundary on which the Dirichlet conditions are prescribed. The test- and trial-functions are given as usual.

fes = H1(mesh, order=2, dirichlet_bbnd="bottom")
u, v = fes.TnT()
print(fes.FreeDofs())

0: 00000000000000000000011111111111111111111111111111
50: 11111111111111111111111111111111111111111111111100
100: 11011011011011011011011011011011011011011011011011
150: 01101101111111111111111111111111111111111111111111
200: 11111111111111111111111111111111111111111111111111
250: 11111111111111111111111111111111111111111111111111
300: 11111111111111111111111111111111111111111111111111
350: 111111111111111111

For the (bi-)linear form we have to take care that we have to define the according integrators on the boundary and that the Trace operator has to be used to obtain the tangential/surface derivative

a = BilinearForm(fes, symmetric=True)
a += grad(u).Trace()*grad(v).Trace()*ds
a.Assemble()

force = sin(x)*y*exp(z)

f = LinearForm(fes)
f += force*v*ds
f.Assemble()

<ngsolve.comp.LinearForm at 0x7ff14469e870>

Solving is done as usual

gfu = GridFunction(fes)
gfu.vec.data = a.mat.Inverse(fes.FreeDofs())*f.vec
Draw(gfu, mesh, "u")

BaseWebGuiScene

Inhomogeneous Dirichlet data

To solve the same problem with non-homogenous Dirichlet data, \(u=u_D\) on \(\Gamma_D\) the same technique as in the volume case is used, where we have to set a function on the BBoundary instead on the boundary

gfu.Set(x, definedon=mesh.BBoundaries("bottom"))
r = f.vec.CreateVector()
r.data = f.vec - a.mat*gfu.vec
gfu.vec.data += a.mat.Inverse(fes.FreeDofs())*r
Draw(gfu)

BaseWebGuiScene

Finite element spaces for surfaces

We differ between two types of finite element spaces for surfaces. The first class consists of spaces, where the restriction of a 3D element to the surface leads to a valid 2D element of the same type: H1, HCurl, HCurlCurl, NumberSpace.

These spaces can directly be used, one has to take care using the Trace operator. Otherwise an exception is thrown during assembling.

a += grad(u)*grad(v)*ds
a.Assemble()

---------------------------------------------------------------------------
NgException                               Traceback (most recent call last)
Cell In[7], line 1
----> 1 a += grad(u)*grad(v)*ds
      2 a.Assemble()

NgException: Trialfunction does not support BND-forms, maybe a Trace() operator is missing, type = grad

The NumberSpace is an exception as it represents only a number, where no Trace operator has to be used.

The second class is given by

Space	Surface Space
L2	SurfaceL2
HDiv	HDivSurface
HDivDiv	HDivDivSurface
FacetFESpace	FacetSurface

Here, a 2D reference element is mapped directly onto the surface. To be consistent, also here the Trace operator has to be used.