This page was generated from unit-3.6-surfacehdg/surfacehdg.ipynb.

3.6 Scalar PDE on surfaces (with HDG)

We are solving the scalar linear transport problem

$${\partial }_{t}u-\epsilon {\mathrm{\Delta }}_{\mathrm{\Gamma }}u+{\mathrm{div}}_{\mathrm{\Gamma }}(\mathbf{w}u)=0\text{ on }\mathrm{\Gamma },$$

where $\mathrm{\Gamma }$ is a closed oriented surface in ${\mathbb{R}}^3$.

from ngsolve import *
from ngsolve.webgui import Draw

Mesh of (only) the surface:

We consider the unit sphere (only surface mesh!)

from netgen.occ import *
sphere = Sphere((0,0,0),1).faces[0]
sphere.name = "sphere"
geo = OCCGeometry(sphere)
mesh = Mesh(geo.GenerateMesh(maxh=0.5))
Draw(mesh)

BaseWebGuiScene

order = 4
if order > 0:
    mesh.Curve(order)
Draw(mesh)

BaseWebGuiScene

Tangential velocity field:

t = Parameter(0.)
b = CoefficientFunction((y,-x,0))
Draw (b, mesh, "b")
eps = 5e-5

some local quantities

normal, tangential, co-normal and mesh size

n = specialcf.normal(mesh.dim)
h = specialcf.mesh_size
local_tang = specialcf.tangential(mesh.dim)
con = Cross(n,local_tang) #co-normal pointing outside of a surface element
bn = InnerProduct(b,con)

Discretization

Surface FESpaces:

fesl2 = SurfaceL2(mesh, order=order)
fesedge = FacetSurface(mesh, order=order)

fes = fesl2*fesedge
u,uE = fes.TrialFunction()
v,vE = fes.TestFunction()

surface gradients, co-normal derivatives and jumps:

gradu = u.Trace().Deriv()
gradv = v.Trace().Deriv()
jumpu = u - uE
jumpv = v - vE
dudn = InnerProduct(gradu,con)
dvdn = InnerProduct(gradv,con)

Diffusion formulation (Hybrid Interior Penalty on Surface):

$$\sum _T{\int }_T\mathrm{\nabla }u\mathrm{\nabla }v-\sum _T{\int }_{\partial T}n\mathrm{\nabla }u(v-\hat{v})-\sum _T{\int }_{\partial T}n\mathrm{\nabla }u(u-\hat{u})+\frac{\alpha p^2}{h}\sum _F{\int }_F(u-\hat{u})(v-\hat{v})$$

a = BilinearForm(fes)
# diffusion

dr = ds(element_boundary=True)

alpha = 10*(order+1)**2
a += eps*gradu * gradv * ds
a += - eps*dudn * jumpv * dr
a += - eps*dvdn * jumpu * dr
a += eps*alpha/h * jumpu * jumpv * dr

Convection formulation (Hybrid Upwinding):

$$\sum _T{\int }_T-u\mathbf{w}\cdot \mathrm{\nabla }v+\sum _T{\int }_{\partial T}(\mathbf{w}\cdot \mathbf{n})u^{hyb-upw}v+\sum _T{\int }_{\partial T_{out}}(\mathbf{w}\cdot \mathbf{n})(\hat{u}-u)\hat{v}$$

with

$$\begin{array}{r}{u}^{hyb-upw}=\{\begin{array}{cll}u& \text{if }\mathbf{w}\cdot \mathbf{n}& (outflow),\\ \hat{u}& \text{else}& (inflow).\end{array}\end{array}$$

outflow	inflow

a = BilinearForm(fes)
# convection (surface version of Egger+Schöberl formulation):
a += - u * InnerProduct(b,gradv) * ds
a += IfPos(bn, bn*u, bn*uE) * v * dr + IfPos(bn, bn, 0) * (uE - u) * vE * dr
a.Assemble()

<ngsolve.comp.BilinearForm at 0x7f95ba36dc30>

Mass matrix:

m = BilinearForm(fes)
m += u*v*ds
m.Assemble()

<ngsolve.comp.BilinearForm at 0x7f95b9938430>

$M^{\ast }$-matrix for implicit Euler

dt = 0.02

mstar = m.mat.CreateMatrix()
mstar.AsVector().data = m.mat.AsVector() + dt * a.mat.AsVector()
mstarinv = mstar.Inverse()

from math import pi
T = 8*pi

Time stepping

Visualization of solution (Deformation in normal direction)

gfu = GridFunction(fes)
gfu.components[0].Set(1.5*exp(-20*(x*x+(y-1)**2+(z)**2)),definedon=mesh.Boundaries("sphere"))

The time loop

sceneun = Draw (gfu.components[0]*n, mesh, "un", deformation = True, autoscale=False)
t.Set(0)
for i in range(int(T/dt)):
    gfu.vec.data = mstarinv @ m.mat * gfu.vec
    sceneun.Redraw()
    t.Set(t.Get()+dt)

Remark:

• We can apply static condensation as in the volume