# 8.2 Integration on level set domains¶

This unit gives a basic introduction to a few fundamental concepts in ngsxfem. Especially, we treat:

• level set geometries and its (P1) approximation

• Cut differential symbols to conveniently define integrals on cut domain

• for higher order level set geometry approximations

• simple example: calculate area of the unit circle using level set representations

• background domain: $$[-1.5, 1.5]^2 \subset \mathbb{R}^2$$.

Let $$\phi$$ be a continuous level set function. We recall:

$\Omega_{-} := \{ \phi < 0 \}, \quad \Omega_{+} := \{ \phi > 0 \}, \quad \Gamma := \{ \phi = 0 \}.$
[1]:

# basic geometry features (for the background mesh)
from netgen.geom2d import SplineGeometry
# ngsolve
from ngsolve import *
# basic xfem functionality
from xfem import *
# for isoparametric mapping
from xfem.lsetcurv import *
# visualisation
from ngsolve.webgui import *
# the contant pi
from math import pi

importing ngsxfem-2.1.2303.dev0


We generate the background mesh of the domain and use a simplicial triangulation

[2]:

square = SplineGeometry()
Draw(mesh)

[2]:

BaseWebGuiScene


On the background mesh we define the level set function. In this example we choose a levelset function which measures the distance to the origin and substracts a constant $$r$$, resulting in an circle of radius $$r$$ as $$\Omega_{-}$$:

[3]:

r = 1
levelset = sqrt(x**2+y**2) - r
DrawDC(levelset, -3.5, 2.5, mesh, 'levelset')

[3]:

BaseWebGuiScene

• Area of the unit circle:

$A = \pi = \int_{\Omega_-(\phi)} 1 \, \mathrm{d}x$

## Piecewise linear geometry approximation¶

• The ngsxfem integration on implicitly defined geometries is only able to handle functions $$\phi$$ which are in $$\mathcal{P}^1(T)$$ for all $$T$$ of the triangulation.

• Therefore we replace the levelset function $$\phi$$ with an approximation $$\phi^\ast \in \mathcal{P}^1(T), T \subset \Omega$$.

[4]:

V = H1(mesh,order=1,autoupdate=True)
lset_approx = GridFunction(V,autoupdate=True)
InterpolateToP1(levelset, lset_approx)
DrawDC(lset_approx, -3.5, 2.5, mesh, 'levelset_p1')

[4]:

BaseWebGuiScene


Notice that this replacement $$\phi \mapsto \phi^\ast$$ will introduce a geometry error of order $$h^2$$. We will investigate this numerically later on in this notebook.

As a next step we need to define the integrand:

[5]:

f = CoefficientFunction(1)


As in NGSolve, we integrate by using the Integrate function. For this to be aware of the fact that we are integration over a level set domain we use the dCut differential symbol, to which we need to pass:

• levelset: level set function which describes the geometry (best choice: P1-GridFunction)

• domain_type: decision on which part of the geometry the integration should take place (NEG/POS/IF)

• order: The order of integration rule to be used

[6]:

order = 2
integral = Integrate(f * dCut(levelset=lset_approx, domain_type=NEG, order=order), mesh=mesh)
error = abs(integral - pi)
print("Result of the integration: ", integral)
print("Error of the integration: ", error)

Result of the integration:  2.8532268778347776
Error of the integration:  0.28836577575501554


### Convergence study¶

To get more accurate results

• save the previous error and execute a refinement

• repeat the steps from above

• Repeatedly execute the lower block to add more refinement steps.

[7]:

errors = []

def AppendResultAndPostProcess(integral):
error = abs(integral - pi)
print("Result of the integration: ", integral)
print("Error of the integration: ", error)

errors.append(error)
eoc = [log(errors[i+1]/errors[i])/log(0.5) for i in range(len(errors)-1)]

print("Collected errors:", errors)
print("experimental order of convergence (L2):", eoc)

AppendResultAndPostProcess(integral)

Result of the integration:  2.8532268778347776
Error of the integration:  0.28836577575501554
Collected errors: [0.28836577575501554]
experimental order of convergence (L2): []

[8]:

scene = Draw(lset_approx, mesh, "lsetp1", min=0, max=0.0)
def RefineAndIntegrate():
# refine cut elements only:
RefineAtLevelSet(gf=lset_approx)
mesh.Refine()

InterpolateToP1(levelset,lset_approx)
scene.Redraw()
integral = Integrate(f * dCut(levelset=lset_approx, domain_type=NEG, order=order), mesh=mesh)
AppendResultAndPostProcess(integral)

[9]:

RefineAndIntegrate()

