This page was generated from unit-8.9-unfmixed/unfmixed.ipynb.

8.9 Unfitted Mixed FEM (using HDiv FE spaces)

Our goal is to consider the mixed Poisson problem with Dirichlet boundary conditions in an unfitted setting: Find \(u: \Omega \rightarrow \mathbb{R}^d\), \(p : \Omega \rightarrow \mathbb{R}\) such that

\begin{alignat}{2} u - \nabla p &= \phantom{-} 0 \quad &&\text{in } \Omega, \tag{C1}\\ \operatorname{div} u &= -f \quad &&\text{in } \Omega, \tag{C2} \\ p &= \phantom{-} p_D \quad &&\text{on } \partial \Omega. \tag{C3} \end{alignat}

Here \(\Omega \subset \mathbb{R}^d\) is a bounded domain that is not parametrized by the computational mesh, but contained in a background domain \(\widetilde{\Omega}\) such that \(\overline \Omega \subseteq \widetilde{\Omega}\) which is triangulated by a simplicial, shape regular and quasi-uniform triangulation \(\widetilde{\mathcal{T}_h}\). We set \begin{align} \mathcal{T}_h &= \big\{T \in \widetilde{\mathcal{T}_h} : \text{meas}(T \cap \Omega) > 0 \big\}, \tag{T1} \\ \mathcal{T}_h^{\text{cut}} &= \big\{T \in \widetilde{\mathcal{T}_h} : \text{meas}(T \cap \partial \Omega) > 0 \big\}, \tag{T2}\\ \mathcal{T}_h^{\text{int}} &= \mathcal{T}_h \setminus \mathcal{T}_h^{\text{cut}}. \tag{T3} \end{align}

The Method

To obtain inf-sup stable unfitted discretizations, one usually add stabilization terms like the ghost penalty stabilization. In the case of mixed problems, however, adding such terms on the \(p\)-\(q\)-coupling block would pollute the mass balance. Instead, one can either use a ghost penalty on the \(\operatorname{div}\)- and \(\nabla\)-block as done e.g. in this preprint(arXiv:2205.12023) or change the domain on which the \(\operatorname{div}\)- and \(\nabla\)-blocks act as we have done in this preprint(arXiv:2306.12722). We showed that one can consider the following weak formulation: Find \((u_h,\bar p_h) \in \mathbb{RT}^k(\mathcal{T}_h) \times \mathbb{P}^k(\mathcal{T}_h)\) such that \begin{alignat}{2} a_h(u_h,v_h) + b_h(v_h,\bar p_h) &= (v_h \cdot n, p_D)_{\partial \Omega} \qquad &&\forall v_h \in \mathbb{RT}^k(\mathcal{T}_h), \tag{M1} \\ b_h(u_h,q_h) &= - (f_h,q_h)_{\mathcal{T}_h} \qquad &&\forall q_h \in \mathbb{P}^k(\mathcal{T}_h), \tag{M2} \end{alignat} where \begin{align} a_h(u_h,v_h) &:= (u_h,v_h)_{\Omega} + \gamma_u j_h (u_h,v_h), \tag{A}\\ b_h(v_h,p_h) &:= (\operatorname{div} v_h,p_h)_{\mathcal{T}_h}. \tag{B} \end{align}

The following code imports the basic functionality of netgen, ngsolve and ngsxfem. Afterwards, we set some general parameters for this notebook.

#basic imports
from netgen.geom2d import SplineGeometry
from ngsolve import *
from ngsolve.internal import *
from xfem import *
from xfem.lsetcurv import *

importing ngsxfem-2.1.2303.dev0

We consider a ring geometry that is contained inside the square \([-1,1]^2\) with inner radius \(r_1 = 0.25\) and outer radius \(r_2 = 0.75\). Then, the geometry is described by the following signed distance function \begin{equation} \Phi(x) = \Big\vert \sqrt{x_1^2 + x_2^2} - \frac{1}{2}(r_1+r_2) \Big\vert - \frac{r_2-r_1}{2}. \end{equation}

square = SplineGeometry()
square.AddRectangle((-1, -1), (1, 1), bc=1)
mesh = Mesh(square.GenerateMesh(maxh=0.2))
r1, r2 = 1/4, 3/4  # inner/ outer radius
rr, rc = (r2 - r1) / 2.0 , (r1 + r2) / 2.0
r = sqrt(x**2 + y**2)
levelset = IfPos(r - rc, r - rc - rr, rc - r - rr)
DrawDC(levelset,-1.0,2.0,mesh,"x")

