Wave equation

from ngsolve import *
from ngstrefftz import *
from netgen.occ import *
from ngsolve.TensorProductTools import MakeTensorProductMesh, SegMesh
from ngsolve import *
from ngsolve.webgui import Draw

We aim to solve the first order wave equation

\[\begin{split}\begin{align} \nabla\cdot\sigma + c^{-2}\partial_t v=0 \\ \nabla v + \partial_t \sigma =0\\ v=v_0,\ \sigma=\sigma_0 \text{ at } t=0\\ v=g_D \text{ on }\Gamma_D \end{align}\end{split}\]

Create tensor-product mesh

A space-time mesh with time along the z-Axis.

meshx = Mesh(unit_square.GenerateMesh(maxh=0.2))
mesht = Mesh(SegMesh(5,0,1,periodic=False) )
mesh = Mesh(MakeTensorProductMesh(meshx,mesht))

Trefftz space

The Trefftz-space for an element \(K\) in the mesh is given by

\[\mathbb T^p(K)=\left\{v\in \mathbb P^p \text{ such that } (c^{-2}\partial_{tt} -\Delta) u=0 \text{ on } K\right\}\]

c=2
order=4
fes = trefftzfespace(mesh,order=order,dgjumps=True,eq="fowave")
fes.SetCoeff(c)

Test case

from numpy import pi
sq = sqrt(2.0)
sig_ex = CoefficientFunction((-cos(pi*x)*sin(pi*y)*sin(pi*z*c*sq)/sq,
    -sin(pi*x)*cos(pi*y)*sin(pi*z*c*sq)/sq))
v_ex = sin(pi*x)*sin(pi*y)*cos(pi*z*c*sq)*c
gD = v_ex

Trefftz-DG formulation

\[\begin{split}\begin{align} \text{Seek}&\quad (v_{hp},\sigma_{hp})\in V_p (\mathcal{T}_h)\nonumber\\ \text{such that}&\quad \mathcal{A}(v_{hp},\sigma_{hp}; w ,\tau )=\ell( w ,\tau )\quad \forall ( w ,\tau )\in V_p (\mathcal{T}_h), \nonumber\\ \text{where}&\nonumber\\ \mathcal{A}(v_{hp},&\sigma_{hp}; w ,\tau ):= %-\sum_{K\in\mathcal{T}_h} \int_K\bigg(v\Big(\nabla\cdot\tau+c^{-2}\partial_t w \Big) +\sigma_{hp}\cdot\Big(\partial_t \tau +\nabla w \Big)\bigg) %\nonumber\\ +& \int_{F^{space}}\big(c^{-2}v_{hp}^-[[{w}]]_t+\sigma_{hp}^-\cdot[[{\tau}]]_t+v_{hp}^-[[{\tau}]]_N+\sigma_{hp}^-\cdot[[{w}]]_N\big) \nonumber \\ +&\int_{F^{time}}\!\! \big( \{\{{v_{hp}}\}\}[[{\tau }]]_N+\{\{{\sigma_{hp}}\}\}\cdot[[{ w }]]_N +\alpha[[{v_{hp}}]]_N\cdot[[{ w }]]_N+ \beta[[{\sigma_{hp}}]]_N[[{\tau }]]_N \big) \nonumber\\ +&\int_{F^T} (c^{-2}v_{hp} w +\sigma_{hp} \cdot\tau ) +\int_{F^D} \big(\sigma\cdot n_\Omega^x\, w +\alpha v_{hp} w \big) \nonumber\\ \ell( w ,&\tau ):= \int_{F^O} ( c^{-2}v_0 w +\sigma_0\cdot \tau ) +\int_{F^D} g_D\big(\alpha w -\tau\cdot n_\Omega^x\big) \nonumber \end{align}\end{split}\]

D = fes.mesh.dim - 1

vsig = fes.TrialFunction()
wtau = fes.TestFunction()
vsigo=vsig.Other()
wtauo=wtau.Other()

v = vsig[D]
sig = CoefficientFunction(tuple([vsig[i] for i in  range(D)]))
w = wtau[D]
tau = CoefficientFunction(tuple([wtau[i] for i in  range(D)]))

vo = vsigo[D]
sigo = CoefficientFunction(tuple([vsigo[i] for i in  range(D)]))
wo = wtauo[D]
tauo = CoefficientFunction(tuple([wtauo[i] for i in  range(D)]))

h = specialcf.mesh_size
n = specialcf.normal(D+1)
n_t = n[D]/Norm(n)
n_x = CoefficientFunction( tuple([n[i]/Norm(n) for i in  range(D)]) )

mean_v = 0.5*(v+vo)
mean_w = 0.5*(w+wo)
mean_sig = 0.5*(sig+sigo)
mean_tau = 0.5*(tau+tauo)

jump_vx = ( v - vo ) * n_x
jump_wx = ( w - wo ) * n_x
jump_sigx = (( sig - sigo ) * n_x)
jump_taux = (( tau - tauo ) * n_x)

jump_vt = ( v - vo ) * n_t
jump_wt = ( w - wo ) * n_t
jump_sigt = ( sig - sigo ) * n_t
jump_taut = ( tau - tauo ) * n_t

timelike = n_x*n_x # n_t=0
spacelike = n_t**2 # n_x=0

#params
alpha=0.5
beta=0.5

a = BilinearForm(fes)
#space like faces, w/o x jump ASSUME TENSOR MESH
a += spacelike * ( pow(c,-2)*IfPos(n_t,v,vo)*jump_wt + IfPos(n_t,sig,sigo)*jump_taut) * dx(skeleton=True)
#time like faces
a += timelike * ( mean_v*jump_taux + mean_sig*jump_wx + alpha*jump_vx*jump_wx + beta*jump_sigx*jump_taux ) * dx(skeleton=True )        #t=T (or *x)
a += ( pow(c,-2)*v*w + sig*tau) *  IfPos(n_t,1,0)*spacelike * ds(skeleton=True)
#dirichlet boundary 'timelike'
a += ( sig*n_x*w + alpha*v*w ) * timelike *   ds(skeleton=True)

f = LinearForm(fes)
f += ( pow(c,-2)*v_ex*w + sig_ex*tau ) * IfPos(n_t,0,1)*spacelike *  ds(skeleton=True)
f += ( gD * (alpha*w - tau*n_x) ) * timelike *   ds(skeleton=True)

with TaskManager():
    a.Assemble()
    f.Assemble()

gfu = GridFunction(fes, name="uDG")
gfu.vec.data = a.mat.Inverse()*f.vec
Draw(gfu,clipping={"z":-1, "dist":-.11})

BaseWebGuiScene

Use clipping to show the solution inside the space-time mesh.

\(L^2\)-error

Integrate((gfu-CoefficientFunction((sig_ex,v_ex)))**2,mesh)

4.0430535968112005e-07