This page was generated from unit8.1basics/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
FESpace
s,integrators
andMatrices
to active parts.
Geometry description with the level set function¶
Most features of the ngsxfem
addon to NGSolve are based on a level set description.
Let \(\phi\) be a continuous level set function.
We define the domains and the interface
Let us import Netgen/NGSolve
stuff first:
[1]:
#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:
[2]:
# basic xfem functionality
from xfem import *
importing ngsxfem2.1.2303
The background mesh¶
We generate the background mesh of the domain
[3]:
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.
[4]:
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...5a11870>, <function _DrawDocu at 0x7faf85bc7490>)
Visualization method for functions that are nonsmooth 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.
[5]:
lsetp1 = GridFunction(H1(mesh,order=1))
InterpolateToP1(levelset,lsetp1)
DrawDC(lsetp1,1,1,mesh,"lsetp1")
[5]:
BaseWebGuiScene
InterpolateToP1
takes vertex values (in contrast toSet
)
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:
[6]:
ci = CutInfo(mesh, lsetp1)
help(CutInfo)
Help on class CutInfo in module xfem.ngsxfem_py:
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.ngsxfem_py.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.ngsxfem_py.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.ngsxfem_py.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.ngsxfem_py.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.ngsxfem_py.CutInfo) > ngsolve.comp.Mesh


 Returns mesh of CutInfo

 Update(...)
 Update(self: xfem.ngsxfem_py.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 spacetime discretizations.

 __init__(...)
 __init__(self: xfem.ngsxfem_py.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 spacetime 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 FalsePOS
: True if \(\phi > 0\) everywhere on this element, else FalseIF
: True if \(\phi = 0\) somewhere on this element, else False
[7]:
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()
[8]:
next(dme)
Drawing DOMAIN_TYPE.IF
[9]:
next(dme)
Drawing DOMAIN_TYPE.NEG
[10]:
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 
[11]:
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
[12]:
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
:
[13]:
levelset = sqrt((x0.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/..)
Facetmarking based on two elementbased 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 ( 
partial result ( 

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 NEGdomain:
[14]:
#help(GetFacetsWithNeighborTypes)
ba_facets = GetFacetsWithNeighborTypes(mesh,
a=ci.GetElementsOfType(NEG),
b=ci.GetElementsOfType(IF))
We can also go from facets to elements
[15]:
#help(GetElementsWithNeighborFacets)
ba_surround_facets = GetElementsWithNeighborFacets(mesh,ba_facets)
[16]:
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")
[17]:
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)
[18]:
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")
[18]:
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/dof
s¶
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 ofdof
s (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:
[19]:
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:
[20]:
VhR = Restrict(Vhbase,ci.GetElementsOfType(HASNEG))
Now, we can compare:
[21]:
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
andRestrict
allow to restrict vectors and matrices to "active" subset of the background FE spaceWhat 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):
[22]:
integrator = VhR.TrialFunction()*VhR.TestFunction() * dx(definedonelements=ci.GetElementsOfType(HASNEG))
[23]:
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 Restrict
ed or Compress
ed 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.
[24]:
#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()
[24]:
<xfem.ngsxfem_py.RestrictedBilinearFormDouble at 0x7faf85790900>
This does not change the dimension:
[25]:
print(a.mat.height,a.mat.width)
print(a_full.mat.height,a_full.mat.width)
56 56
56 56
But the number on nonzero entries:
[26]:
print(len(a_full.mat.AsVector()))
print(len(a.mat.AsVector()))
338
197