BaseWebGuiScene

Discretization parameters

We choose discretization order and ghost penalty stabilization parameter (0 for no stabilization):

order = 2
gamma_stab = 0 # 0 if no ghost penalty is to be used

Geometry approximation

Next, we construct the mesh deformation for higher order accuracy (see the intlset-tutorial for further details):

lsetmeshadap = LevelSetMeshAdaptation(mesh, order=order+1, threshold=0.1,discontinuous_qn=True)
deformation = lsetmeshadap.CalcDeformation(levelset)
lsetp1 = lsetmeshadap.lset_p1

We will need several different element and facet markers for the setup of the method:

ci = CutInfo(mesh, lsetp1)
hasneg = ci.GetElementsOfType(HASNEG) # elements with negative level set (root elements)
pos = ci.GetElementsOfType(POS)       # outside elements
hasif = ci.GetElementsOfType(IF)      # cut elements
hasany = ci.GetElementsOfType(ANY)    # all elements (trivial array)

We define patches (for minimizing GP stabilization (if used) and postprocessings). Here, we mark all cut elements as "bad", i.e. they need to be supported.

EA = ElementAggregation(mesh)
EA.Update(ci.GetElementsOfType(NEG), hasif)
patch_number_field = GridFunction(L2(mesh))
patch_number_field.vec.FV()[:] = EA.element_to_patch
Draw(patch_number_field, mesh, "patch_number_field") #, deformation=CF((0,0,0.02*patch_number_field)))

BaseWebGuiScene

Solving the problem

We consider the exact solution \(p = \sin(x)\) and derive \(u = \nabla p\) and \(f = - \Delta p\). For simplicity, we assume here that we are given a proper extension of \(f\) from \(\Omega\) to the active mesh by simply evaluating the same closed form expression. This can be generalized to deal with functions only defined on \(\Omega\). This is discussed in a separate section below.

p_exact = sin(x)
coeff_f = - (p_exact.Diff(x).Diff(x) + p_exact.Diff(y).Diff(y)).Compile()
u_exact = CF((p_exact.Diff(x),p_exact.Diff(y)))

Finite element spaces

We use a Raviart-Thomas finite element space for the flux and a standard L2 space for the scalar variable, both are restricted to the active mesh. Note, that if we use ghost penalty stabilization we have to pay with additional couplings on cut elements, i.e. we need to set dgjumps=True in this case.

#FESpaces
Shsbase = HDiv(mesh, order=order, dirichlet=[], dgjumps=(gamma_stab > 0), RT=True)
Vhbase = L2(mesh, order=order, dirichlet=[])
Vh = Restrict(Vhbase, hasneg)
Shs = Restrict(Shsbase, hasneg)
Wh = Shs*Vh

GridFunction, trial- and test functions, and some standard coefs

gfw = GridFunction(Wh)
gfu, gfpT = gfw.components[0:2]
(uh,pT), (vh,qT) = Wh.TnT()

h = specialcf.mesh_size     # mesh size
n = Normalize(grad(lsetp1)) # normal to level set

Differential symbols for different integration domains:

# active mesh integrals:
dxbar = dx(definedonelements=ci.GetElementsOfType(HASNEG), deformation=deformation)
# uncut mesh integrals:
dxinner = dx(definedonelements=ci.GetElementsOfType(NEG), deformation=deformation)
# cut element integrals (full elements, no cut integral):
dxcut = dx(definedonelements=ci.GetElementsOfType(IF), deformation=deformation)
# integral on zero level set:
ds = dCut(lsetp1, IF, definedonelements=hasif, deformation=deformation)
# integral on Omega (physical domain, this is a cut integral):
dX = dCut(lsetp1, NEG, definedonelements=hasneg, deformation=deformation)

If ghost penalties are used, we further need a facet patch integral on some facet patches. Here, we only stabilize on facets that lie within a stabilizing patch, cf. aggregation tutorial for details.

if gamma_stab > 0:
    #Integration domain for ghost penalty
    dP = dFacetPatch(definedonelements=EA.patch_interior_facets, deformation=deformation)

Next, we create the BilinearForm and again distinguish the ghost penalty case from the unstabilized case, because the GP case require a larger sparsity pattern, however only on the stabilizing patches:

if gamma_stab > 0:
    a = RestrictedBilinearForm(Wh, element_restriction=hasneg,
                               facet_restriction=EA.patch_interior_facets, symmetric=True)
else:
    a = BilinearForm(Wh, symmetric=True,condense=False)
f = LinearForm(Wh)

Definition of variational formulation:

Now, we implement the terms in \((M1)\) and \((M2)\)

a += uh*vh * dX
a += (div(uh) * qT + div(vh) * pT) * dxbar
if gamma_stab > 0: # ghost penalty
    a += gamma_stab * (uh - uh.Other()) * (vh - vh.Other()) * dP
f += -coeff_f * qT * dxbar
f += p_exact * vh * n * ds

Assembly and solution

a.Assemble()
f.Assemble()
gfw.vec.data = a.mat.Inverse(Wh.FreeDofs()) * f.vec

Draw(BitArrayCF(hasneg)*gfw.components[0],mesh,"u",deformation=deformation, min=0, max=1.2) # flux variable on active mesh

BaseWebGuiScene

DrawDC(lsetp1,gfw.components[1],-1,mesh,"p",deformation=deformation, min=-1, max=1)        # p on physical domain

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

BaseWebGuiScene

We can now measure the error and considered different measures:

• p_l2error = \(\Vert p - p_h \Vert_{\Omega}\) (the error on the physical domain)

• p_inner_l2error = \(\Vert p - p_h \Vert_{\mathcal{T}_h^{\text{int}}}\) (the error only on uncut elements)

• u_l2error = \(\Vert u - u_h \Vert_{\Omega}\) (the error on the physical domain)

• u_l2error_bar = \(\Vert u - u_h \Vert_{\mathcal{T}_h}\) (the error on the whole active mesh)

p_l2error = sqrt(Integrate((gfpT - p_exact)**2*dX.order(2*order+3), mesh))
p_inner_l2error = sqrt(Integrate((gfpT - p_exact)**2*dxinner, mesh))
u_l2error = sqrt(Integrate((gfu - u_exact)**2*dX.order(2*order+3), mesh))
u_l2error_bar = sqrt(Integrate((gfu - u_exact)**2*dxbar, mesh))
print("p_l2error = ", p_l2error)
print("p_inner_l2error = ", p_inner_l2error)
print("u_l2error = ", u_l2error)
print("u_l2error_bar = ", u_l2error_bar)

p_l2error =  0.06454534811028939
p_inner_l2error =  9.366546623268328e-05
u_l2error =  0.00017930004860997905
u_l2error_bar =  0.03653877600813849

We observe several (not unexpected) effects:

1. The error for \(u\) is good on the physical domain, but much worse on the active mesh

2. The error for \(p\) is not even good on the physical domain, but only on the uncut elements

The first effect can be cured when applying a ghost penalty stabilization which ensures a smooth extension of the discrete solution onto the whole active mesh.

The second effect comes from the inconsistency that we introduced when changing the integration domain from \(\Omega\) to the active mesh. The latter problem can be cured by reconstructing a better, consistent approximation of \(p\) based on good interior values and a good approximation of \(u\) by postprocessing, which we discuss next.

Patchwise Post-processing

In the following post-processing we exploit that \(u\) is (optimally) accurate on \(\Omega\) (independent of GP choice) and that \(p\) is accurate in the interior. We then apply standard post-processing strategies on uncut elements and use an adapted patch-wise version for cut elements.

On each patch \(\omega\) with subtriangulation \(\mathcal{T}_\omega\), we want to find \(p_h^\ast \in \mathbb{P}^{k+1}(\mathcal{T}_\omega)\) so that \begin{align} (\nabla p_h^\ast, \nabla q_h^\ast)_{\Omega \cap \omega} + j_h^\omega(p_h^\ast,q_h^\ast) &= (u_h,\nabla q_h^\ast)_{\Omega \cap \omega} \qquad \forall q_h^\ast \in \mathbb{P}^{k+1}(\mathcal{T}_\omega) \setminus \mathbb{R} \tag{P1-local} \\ (p_h^\ast,1)_{\omega \cap \Omega^{\text{int}}} &= (\bar p_h,1)_{\omega \cap \Omega^{\text{int}}}. \tag{P2-local} \end{align} Here \(j_h^{\omega}(\cdot,\cdot)\) is a patch-local ghost penalty with an \(H^1\) (instead of \(L^2\)) scaling.

