Unfitted Trefftz-DG method

We are looking to solve

\[\begin{split}\begin{align*} \begin{split} \begin{cases} -\Delta u = 0 &\text{ in } \Omega, \\ u=g &\text{ on } \partial \Omega, \end{cases} \end{split} \end{align*}\end{split}\]

using a Trefftz-DG formulation on an unfitted domain. We use ngsxfem for the unfitted finite element discretization.

from ngsolve import *
from ngsolve import dx as ngsdx
from xfem import *
from xfem.lsetcurv import *
from ngstrefftz import *
from netgen.occ import *
order=4

importing ngsxfem-2.1.dev

Level set

We setup the background mesh and describe the geometry using a levelset function

levelset = (sqrt(((x-0.5)*3)**2+((y-0.5)*3)**2)) - 1.0
mesh = Mesh(unit_square.GenerateMesh(maxh=.1))

lsetmeshadap = LevelSetMeshAdaptation(mesh, order=order, threshold=0.1,
                                      discontinuous_qn=True,heapsize=int(1e8))
deformation = lsetmeshadap.CalcDeformation(levelset)
mesh.SetDeformation(deformation)
lsetp1 = lsetmeshadap.lset_p1
DrawDC(lsetp1, 1, 0, mesh, "levelset", min=0)

BaseWebGuiScene

Trefftz FE space

We use the Trefftz space of harmonic polynomials and \(\texttt{Restrict}\) to remove unused dofs

ci = CutInfo(mesh, lsetp1)
els_hasneg = ci.GetElementsOfType(HASNEG)
els_if = ci.GetElementsOfType(IF)
# Facets with parts inside the domain
facets_dg = GetFacetsWithNeighborTypes(mesh, a=els_hasneg, b=els_hasneg)

Vhbase = trefftzfespace(mesh, order=order, eq='laplace', dgjumps=True)
Vh = Restrict(Vhbase, els_hasneg)
u, v = Vh.TnT()
gfu = GridFunction(Vh)

DG formulation

