This page was generated from unit-8.7-mlset_pde/mlset_pde.ipynb.

# 8.7 Unfitted FEM for domains described by multiple level sets¶

In order to solve a partial differential equation on a stationary domain defined by multiple level set functions, we proceed similarly to the demo cutfem.ipynb, with a few key differences in order to deal with corners.

For simplicity we consider a Poisson problem on the unit square in order to show the differences between the single and multiple level set setting.

The libraries are imported as usual:

```
[1]:
```

```
# Basic NGSolve things
from netgen.geom2d import SplineGeometry
from ngsolve import *
from ngsolve.solvers import PreconditionedRichardson as PreRic
# ngsxfem and the mlset convenience layer
from xfem import *
from xfem.mlset import *
# Visualisation
from ngsolve.webgui import *
DrawDC = MakeDiscontinuousDraw(Draw)
```

```
importing ngsxfem-2.1.2303.dev0
```

We want to solve the problem

where the domain \(\Omega\subset\widetilde\Omega\) is a subset of the meshed background domain \(\widetilde\Omega\).

```
[2]:
```

```
geo = SplineGeometry()
geo.AddRectangle((-0.2, -0.2), (1.2, 1.2), bcs=("bottom", "right", "top", "left"))
mesh = Mesh(geo.GenerateMesh(maxh=0.1))
```

We define the domain as the set \(\Omega := \bigcap_i \{\phi_i < 0\}\) with \(\{\phi_i\} = \{ -y, x - 1, y - 1, -x \}\) and use a `DomainTypeArray`

as a container for this set and to generate the boundary.

```
[3]:
```

```
level_sets = [-y, x - 1, y - 1, -x]
level_sets_p1 = tuple(GridFunction(H1(mesh, order=1)) for i in range(len(level_sets)))
for i, lsetp1 in enumerate(level_sets_p1):
InterpolateToP1(level_sets[i], lsetp1)
DrawDC(lsetp1, -3.5, 2.5, mesh, "lsetp1_{}".format(i))
square = DomainTypeArray((NEG, NEG, NEG, NEG))
boundary = square.Boundary()
```

Finite Element spaces are defined as usual:

```
[4]:
```

```
k = 3
V = H1(mesh, order=k, dgjumps=True)
gfu = GridFunction(V)
```

## Multiple level set cut-information¶

As before, we use standard FE spaces restricted to the active part of the mesh. In order to handle multiple level sets, we need the `MultiLevelsetCutInfo`

class:

```
[5]:
```

```
mlci = MultiLevelsetCutInfo(mesh, level_sets_p1)
```

### Marking the correct elements¶

Marking the correct active elements needs to be handled differently in the multipe level set setting, since elements of the type `(IF,IF)`

are sometimes not part of the active domain, as can be seen in the following image:

In order to select the correct active elements, the `MultiLevelsetCutInfo`

class has a member function `GetElementsWithContribution()`

. This takes in a `tuple(ENUM)`

, a list of such tuples or a `DomainTypeArray`

```
[6]:
```

```
els_hasneg = mlci.GetElementsWithContribution(square)
els_if = mlci.GetElementsWithContribution(boundary)
Draw(BitArrayCF(els_hasneg), mesh, "els_hasneg")
Draw(BitArrayCF(els_if), mesh, "els_if")
```

```
[6]:
```

```
BaseWebGuiScene
```

To deal with arbitraty mesh-interface cuts we use ghost-penalty stabilisation. The elements needed for this can be obtained in the same manner as before and similarly, the free degrees of freedom can also be determined as before:

```
[7]:
```

```
facets_gp = GetFacetsWithNeighborTypes(mesh, a=els_hasneg, b=els_if,
use_and=True)
freedofs = GetDofsOfElements(V, els_hasneg) & V.FreeDofs()
```

We implement the homogeneous Dirichlet boundary conditions using Nitsche’s trick. As each boundaty section has a different unit-normal vector, we have to decompose the Nitsche term into the different segments. Each will be added as a separate weak form, so we also require element markers for each side seperately.

Since the easient and most flexible nomenclature for each segment is the tuple of `DOMAIN_TYPE`

s, we store these element markers in a dictionary using the tuple as the key:

```
[8]:
```