We can reformulate this as a global problem on the whole mesh, especially the integral constraint, with the help of a Lagrange multiplier \(\lambda \in \mathbb{P}^0(\mathcal{T}_h^{\text{int}})\): Find \(p_h^s \in \mathbb{P}^{k+1}(\mathcal{T}_h)\) and \(\lambda \in \mathbb{P}^0(\mathcal{T}_h^{\text{int}})\) so that \begin{align} (\nabla p_h^\ast, \nabla q_h^\ast)_{\Omega} + j_h^\omega(p_h^\ast,q_h^\ast) + (q_h^\ast,\lambda)_{\mathcal{T}_h^{\text{int}}} &= (u_h,\nabla q_h^\ast)_{\Omega} \qquad \forall q_h^\ast \in \mathbb{P}^{k+1}(\mathcal{T}_h) \tag{P1} \\ (p_h^\ast,\mu)_{\mathcal{T}_h^{\text{int}}} &= (p_h^\ast,\mu)_{\mathcal{T}_h^{\text{int}}} \qquad \forall \mu \in \mathbb{P}^0(\mathcal{T}_h^{\text{int}}). \tag{P2} \end{align}

Let’s start with the corresponding definition of the spaces (and corresponding trial/test functions and a solution GridFunction):

#FE Spaces for el-wise post-processing
Vhpbase = L2(mesh, order=order+1, dgjumps=True, dirichlet=[]) # for p
Lhbase = L2(mesh, order=0, dgjumps=True, dirichlet=[])        # for integral constraints
Vhp = Restrict(Vhpbase, hasneg)
Lh = Restrict(Lhbase, ci.GetElementsOfType(NEG)) # one constraint per patch / root element
Zh = Vhp * Lh
(ps,lam),(vs,mu) = Zh.TnT()

gfz = GridFunction(Zh)
gfps, gflam = gfz.components

Then, we setup the formulation \((P1)\)-\((P2)\):

#Integration domain for ghost penalty
dP = dFacetPatch(definedonelements=EA.patch_interior_facets, deformation=deformation)

lhs_integrals = [grad(ps) * grad(vs) * dX,
                 (lam * vs + mu * ps) * dxinner,
                 1/h**2 * (ps - ps.Other()) * (vs - vs.Other()) * dP]
rhs_integrals = [gfu * grad(vs) * dX,
                 gfpT * mu * dxinner]

Now, we could put this into a global Matrix and vector, as usual, or we exploit the a-priori knowledge that the corresponding linear system decouples into the patches. For this we can use the PatchwiseSolve:

help(PatchwiseSolve)

Help on built-in function PatchwiseSolve in module xfem.xfem:

PatchwiseSolve(...) method of builtins.PyCapsule instance
    PatchwiseSolve(elagg: xfem.xfem.ElementAggregation, fes: ngsolve.comp.FESpace, bf: ngsolve.comp.SumOfIntegrals, lf: ngsolve.comp.SumOfIntegrals, heapsize: int = 1000000) -> ngsolve.la.BaseVector


    Solve patch-wise problem based on the patches provided by the element aggregation input.

    Parameters

    elagg: ElementAggregatetion
      The instance defining the patches

    fes : FESpace
      The finite element space on which the local solve is performed.

    bf : SumOfIntegrals
      Integrators defining the matrix problem.

    lf : SumOfIntegrals
      Integrators defining the right-hand side.

gfz.vec.data = PatchwiseSolve(EA,Zh,sum(lhs_integrals),sum(rhs_integrals))

Let’s take a look at the postprocessed solution:

DrawDC(lsetp1,gfps,-1,mesh,"p",deformation=deformation, min=-1, max=1)        # p on physical domain

BaseWebGuiScene

\(p_h^\ast\) looks much better now. Let’s also take a look at the error:

ps_l2error = sqrt(Integrate((gfps - p_exact)**2*dX.order(2*order+3), mesh))
print("ps_l2error = ", ps_l2error)

ps_l2error =  0.0006351305619653827

And also the error is much improved by the post-processing. This concludes the main part of this tutorial.

What’s left?

There are several follow-ups. Below, in the appendix we discuss some variants:

1. How to deal with the r.h.s. if it is not provided on the active mesh, but only on \(\Omega\)

2. An alternative post-processing that work element-by-element in case that a ghost penalty stabilization has been used before

Several other tasks are left for you to play around with:

• What is the convergence order for the considered fields in the discretization?

