This page was generated from unit-8.1-basics/basics.ipynb.

8.1 Fundamental concepts of ngsxfem

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

• level set geometries and its (P1) approximation (for higher order approximations see intlset.ipynb)

• Visualization of discontinuous (across level set zeros) functions

• handling of cut configurations inside the mesh (active elements / facets / dofs)

• Restriction of FESpaces, integrators and Matrices to active parts.

Geometry description with the level set function

Most features of the ngsxfem add-on to NGSolve are based on a level set description.

Let \(\phi\) be a continuous level set function.

We define the domains and the interface

\[\begin{split}\begin{array}{rccc} \text{domain:} & \Omega_{-} := \{ \phi < 0 \}, & \quad \Omega_{+} := \{ \phi > 0 \}, & \quad \Gamma := \{ \phi = 0 \}. \\ \texttt{ngsxfem} \text{ keyword:} & \texttt{NEG} & \texttt{POS} & \texttt{IF} \end{array}\end{split}\]

Let us import Netgen/NGSolve stuff first:

#import netgen.gui
# ngsolve stuff
from ngsolve import *
from ngsolve.webgui import *
# basic geometry features (for the background mesh)
from netgen.geom2d import SplineGeometry
# visualization stuff
from ngsolve.internal import *

We add the basic functionality of the ngsxfem package:

# basic xfem functionality
from xfem import *

importing ngsxfem-2.1.2303.dev0

The background mesh

We generate the background mesh of the domain

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

On the background mesh we define the level set function. In the visualization – using DrawDC – we distinguish the negative and the positive part of the level set values corresponding to the subdomains that are described.

levelset = (sqrt(sqrt(x**4+y**4)) - 1.0)
DrawDC(levelset, -1.0, 2.0, mesh, "x")
help(DrawDC)

Help on partial in module functools:

functools.partial(<function DrawDiscontinuous_we...f0c3520>, <function _DrawDocu at 0x7ff12f0c0ee0>)
    Visualization method for functions that are non-smooth across
    level set interfaces. Effectively calls netgen.webgui.Draw with
    a few manipulations.

    Parameters
    ----------
    levelset : CoefficientFunction
        (scalar) CoefficientFunction that describes the (implicit) geometry
    fneg : CoefficientFunction
        CoefficientFunction that is evaluated where the level set function is negative
    fpos : CoefficientFunction
        CoefficientFunction that is evaluated where the level set function is positive

    deformation : deformation (optional)
        (vectorial) CoefficientFunction that describes the deformation of the background mesh
    *remainder* : *
        all remainder arguments are passed to netgen.webgui.Draw

Geometry approximation through level set approximation

• Arbitrary level set functions can be difficult to handle (cut topology)

• To obtain something feasible, we interpolate into the space of piecewise linear functions

• The thusly obtained approximation is straight within each (simplex) element

• Simplifies numerical integration

• We will discuss enhancements later.

lsetp1 = GridFunction(H1(mesh,order=1))
InterpolateToP1(levelset,lsetp1)
DrawDC(lsetp1,-1,1,mesh,"lsetp1")

BaseWebGuiScene

• InterpolateToP1 takes vertex values (in contrast to Set)

Cut Information

In unfitted FEM it is important to know which elements are cut by the interface \(\{ \phi = 0 \}\). These informations are gathered in the cut information class:

ci = CutInfo(mesh, lsetp1)
help(CutInfo)

Help on class CutInfo in module xfem.xfem:

class CutInfo(pybind11_builtins.pybind11_object)
 |  A CutInfo stores and organizes cut informations in the mesh with respect to a level set function.
 |  Elements (BND / VOL) and facets can be either cut elements or in the positive (POS) or negative
 |  (NEG) part of the domain. A CutInfo provides information about the cut configuration in terms of
 |  BitArrays and Vectors of Ratios. (Internally also domain_types for different mesh nodes are stored.)
 |
 |  Method resolution order:
 |      CutInfo
 |      pybind11_builtins.pybind11_object
 |      builtins.object
 |
 |  Methods defined here:
 |
 |  GetCutRatios(...)
 |      GetCutRatios(self: xfem.xfem.CutInfo, VOL_or_BND: ngsolve.comp.VorB = <VorB.VOL: 0>) -> ngsolve.la.BaseVector
 |
 |
 |      Returns Vector of the ratios between the measure of the NEG domain on a (boundary) element and the
 |      full (boundary) element
 |
 |  GetElementsOfType(...)
 |      GetElementsOfType(self: xfem.xfem.CutInfo, domain_type: object = <DOMAIN_TYPE.IF: 2>, VOL_or_BND: ngsolve.comp.VorB = <VorB.VOL: 0>) -> pyngcore.pyngcore.BitArray
 |
 |
 |      Returns BitArray that is true for every element that has the
 |      corresponding combined domain type
 |      (NO/NEG/POS/UNCUT/IF/HASNEG/HASPOS/ANY)
 |
 |  GetElementsWithThresholdContribution(...)
 |      GetElementsWithThresholdContribution(self: xfem.xfem.CutInfo, domain_type: object = <DOMAIN_TYPE.NEG: 0>, threshold: float = 1.0, VOL_or_BND: ngsolve.comp.VorB = <VorB.VOL: 0>) -> pyngcore.pyngcore.BitArray
 |
 |
 |      Returns BitArray marking the elements where the cut ratio is greater or equal to the given
 |      threshold.
 |
 |      Parameters
 |
 |      domain_type : ENUM
 |          Check POS or NEG elements.
 |
 |      threshold : float
 |          Mark elements with cut ratio (volume of domain_type / volume background mesh) greater or equal to threshold.
 |
 |      VOL_or_BND : ngsolve.comp.VorB
 |          input VOL, BND, ..
 |
 |  GetFacetsOfType(...)
 |      GetFacetsOfType(self: xfem.xfem.CutInfo, domain_type: object = <DOMAIN_TYPE.IF: 2>) -> pyngcore.pyngcore.BitArray
 |
 |
 |      Returns BitArray that is true for every facet that has the
 |      corresponding combined domain type
 |      (NO/NEG/POS/UNCUT/IF/HASNEG/HASPOS/ANY)
 |
 |  Mesh(...)
 |      Mesh(self: xfem.xfem.CutInfo) -> ngsolve.comp.Mesh
 |
 |
 |      Returns mesh of CutInfo
 |
 |  Update(...)
 |      Update(self: xfem.xfem.CutInfo, levelset: ngsolve.fem.CoefficientFunction, subdivlvl: int = 0, time_order: int = -1, heapsize: int = 1000000) -> None
 |
 |
 |      Updates a CutInfo based on a level set function.
 |
 |      Parameters
 |
 |      levelset : ngsolve.CoefficientFunction
 |        level set function w.r.t. which the CutInfo is generated
 |
 |      subdivlvl : int
 |        subdivision for numerical integration
 |
 |      time_order : int
 |        order in time that is used in the integration in time to check for cuts and the ratios. This is
 |        only relevant for space-time discretizations.
 |
 |  __init__(...)
 |      __init__(self: xfem.xfem.CutInfo, mesh: ngsolve.comp.Mesh, levelset: object = None, subdivlvl: int = 0, time_order: int = -1, heapsize: int = 1000000) -> None
 |
 |
 |      Creates a CutInfo based on a level set function and a mesh.
 |
 |      Parameters
 |
 |      mesh : Mesh
 |
 |      levelset : ngsolve.CoefficientFunction / None
 |        level set funciton w.r.t. which the CutInfo is created
 |
 |      time_order : int
 |        order in time that is used in the integration in time to check for cuts and the ratios. This is
 |        only relevant for space-time discretizations.
 |
 |  ----------------------------------------------------------------------
 |  Static methods inherited from pybind11_builtins.pybind11_object:
 |
 |  __new__(*args, **kwargs) from pybind11_builtins.pybind11_type
 |      Create and return a new object.  See help(type) for accurate signature.

Marking elements

DOMAIN_TYPE

We can ask for BitArrays corresponding to the differently marked elements

• NEG : True if \(\phi < 0\) everywhere on this element, else False

• POS : True if \(\phi > 0\) everywhere on this element, else False

• IF : True if \(\phi = 0\) somewhere on this element, else False

dts = [IF,NEG,POS]
def DrawMarkedElements():
    for dt in dts:
        print("Drawing", dt);
        Draw(BitArrayCF(ci.GetElementsOfType(dt)),mesh,"elements_"+str(dt))
        yield
