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.

[1]:

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

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[1], line 3
      1 from ngsolve import *
      2 from ngsolve import dx as ngsdx
----> 3 from xfem import *
      4 from xfem.lsetcurv import *
      5 from ngstrefftz import *

ModuleNotFoundError: No module named 'xfem'

Level set

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

[2]:

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

[3]:

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)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[3], line 1
----> 1 lsetmeshadap = LevelSetMeshAdaptation(mesh, order=order, threshold=0.1,
      2                                       discontinuous_qn=True,heapsize=int(1e8))
      3 deformation = lsetmeshadap.CalcDeformation(levelset)
      4 mesh.SetDeformation(deformation)

NameError: name 'LevelSetMeshAdaptation' is not defined

Trefftz FE space

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

[4]:

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)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[4], line 1
----> 1 ci = CutInfo(mesh, lsetp1)
      2 els_hasneg = ci.GetElementsOfType(HASNEG)
      3 els_if = ci.GetElementsOfType(IF)

NameError: name 'CutInfo' is not defined

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

[5]:

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

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[5], line 3
      1 lambda_dg = lambda_nitsche = 10
      2 h = specialcf.mesh_size
----> 3 nh = 1.0 / Norm(grad(lsetp1)) * grad(lsetp1)
      4 nF = specialcf.normal(mesh.dim)
      5 flux_u = -0.5 * (grad(u) + grad(u.Other())) * nF

NameError: name 'lsetp1' is not defined

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*}\]

[6]:

# 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

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[6], line 2
      1 # Facets for ghost penalty stabilization
----> 2 facets_gp = GetFacetsWithNeighborTypes(mesh, a=els_hasneg, b=els_if)
      3 gamma_stab = 10
      4 dw = dFacetPatch(definedonelements=facets_gp)

NameError: name 'GetFacetsWithNeighborTypes' is not defined

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

[7]:

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)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[7], line 4
      1 exact = exp(x)*sin(y)
      2 coef_g = exact
----> 4 f = LinearForm(Vh)
      5 f += (- grad(v) * nh * coef_g + lambda_nitsche / h * coef_g * v) * ds
      7 f.Assemble()

NameError: name 'Vh' is not defined

[8]:

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

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[8], line 1
----> 1 DrawDC(lsetp1, gfu, -1, mesh, "u",min=0,max=1.5)

NameError: name 'DrawDC' is not defined