This page was generated from unit-8.2-intlset/intlset.ipynb.

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

# 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

square = SplineGeometry()
square.AddRectangle([-1.5,-1.5], [1.5,1.5], bc=1)
mesh = Mesh(square.GenerateMesh(maxh=0.8, quad_dominated=False))
Draw(mesh)

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_{-}\):

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

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\).

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')

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:

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

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.

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): []

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)

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

%%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

mesh = Mesh(square.GenerateMesh(maxh=0.4, quad_dominated=False))
lsetmeshadap = LevelSetMeshAdaptation(mesh, order=2, threshold=1000, discontinuous_qn=True)
deformation = lsetmeshadap.CalcDeformation(levelset)
Draw(deformation,mesh,"deformation")

BaseWebGuiScene

We can observe the geometrical improvement in the following sequence:

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

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

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:

order = 2
mesh = Mesh (square.GenerateMesh(maxh=0.7, quad_dominated=False))

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

lsetmeshadap = LevelSetMeshAdaptation(mesh, order=order, threshold=0.2, discontinuous_qn=True)
lsetp1 = lsetmeshadap.lset_p1
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:

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
        deformation = lsetmeshadap.CalcDeformation(levelset)

        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:
    RefineAtLevelSet(gf=lsetmeshadap.lset_p1)

Convergence rates:

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.112397148263284e-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.4629679270340263, 2.0696189609453226, 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.