• How can you implement Neumann boundary conditions (with or without polluting the mass balance)

• If no GP stabilization is applied the coupling structure is the same as for body-fitted mixed methods. Hence, you can also hybridize the system yielding a symmetric positive definite system for the hybrid variables after static condensation. Try it out!

Further explanations on the method, the above special cases and a demo file are provided through our preprint(arXiv:2306.12722). Have a look!

Appendix: extensions

Instead of assuming \(f\) to be given on the whole active mesh, we can also use the ghost penalty stabilization to construct an approximation \(f_h\) of \(f\) by a discrete extension. For \(\gamma_f > 0\), we find \(f_h \in \mathbb{P}^{k_f}(\mathcal{T}_h)\), \(k_f \in \mathbb{N}_0\), by solving \begin{equation} (f_h,q_h)_{\Omega} + \gamma_f j_h (f_h,q_h) = (f,q_h)_{\Omega} \qquad \forall q_h \in \mathbb{P}^{k_f}(\mathcal{T}_h). \end{equation}

This problem can be solved patch-by-patch if \(j_h(\cdot,\cdot)\) decouples across patches. The following code implements this approximation using patch-wise solves.

kf = order
gfh = GridFunction(L2(mesh,order=kf))
fh, wh = gfh.space.TnT()
lhs_integrals = fh * wh * dX + 0.01*(fh-fh.Other())*(wh-wh.Other()) * dFacetPatch(deformation=deformation)
rhs_integrals = coeff_f * wh * dX
gfh.vec.data = PatchwiseSolve(EA,gfh.space,lhs_integrals,rhs_integrals)

If \(u\) is well approximated not only on \(\Omega\), but on the whole active mesh, i.e. if GP stabilization has been used, a simpler element-by-element post-processing scheme is possible.

By making use of the relation \(u=\nabla p\), we can repair the inconsistencies of the approximation of \(p\) and achieve higher order convergence. For each element \(T \in \mathcal{T}_h\), we want to find \(p_h^\ast \in \mathcal{P}^{k+1}(T)\) such that \begin{alignat}{2} (\nabla p_h^\ast, \nabla q_h^\ast)_T &= (u_h, \nabla q_h^\ast)_T \qquad &&\forall q_h^\ast \in \mathcal{P}^{k+1}(T)\\ (p_h^\ast,1)_T &= (\bar p_h,1)_T \qquad &&\text{if } T \in \mathcal{T}_h^{\text{int}}\\ (p_h^\ast,1)_{T \cap \partial \Omega} &= (p_D,1)_{T \cap \partial \Omega} \qquad &&\text{if } T \in \mathcal{T}_h^{\text{cut}}\\ \end{alignat}

#FE Spaces for el-wise post-processing
Vhpbase = L2(mesh, order=order+1, dgjumps=False, dirichlet=[]) # order+1
Lhbase = L2(mesh, order=0, dgjumps=False, dirichlet=[])

Vhp = Restrict(Vhpbase, hasneg)
Lh = Restrict(Lhbase, hasneg)
Zh = Vhp * Lh

gflh = GridFunction(Lh)

gfz = GridFunction(Zh)
gfps, gflam = gfz.components

#Test- & Trialfunction
(ps,lam),(vs,mu) = Zh.TnT()

#Bilinear Form
p = RestrictedBilinearForm(Zh, symmetric=False)
p += grad(ps) * grad(vs) * dxbar
p += (lam * vs + mu * ps) * dxinner
p += (lam * vs + mu * ps) * ds

# R.h.s. term:
pf = LinearForm(Zh)
pf += gfu * grad(vs) * dxbar

pf += p_exact * mu * ds
pf += gfpT * mu * dxinner

#Assembly
p.Assemble()
pf.Assemble()

#Solving the system
gfz.vec.data = p.mat.Inverse(Zh.FreeDofs(),inverse="umfpack") * pf.vec

#For visualization:
gflh.Set(p_exact, definedonelements=ci.GetElementsOfType(HASNEG))
gflam.Set(gfps, definedonelements=ci.GetElementsOfType(HASNEG))

DrawDC(lsetp1,gfps,-1,mesh,"p",deformation=deformation, min=-1, max=1)        # p on physical domain

BaseWebGuiScene

ps_l2error = sqrt(Integrate((gfps - p_exact)**2*dX.order(2*order+3), mesh))
print("ps_l2error = ", ps_l2error)

ps_l2error =  3.8321764785411405e-06