# 3.5 Scalar PDE on surfaces (with HDG)¶

We are solving the scalar linear transport problem

$\partial_t u - \varepsilon \Delta_\Gamma u + \operatorname{div}_\Gamma ( \mathbf{w} u ) = 0 \quad \text{ on } \Gamma,$

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

[1]:

from ngsolve import *
from ngsolve.webgui import Draw
from math import pi


## Mesh of (only) the surface:¶

We consider the unit sphere (only surface mesh!)

[2]:

from netgen.csg import *
geo = CSGeometry()

[2]:

0

[3]:

from netgen.meshing import MeshingParameters
from netgen.meshing import MeshingStep
mp = MeshingParameters(maxh=0.5, perfstepsend = MeshingStep.MESHSURFACE)

mesh = Mesh(geo.GenerateMesh(mp=mp))

[4]:

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

[4]:

BaseWebGuiScene


### Tangential velocity field:¶

[5]:

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

[6]:

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:¶

[7]:

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:¶

[8]:

gradu = u.Trace().Deriv()
jumpu = u - uE
jumpv = v - vE


### Diffusion formulation (Hybrid Interior Penalty on Surface):¶

$\sum_T \int_T \nabla u \nabla v - \sum_T \int_{\partial T} n \nabla u (v-\widehat v) - \sum_T \int_{\partial T} n \nabla u (u-\widehat u) + \frac{\alpha p^2}{h} \sum_F \int_F (u-\widehat u)(v-\widehat v)$
[9]:

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 \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{split}u^{hyb-upw} = \left\{ \begin{array}{cll} u & \text{if } \mathbf{w} \cdot \mathbf{n} & (outflow), \\ \hat{u} & \text{else} & (inflow) . \end{array} \right.\end{split}$

outflow

inflow

[10]:

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

[10]:

<ngsolve.comp.BilinearForm at 0x7facfed21430>


### Mass matrix:¶

[11]:

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

[11]:

<ngsolve.comp.BilinearForm at 0x7facfed22570>


### $$M^\ast$$-matrix for implicit Euler¶

[12]:

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

[13]:

gfu = GridFunction(fes)
gfu.components[0].Set(1.5*exp(-20*(x*x+(y-1)**2+(z)**2)),definedon=mesh.Boundaries("sphere"))
sceneu = Draw (gfu.components[0], mesh, "u")
sceneun = Draw (gfu.components[0]*n, mesh, "un", order=3, autoscale=False)
res = gfu.vec.CreateVector()
SetVisualization(deformation=True)


### The time loop¶

[14]:

t.Set(0)
for i in range(int(T/dt)):
res.data = m.mat * gfu.vec
gfu.vec.data = mstarinv * res
sceneu.Redraw()
sceneun.Redraw()
t.Set(t.Get()+dt)


### Remark:¶

• We can apply static condensation as in the volume