\(\newcommand{\nF}{\nb_F} \newcommand{\dnF}{\nF\cdot\nabla} \newcommand{\dn}{\nb\cdot\nabla} \newcommand{\jump}[1]{[\![#1]\!]} \newcommand{\mean}[1]{\{\!\!\{#1\}\!\!\}} \newcommand{\Thact}{\mathcal{T}_h^{\text{act}}} \newcommand{\Fhact}{\mathcal{F}_h^{\text{act}}} \newcommand{\ThG}{\Th^{\Gamma}} \newcommand{\Fhdg}{\Fh(\Th)} \newcommand{\Fhdgact}{\Fh(\Thact)} \newcommand{\Fhgp}{\Fh^\text{gp}} \newcommand{\Fhgpmin}{\Fh^{\text{gp},\min}} \newcommand{\Fhgpstar}{\Fh^{\text{gp}\star}} \newcommand{\xb}{\bm{x}} \newcommand{\xbh}{\hat{\xb}} \newcommand{\nb}{\mathbf{n}} \newcommand{\nbF}{\nb_F} \newcommand{\nbFh}{\hat{\nb}_F} \newcommand{\Ex}{\mathcal{E}} \newcommand{\inter}{\pi_h} \newcommand{\interE}{\pi_h^e} \newcommand{\PP}{\mathbb{P}} \newcommand{\Th}{\mathcal{T}_h} \newcommand{\Sh}{\mathcal{S}_h} \newcommand{\Fh}{\mathcal{F}_h} \newcommand{\Sh}{\mathcal{S}_h} \DeclareMathOperator{\Int}{Int} \newcommand{\Ah}{\mathcal{A}_h}\) We introduce the bilinear form \(\Ah\) given by

\[\begin{split}\begin{equation} \begin{aligned} \Ah(u,v) := &\ (\nabla u, \nabla v)_\Omega - (\dn u, v)_\Gamma - (u, \dn v)_\Gamma + \beta (h^{-1} u,v)_\Gamma \\ &- (\mean{\dnF u }, \jump{ v })_{\Fh \cap \Omega} - (\jump{u}, \mean{\dnF v } )_{\Fh \cap \Omega} + \beta (h^{-1} \jump{u}, \jump{v})_{\Fh \cap \Omega}, \end{aligned} \end{equation}\end{split}\]

and setup the necessary differential symbols for cut elements

lambda_dg = lambda_nitsche = 10
h = specialcf.mesh_size
nh = 1.0 / Norm(grad(lsetp1)) * grad(lsetp1)
nF = specialcf.normal(mesh.dim)
flux_u = -0.5 * (grad(u) + grad(u.Other())) * nF
flux_v = -0.5 * (grad(v) + grad(v.Other())) * nF
jump_u = u - u.Other()
jump_v = v - v.Other()

dX = ngsdx(definedonelements=els_hasneg)
dx = dCut(lsetp1, NEG, definedonelements=els_hasneg)
dk = dCut(lsetp1, NEG, skeleton=True, definedonelements=facets_dg)
ds = dCut(lsetp1, IF, definedonelements=els_if)

a = RestrictedBilinearForm(Vh, element_restriction=els_hasneg,
                           facet_restriction=facets_dg,
                           check_unused=True)
a += (grad(u) * grad(v)) * dx
a += (lambda_dg / h * jump_u * jump_v + flux_u * jump_v + flux_v * jump_u) * dk
a += (-grad(u) * nh * v - grad(v) * nh * u + lambda_nitsche / h * u * v) * ds

To deal with shape-irregular cuts we add a ghost penalty stabilization term. Different realisations are possible, here we will implement

\[\begin{align}\label{eq:sh1} \Sh(u,v) = \sum_{F \in \Fhgp} \gamma ~ ({h_F^{-2}} \jump{ \Pi_{\Th} u}_{\omega_F}, \jump{\Pi_{\Th} v}_{\omega_F})_{\omega_F}. \end{align}\]

where \(\gamma > 0\) is a corresponding stabilisation parameter.

Here, \(\Pi_{\Th}\) denotes the element-wise \(L^2\) projection onto \(\PP^k(\Th)\), \(\omega_F = \Int( \overline{T}_1 \cup \overline{T}_2\)) denotes the element aggregation to a facet \(F \in \Fhgp\), \(F= \partial T_1 \cap \partial T_2\). The facet (volumetric) patch jump \(\jump{u}_{\omega_F}\) of a polynomial \(u_h\in\PP^k(\Th)\) is defined as

\[\begin{equation*} \jump{u_h}_{\omega_F}|_{T_i} = u_h|_{T_i} - \Ex^P(u_h|_{T_j}),\quad\text{for }i,j\in\{1,2\}\text{ and } i\neq j, \end{equation*}\]

where \(\Ex^P\) denotes the canonical extension of a polynomial from \(T\) to \(\tilde \Omega\), i.e.

\[\begin{equation*} \Ex^P:\PP^k(T)\rightarrow \PP^k(\tilde\Omega)\ \text{such that}\ \Ex^Pv|_T=v \ \text{for} \ v \in \PP^k(T). \end{equation*}\]

# Facets for ghost penalty stabilization
facets_gp = GetFacetsWithNeighborTypes(mesh, a=els_hasneg, b=els_if)
gamma_stab = 10
dw = dFacetPatch(definedonelements=facets_gp)

a += (gamma_stab / h**2 * (u - u.Other()) * (v - v.Other())) * dw

We setup the right hand side for the exact solution \(u=\exp(x) \sin(y)\)

exact = exp(x)*sin(y)
coef_g = exact

f = LinearForm(Vh)
f += (- grad(v) * nh * coef_g + lambda_nitsche / h * coef_g * v) * ds

f.Assemble()
a.Assemble()

nd = Vh.ndof
gfu.vec.data = a.mat.Inverse(Vh.FreeDofs(), 'sparsecholesky') * f.vec

l2error = sqrt(Integrate((gfu - exact)**2 * dx, mesh))
print(l2error)
print(mesh.ne)

1.3179411097355532e-08
226

DrawDC(lsetp1, gfu, -1, mesh, "u",min=0,max=1.5)

/home/runner/.local/lib/python3.10/site-packages/netgen/webgui.py:17: RuntimeWarning: overflow encountered in cast
  values = np.array(data.flatten(), dtype=dtype)

BaseWebGuiScene