dme = DrawMarkedElements()

next(dme)

Drawing DOMAIN_TYPE.IF

next(dme)

Drawing DOMAIN_TYPE.NEG

next(dme)

Drawing DOMAIN_TYPE.POS

COMBINED_DOMAIN_TYPE

There are also some predefined combinations of the previous markers:

IF	POS	NEG	combined flag	description
0	0	0	NO	False on all elements
0	0	1	NEG	True if \(\phi < 0\) everywhere on this element, else False
0	1	0	POS	True if \(\phi > 0\) everywhere on this element, else False
0	1	1	UNCUT	True if \(\phi \neq 0\) everywhere on this element, else False
1	0	0	IF	True if \(\phi = 0\) somewhere on this element, else False
1	0	1	HASNEG	True if \(\phi \leq 0\) somewhere on this element, else False
1	1	0	HASPOS	True if \(\phi \geq 0\) somewhere on this element, else False
1	1	1	ANY	True on all elements

cdts = [NO,CDOM_NEG,CDOM_POS,UNCUT,CDOM_IF,HASNEG,HASPOS,ANY]
for dt in cdts: print(int(dt),": {0:03b} ,".format(int(dt)),str(dt) )

0 : 000 , COMBINED_DOMAIN_TYPE.NO
1 : 001 , COMBINED_DOMAIN_TYPE.CDOM_NEG
2 : 010 , COMBINED_DOMAIN_TYPE.CDOM_POS
3 : 011 , COMBINED_DOMAIN_TYPE.UNCUT
4 : 100 , COMBINED_DOMAIN_TYPE.CDOM_IF
5 : 101 , COMBINED_DOMAIN_TYPE.HASNEG
6 : 110 , COMBINED_DOMAIN_TYPE.HASPOS
7 : 111 , COMBINED_DOMAIN_TYPE.ANY

for dt in cdts: print(dt);Draw(BitArrayCF(ci.GetElementsOfType(dt)),mesh,"elements_"+str(dt),min=0,max=1)

COMBINED_DOMAIN_TYPE.NO

COMBINED_DOMAIN_TYPE.CDOM_NEG

COMBINED_DOMAIN_TYPE.CDOM_POS

COMBINED_DOMAIN_TYPE.UNCUT

COMBINED_DOMAIN_TYPE.CDOM_IF

COMBINED_DOMAIN_TYPE.HASNEG

COMBINED_DOMAIN_TYPE.HASPOS

COMBINED_DOMAIN_TYPE.ANY

Updating the CutInformation

• level set function can be changed with Update:

levelset = sqrt((x-0.5)**2+y**2)-1
scene = Draw(levelset,mesh,"levelset") # overwrite visualization
InterpolateToP1(levelset,lsetp1)
ci.Update(lsetp1);scene.Redraw()

Marking with facets

• GetFacetsWithNeighborTypes

• GetElementsWithNeighborFacets

• Sometimes want to mark some facets depending on their neighbor states (NEG/POS/IF)

• Needed for some stabilizations in unfitted FEM (CutFEM/XFEM/..)

• Facet-marking based on two element-based BitArrays.

• There are two options to connect two BitArrays A and B:

• use_and=True (default): the result is True for a facet if (A(left) and B(right)) or (B(left) and A(right))

• use_and=False (or): the result is True for a facet if (A(left) or B(right)) or (B(left) or A(right))

A(left)	B(right)	partial result (use_and=True)	partial result (use_and=False)
False	False	False	False
False	True	False	True
True	False	False	True
True	True	True	True

• with bnd_val_a and bnd_val_b (default: True) we can prescribe how the BitArray should be interpreted if there is only one neighbor (boundary).

Example 1:

• All facets that are on the outer boundary of the NEG-domain:

#help(GetFacetsWithNeighborTypes)
ba_facets = GetFacetsWithNeighborTypes(mesh,
                                       a=ci.GetElementsOfType(NEG),
                                       b=ci.GetElementsOfType(IF))

• We can also go from facets to elements

#help(GetElementsWithNeighborFacets)
ba_surround_facets = GetElementsWithNeighborFacets(mesh,ba_facets)