```
els_if_single = {}
for i, dtt in enumerate(boundary):
els_if_single[dtt] = mlci.GetElementsWithContribution(dtt)
Draw(BitArrayCF(els_if_single[dtt]), mesh, "els_if_singe"+str(i))
```

## Defining the weak form¶

As in any other unfitted problem, we define some parameters and get the trial and test functions and mesh-size coefficient functions as usual.

```
[9]:
```

```
gamma_n = 10
gamma_s = 0.5
u, v = V.TnT()
h = specialcf.mesh_size
```

For the normal vectors, we use the `DomainTypeArray.GetOuterNormals()`

functionality, see mlset_basic.ipynb

```
[10]:
```

```
normals = square.GetOuterNormals(level_sets_p1)
```

We then define the cut `SymbolicDifferentialSymbols`

for the volume, boundary and ghost-penalty facet-patch integrators, see also mlset_basic.ipynb. Note that we pass the `BitArrays`

on which the integrators are defined to the `SymbolicDiffeentialSymbol`

```
[11]:
```

```
dx = dCut(level_sets_p1, square, definedonelements=els_hasneg)
ds = {dtt: dCut(level_sets_p1, dtt, definedonelements=els_if_single[dtt])
for dtt in boundary}
dw = dFacetPatch(definedonelements=facets_gp)
```

### Integrators¶

As we have an unfitted problem, we use a `RestrictedBilinearForm`

```
[12]:
```

```
a = RestrictedBilinearForm(V, element_restriction=els_hasneg, facet_restriction=facets_gp, check_unused=False)
```

The diffusion and ghost_penalty forms can then be added normally

```
[13]:
```

```
a += InnerProduct(Grad(u), Grad(v)) * dx
a += gamma_s / (h**2) * (u - u.Other()) * (v - v.Other()) * dw
```

For the Nitsche terms, we loop over the normals dictionary. This gives us acess to both the boundary domain-tuples and the normal vectors:

```
[14]:
```

```
for bnd, n in normals.items():
a += -InnerProduct(Grad(u) * n, v) * ds[bnd]
a += -InnerProduct(Grad(v) * n, u) * ds[bnd]
a += (gamma_n * k * k / h) * InnerProduct(u, v) * ds[bnd]
```

We want the exact solution to be \(u = 16 x (1 - x) y (1 - y)\). To add the appropriate forcing term to the right-hand side, we do not have anything special

```
[15]:
```

```
f = LinearForm(V)
f += 32 * (y * (1 - y) + x * (1 - x)) * v * dx
```

## Solving the system¶

Now the system has been set up, we can assemble the systrem and sove in the usual manner

```
[16]:
```

```
f.Assemble()
a.Assemble()
inv = a.mat.Inverse(freedofs=freedofs, inverse="sparsecholesky")
gfu.vec.data = PreRic(a=a, rhs=f.vec, pre=inv, freedofs=freedofs)
```

```
it = 0 ||res||_2 = 1.0491905634098577
it = 1 ||res||_2 = 9.62273091012537e-15
```

We can now visualise the solution:

```
[17]:
```

```
u_ex = 16 * x * (1 - x) * y * (1 - y)
grad_u_ex = CoefficientFunction((16 * (1 - 2 * x) * y * (1 - y),
16 * x * (1 - x) * (1 - 2 * y)))
DrawDC(square.Indicator(level_sets_p1), 0, gfu, mesh, "solution")
```

```
[17]:
```

```
BaseWebGuiScene
```

And compute the error. To integrate with quadrature rules of the correct order we generate new cut `SumbolicDifferentialSymbol`

s using the `.order(k)`

method with returns a new `DifferentialSymbol`

which tells the `Integrate`

to integrate with order `k`

```
[18]:
```

```
err_l2 = sqrt(Integrate((gfu - u_ex)**2 * dx.order(2 * k), mesh=mesh))
err_h1 = sqrt(Integrate((Grad(gfu) - grad_u_ex)**2 * dx.order(2 * (k - 1)), mesh=mesh))
print("L2 error = {:1.3e}".format(err_l2))
print("H1 error = {:1.3e}".format(err_h1))
```

```
L2 error = 3.497e-06
H1 error = 3.809e-04
```