This page was generated from unit8.8aggregation/aggregation.ipynb.
8.8 Cell and basis aggregation in ngsxfem
¶
This unit gives a basic introduction to basis aggregation in ngsxfem
. Especially, we treat:
definition of aggregation patches
extension of basis functions on aggregation patches
[1]:
from ngsolve import *
from xfem import *
from ngsolve.webgui import *
importing ngsxfem2.1.dev
Introduction: Cut cells' best friends: their neighbors¶
Let’s start with a simple example (structured mesh, planar domain boundary):
[2]:
from ngsolve.meshes import MakeStructured2DMesh
mesh = MakeStructured2DMesh(quads=False, nx=3, ny=2)
levelset = x  0.8
gflset = GridFunction(H1(mesh))
InterpolateToP1(levelset, gflset)
DrawDC(gflset, 1.0, 2.0, mesh, "x")
[2]:
BaseWebGuiScene
In this simple configuration we have the right block of elements cut by the zero level set.
The cut configuration on cut elements can become arbitrarily bad so that basis functions on the cut elements may have degenerated support
\(\leadsto\) can lead to stability or (at least) conditioning issues.
To stabilize those cases we often seek for help in the neighborhood. This is typically done either
through ghost penalty stabilization where essentially functions are glued together weakly in a neighborhood or
by gluing basis functions on cut elements to basis functions on uncut elements in the virtue of an extrapolation
Patches¶
To define proper neighborhood regions where support is given from "good" to "bad" elements, we will define patches. First, we will need to define: * which elements are "good" and can serve as "root" elements * and which elements are "bad" and need support. Here, we will make the simple distinction:
[3]:
ci = CutInfo(mesh, gflset)
roots = ci.GetElementsOfType(NEG)
bads = ci.GetElementsOfType(IF)
print("bad element array: ", bads)
print("bad elements: ", [i for i in range(len(bads)) if bads[i]])
Draw(BitArrayCF(bads), mesh, "bad_elements")
bad element array: 0: 000011000011
bad elements: [4, 5, 10, 11]
[3]:
BaseWebGuiScene
ElementAggregation
¶
Next, we introduce a new class ElementAggregation
that organizes the 'neighborhood support':
Let’s have a brief look at its functionality:
[4]:
help(ElementAggregation)
Help on class ElementAggregation in module xfem.ngsxfem_py:
class ElementAggregation(pybind11_builtins.pybind11_object)
 ElementAggregation does the following:
 It collects patches of elements that allow to stabilize bad cut elements by at least one
 good element (the root element).
 )

 Method resolution order:
 ElementAggregation
 pybind11_builtins.pybind11_object
 builtins.object

 Methods defined here:

 Update(...)
 Update(self: xfem.ngsxfem_py.ElementAggregation, root_elements: pyngcore.pyngcore.BitArray, bad_elements: pyngcore.pyngcore.BitArray, heapsize: int = 1000000) > None


 Updates a Element Aggregation based ...

 __init__(...)
 __init__(self: xfem.ngsxfem_py.ElementAggregation, mesh: ngsolve.comp.Mesh, root_elements: object = None, bad_elements: object = None, heapsize: int = 1000000) > None


 Creates a ElementAggregation based on a mesh, a list of root and a list of bad elements.

 
 Readonly properties defined here:

 element_to_patch
 vector mapping elements to (nontrivial) patches

 els_in_nontrivial_patch
 BitArray that is true for every element that is part of a (nontrivial) patch

 els_in_trivial_patch
 BitArray that is true for every element that is not part of a (nontrivial) patch

 facet_to_patch
 vector mapping facets to (nontrivial) patches

 patch_interior_facets
 BitArray that is true for every facet that is *inside* an aggregation cluster

 
 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.