Result of the integration:  3.033381503911132
Error of the integration:  0.1082111496786613
Collected errors: [0.28836577575501554, 0.1082111496786613]
experimental order of convergence (L2): [1.4140507939031013]


We observe only a second order convergence which results from the piecewise linear geometry approximation.

## Higher order geometry approximation with an isoparametric mapping¶

In order to get higher order convergence we can use the isoparametric mapping functionality of xfem.

We apply a mesh transformation technique in the spirit of isoparametric finite elements:

Video of the mesh transformation

[10]:

%%html
<iframe width=100% height=600px src="https://www.youtube.com/embed/Mst_LvfgPCg?rel=0" frameborder="0" allow="autoplay; encrypted-media" allowfullscreen></iframe>


To compute the corresponding mapping we use the LevelSetMeshAdaptation class

[11]:

mesh = Mesh(square.GenerateMesh(maxh=0.4, quad_dominated=False))
Draw(deformation,mesh,"deformation")

[11]:

BaseWebGuiScene


We can observe the geometrical improvement in the following sequence:

[12]:

scene1 = Draw(deformation, mesh, "deformation", autoscale=False, min=0.0, max=0.028)

[13]:

scene2 = DrawDC(lsetmeshadap.lset_p1, -3.5, 2.5, mesh, "lsetp1_ho", deformation=deformation)

[14]:

N=100

deformation.vec.data = 1.0/N*deformation.vec
for i in range (1, N+1):
deformation.vec.data = (i+1)/i * deformation.vec
scene1.Redraw()
scene2.Redraw()


### Convergence study with isoparametric mapping applied¶

We now apply the mesh deformation (2nd order mapping) to improve the accuracy of the integration problem.

error tables:

[15]:

order = 2

levelset = sqrt(x*x+y*y)-1
referencevals = { POS : 9-pi, NEG : pi, IF : 2*pi }

errors_uncurved = dict()
errors_curved = dict()
eoc_uncurved = dict()
eoc_curved = dict()

for key in [NEG,POS,IF]:
errors_curved[key] = []
errors_uncurved[key] = []
eoc_curved[key] = []
eoc_uncurved[key] = []


refinements:

[16]:

f = CoefficientFunction (1.0)
refinements = 5
for reflevel in range(refinements):
if(reflevel > 0):
mesh.Refine()

for key in [NEG,POS,IF]:
# Applying the mesh deformation

integrals_uncurved = Integrate(f * dCut(levelset=lsetp1, domain_type=key, order=order), mesh=mesh)
# dCut deals with the mesh deformation internally
integrals_curved = Integrate(f * dCut(levelset=lsetp1, domain_type=key, order=order, deformation=deformation), mesh=mesh)

errors_curved[key].append(abs(integrals_curved - referencevals[key]))
errors_uncurved[key].append(abs(integrals_uncurved - referencevals[key]))
# refine cut elements:


Convergence rates:

[17]:

for key in [NEG,POS,IF]:
eoc_curved[key] = [log(a/b)/log(2) for (a,b) in zip (errors_curved[key][0:-1],errors_curved[key][1:]) ]
eoc_uncurved[key] = [log(a/b)/log(2) for (a,b) in zip (errors_uncurved[key][0:-1],errors_uncurved[key][1:]) ]

print("errors (  curved):  \n{}\n".format(  errors_curved))
print("   eoc (  curved):  \n{}\n".format(     eoc_curved))

errors (  curved):
{<DOMAIN_TYPE.NEG: 0>: [0.001710337974651921, 8.945537152715843e-05, 8.112397147819195e-06, 1.932554796191255e-06, 1.0494129565330468e-07], <DOMAIN_TYPE.POS: 1>: [0.11418698405525696, 8.945537152715843e-05, 8.112397149595552e-06, 1.9325547979676116e-06, 1.049413373976904e-07], <DOMAIN_TYPE.IF: 2>: [0.0009223373587694539, 0.0003690527231654528, 2.455580523097467e-05, 7.995576201125232e-07, 1.6669838220906286e-07]}

eoc (  curved):
{<DOMAIN_TYPE.NEG: 0>: [4.256969515213352, 3.4629679271130023, 2.0696189608663467, 4.202854908813071], <DOMAIN_TYPE.POS: 1>: [10.317942475320116, 3.4629679267970985, 2.069618959856161, 4.202854336252534], <DOMAIN_TYPE.IF: 2>: [1.321467599035857, 3.9096908984681504, 4.940718316506891, 2.261961902073998]}



Note: The integration on level set domains (and the mapping) can be extended to bilinear and linear form integrators with the same syntax.