def SelectFacetsToDraw(dt1,dt2):
    ba_facets = GetFacetsWithNeighborTypes(mesh,a=ci.GetElementsOfType(dt1), b=ci.GetElementsOfType(dt2))
    ba_surround_facets = GetElementsWithNeighborFacets(mesh,ba_facets)
    Draw(BitArrayCF(ba_surround_facets), mesh, "surrounding_facets")

SelectFacetsToDraw(CDOM_IF,CDOM_NEG)

Example 2:

We can successively apply this to extend from some part of the domain by direct neighbors, etc ….

(execute the next cell several times)

ba_facets = GetFacetsWithNeighborTypes(mesh,a=ba_surround_facets,b=ba_surround_facets,
                                           bnd_val_a=False,bnd_val_b=False,use_and=False)
ba_surround_facets = GetElementsWithNeighborFacets(mesh,ba_facets)
Draw(BitArrayCF(ba_surround_facets), mesh, "ba_surround")

BaseWebGuiScene

We now learned how to mark elements and facets. Next, we want to apply this to adapt finite element formulations.

Restriction to active mesh/dofs

• Unfitted FEM: Solve a PDE problem restricted to a subdomain, e.g.

  • a fictitious domain method on \(\Omega_{-} = \{ \phi < 0 \}\) or

  • a trace finite element on \(\Gamma = \{ \phi = 0\}\)

• Standard FESpace but only a selection of dofs (active_dofs)

• We have helper functions for that as well

Marking dofs (GetDofsOfElements)

With GetDofsOfElements we can select a subset of dofs corresponding to a set of elements. With Compress we can reduce a finite element space to the corresponding dofs:

Vhbase = H1(mesh, order=1, dirichlet=[])
VhC = Compress(Vhbase,GetDofsOfElements(Vhbase,ci.GetElementsOfType(HASNEG)))

Alternatively, you can directly Restrict a finite element space to a set of elements:

VhR = Restrict(Vhbase,ci.GetElementsOfType(HASNEG))

Now, we can compare:

print(Vhbase.ndof, VhC.ndof, VhR.ndof)

56 30 30

Compress and Restrict yield the same amount of dofs. However, there is a subtle difference: VhR is not defined away from the marked elements whereas VhC can be seen as the background FESpace with the inactive dofs set to zero.

Next:

• Compress and Restrict allow to restrict vectors and matrices to "active" subset of the background FE space

• What to do if you want to set up the finite element formulation only on a selection of elements / facets / unknowns?

Marking integrators (definedonelements)

We can define integrators on a set of elements only (this is a standard NGSolve feature):

integrator = VhR.TrialFunction()*VhR.TestFunction() * dx(definedonelements=ci.GetElementsOfType(HASNEG))

print(ci.GetElementsOfType(HASNEG))

0: 00000000000111111010000000000000000000000101111111
50: 111011100000001110101111111111111111

Note that the BitArray taken in the definedonelements argument is taken by reference, i.e. later changes to the CutInfo ci will automatically change the set of elements the integrator will act on.

Restricted Bilinearforms

Furthermore, to avoid that matrices are set up for the full standard FESpace we can restrict matrix couplings to elements (and facets for dg terms). This may be attractive if you do not work with Restricted or Compressed finite element spaces (e.g. because of changing domains) but want the linear systems to be of reduced size. Further, in unfitted FEM often stabilizations of "Ghost Penalty" type are applied which may lead to dgjumps, i.e. couplings of dofs over facets. To avoid that these couplings are introduced on the whole domain, you may want to restrict the facet couplings explicitly.

#help(RestrictedBilinearForm)
a_full = BilinearForm(Vhbase, check_unused=False)
a = RestrictedBilinearForm(Vhbase,element_restriction=ci.GetElementsOfType(IF),
                           facet_restriction=None, check_unused=False)
a_full.Assemble()
a.Assemble()

<xfem.xfem.RestrictedBilinearFormDouble at 0x7ff0b624b1a0>

• This does not change the dimension:

print(a.mat.height,a.mat.width)
print(a_full.mat.height,a_full.mat.width)

56 56
56 56

• But the number on non-zero entries:

print(len(a_full.mat.AsVector()))
print(len(a.mat.AsVector()))

338
197