Let’s create an ElementAggregation
object and obtain a set of patches:
[5]:
EA = ElementAggregation(mesh, roots, bads)
patch_number_field = GridFunction(L2(mesh))
patch_number_field.vec.FV()[:] = EA.element_to_patch
Draw(patch_number_field, mesh, "patch_number_field", deformation=CF((0,0,0.02*patch_number_field)))
[5]:
BaseWebGuiScene
Obviously not all elements are in patches. Indeed the root elements that are not involved in supporting a bad element form "trivial patches":
[6]:
print("trivial root elements:", EA.els_in_trivial_patch)
print("nontrivial elements:", EA.els_in_nontrivial_patch)
trivial root elements: 0: 111000111000
nontrivial elements: 0: 000111000111
Similarly we can also relate facets that lie within a patch to the corresponding patch:
[7]:
print("EA.patch_interior_facets: ", EA.patch_interior_facets)
print("EA.facet_to_patch: ", EA.facet_to_patch)
EA.patch_interior_facets: 0: 00000001100000000110000
EA.facet_to_patch: [1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
So, now we have some nice patches. These can be used to tune ghost penalty terms (reduce the number of facets where GP is applied) or to change the FEM basis which we will discuss next.
Basis aggregation¶
On each of the patches we now want to get rid of the basis functions that are not supported on "root" elements. Let
be your usual development in your basis \(\{\varphi\}_i\). Now, we distinguish \(\mathcal{I}_{\text{root}}\) and \(\mathcal{I}_{\text{bad}}\), the index sets of basis functions with/without support on root elements. Then, for \(i \in \mathcal{I}_{\text{bad}}\) we want to express \(c_i\) as a (linear) function of \((c_j)_{j \in \mathcal{I}_{\text{root}}}\):
We start with
and obtain
Here \(\{\psi_i\}_{i \in \mathcal{I}_{\text{root}}}\) is the new basis obtained after aggregation. We will however characterize it always through the embedding matrix \(E\).
This can now be achieved by the method AggEmbedding
which yields the embedding matrix \(E\) (we explain how \(E\) is obtained below).
[8]:
help(AggEmbedding)
Help on function AggEmbedding in module xfem:
AggEmbedding(EA, fes, deformation=None, heapsize=1000000)
Computes and returns embedding matrix for a patchwise polynomial extension
(realized through ghost penalties),
followed by averaging if some dofs are shared by multiple patches.
Parameters
elagg : ElementAggregation
ElementAggregation instace defining the patches which are aggregated into a single element
fes : ngsolve.FESpace
The finite element space which is aggregated.
deformation : ngsolve.GridFunction [mesh.dim]
The mesh deformation (needed for Ghost penalty assembly)
heapsize : int
heapsize for local computations.
L2
low order¶
Let us try it out for the space of piecewise constants first:
[9]:
fes = L2(mesh, order=0)
E = AggEmbedding(EA, fes)
from helper import ShowPattern
ShowPattern(E)
print(E)
/usr/lib/python3/distpackages/scipy/__init__.py:146: UserWarning: A NumPy version >=1.17.3 and <1.25.0 is required for this version of SciPy (detected version 1.26.3
warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"
Row 0: 0: 1
Row 1: 1: 1
Row 2: 2: 1
Row 3: 3: 1
Row 4: 3: 1
Row 5: 3: 1
Row 6: 4: 1
Row 7: 5: 1
Row 8: 6: 1
Row 9: 7: 1
Row 10: 7: 1
Row 11: 7: 1
Let us take a look at the basis functions:
[10]:
def ExtendedBasisFunctionsAsMultiDim(EA, fes):
E = AggEmbedding(EA, fes)
gfshow = GridFunction(fes, multidim=0)
gf = GridFunction(fes)
coefvec = E.CreateRowVector()
for i in range(E.width):
coefvec[:] = 0
coefvec[i] = 1
gf.vec.data = E * coefvec
gfshow.AddMultiDimComponent(gf.vec)
return gfshow
[11]:
gfshow = ExtendedBasisFunctionsAsMultiDim(EA,fes)
Draw (gfshow, mesh, interpolate_multidim=False, animate=False, autoscale=True);
Use the multidimslider in the webgui to change the (interior) basis function and see the extended basis functions.
L2
higher order¶
Now, higher order:
[12]:
fes = L2(mesh, order=2)
ShowPattern(AggEmbedding(EA, fes))
gfshow = ExtendedBasisFunctionsAsMultiDim(EA,fes)
Draw (gfshow, mesh, interpolate_multidim=False, animate=False, autoscale=True);
H1
¶
When considering spaces with (partial) continuity, some dofs may appear in several patches but not on a root element. In this case the extensions are average:
[13]:
fes = H1(mesh, order=1)
ShowPattern(AggEmbedding(EA, fes))
gfshow = ExtendedBasisFunctionsAsMultiDim(EA,fes)
Draw (gfshow, mesh, interpolate_multidim=False, animate=False, autoscale=True);
Vectorvalued spaces¶
And now, HDiv
[14]:
fes = HDiv(mesh, order=1)
ShowPattern(AggEmbedding(EA, fes))
gfshow = ExtendedBasisFunctionsAsMultiDim(EA,fes)
Draw (gfshow, mesh, interpolate_multidim=False, animate=False, autoscale=True);
How it works¶
The mechanism behind the patchwise extension is very generic and obviously works for L2
, H1
, HDiv
, … in a generic way. But how?
A formulation as a harmonic extension¶
Let \(w_h(x) = \sum_{i \in \mathcal{I}_{\text{bad}}} c_i \varphi_i(x)\) be a finite element function only supported on bad elements and \(W_h\) the corresponding subspace of the underlying FE space \(V_h\) and \(X_h = W_h \setminus V_h\).
The idea is now that for every basis function \(\varphi_i\) with \(i \in \mathcal{I}_{\text{root}}\) we define its extension into \(W_h\) through
where \(\vert \cdot \vert_\ast\) is a smoothnessmeasuring seminorm. \(w_h\) is like a harmonic extension from the "root" into the "bad" elements.
A formulation as a patchwise harmonic extension¶
To localize this problem we proceed patchbypatch:
For every patch \(\omega\) we define \(W_h^{\omega}\) through functions that are supported on the patch, but not on a root element. Note that some basis function may appear in different patches now.
Then, we solve the local patch problems: For every basis function \(\varphi_i\) with \(i \in \mathcal{I}_{\text{root}}\) and support on the patch \(\omega\) we define its extension into \(W_h^\omega\) through
where \(\vert \cdot \vert_{\ast,\omega}\) is a smoothnessmeasuring seminorm.
The thusly obtained extended functions are afterwards averaged.
A generic seminorm \(\vert \varphi_i + w_h \vert_{\ast,\omega}\)¶
It now remains to define a proper seminorm.
When using AggEmbedding
as above the chosen seminorm is defined through a ghost penalty term inside all patches (Note that patchwise polynomials form the kernel of that operator):
def AggEmbedding(EA, fes, deformation=None, heapsize=1000000):
"""..."""
u,v = fes.TnT()
ghost_penalty = (u  u.Other()) * (v  v.Other()) * dFacetPatch(deformation=deformation)
return ExtensionEmbedding(EA, fes, ghost_penalty, heapsize=heapsize)
If you want to use a different seminorm for the extension, you can use ExtensionEmbedding
accordingly.
Applications of patch aggregations¶
We now have learned about how to form patches and extend basis functions. This can be exploited in several ways:
To solve CutFEM problems where the resulting system is reduced to the root unknowns by the embedding, cf. the
fictdom_aggfem.py
demo ofngsxfem
.Patches can be used to minimize the facets that are used in ghost penalty stabilizations, cf. e.g. unfmixed example.
Patches can be used to solve local patchwise problems, e.g. for postprocessings in unfitted mixed methods, cf. unfmixed example.