This page was generated from unit-2.10-dualbasis/dualbasis.ipynb.

2.10 Dual basis functions

We use dual basis functions to define interpolation operators, define transfer operators between different finite element spaces, and auxiliary space preconditioners.

Canonical interpolation

The canonical finite element interpolation operator is defined by specifying the degrees of freedom. For low order methods these are typically nodal values, while for high order methods these are most often moments. For example, the interpolation of a function \(u\) onto the \(p^{th}\) order triangle is given by: find \(u_{hp} \in V_{hp}\) such that

\begin{eqnarray*} u_{hp} (V) & = & u(V) \quad \forall \text{ vertices } V \\ \int_E u_{hp} q & = & \int_E u q \quad \forall q \in P^{p-2}(E) \; \forall \text{ edges } E \\ \int_T u_{hp} q & = & \int_T u q \quad \forall q \in P^{p-3}(T) \; \forall \text{ triangles } T \end{eqnarray*}

from ngsolve import *
from ngsolve.webgui import Draw

import matplotlib.pylab as plt
mesh = Mesh(unit_square.GenerateMesh(maxh=2))

The NGSolve 'Set' function does local projection, and simple averaging. In particular, this does not respect point values in mesh vertices.

fes = H1(mesh, order=3, low_order_space=False)

func = x*x*x*x
gfu = GridFunction(fes)
gfu.Set(func)
Draw (gfu)
print (gfu.vec)

 -0.0223792
 0.991843
 0.977621
 -0.00815655
 3.34984
 2.19912
 3.34984
       2
 0.00660131
 1.39013e-14
 0.00660131
 -3.19744e-14
 3.34984
 -1.80088
 -2.17601
 1.82399

Most NGSolve finite element spaces provide now a "dual" operator, which delivers the moments (i.e. the dual space basis functions) instead of function values. The integrals over faces, edges and also vertices are defined by co-dimension 1 (=BND), co-dimension 2 (=BBND) or co-dimension 3 (=BBBND) integrals over the volume elements. We define a variational problem for canonical interpolation:

u,v = fes.TnT()
vdual = v.Operator("dual")

a = BilinearForm(fes)
a += u*vdual*dx + u*vdual*dx(element_vb=BND) + \
    u*vdual*dx(element_vb=BBND)
a.Assemble()

f = LinearForm(fes)
f += func*vdual*dx + func*vdual*dx(element_vb=BND) + \
    func*vdual*dx(element_vb=BBND)
f.Assemble()

# interpolation in vertices preserves values 0 and 1
gfu.vec.data = a.mat.Inverse() * f.vec
print (gfu.vec)
Draw (gfu);

       0
       1
       1
       0
     3.6
       2
     3.6
       2
       0
       0
 -6.93335e-32
 -1.66533e-15
     3.6
      -2
      -2
       2

The vertex degrees of freedom vanish for edge and element basis functions, and the edge degrees of freedom vanish for element basis functions, but not vice-versa. Thus, the obtained matrix A is block-triangular:

import scipy.sparse as sp
A = sp.csr_matrix(a.mat.CSR())
plt.rcParams['figure.figsize'] = (4,4)
plt.spy(A)
plt.show()

We can use proper block Gauss-Seidel smoothing for solving with that block triangular matrix by blocking the dofs for the individual vertices, edges and elements. Since the NGSolve Gauss-Seidel smoother reorders the order of smoothing blocks for parallelization, we have to take care to first compute vertex values, then edge values, and finally element values by running three different Gauss-Seidel sweeps.

vblocks = [fes.GetDofNrs(vertex) for vertex in mesh.vertices]
eblocks = [fes.GetDofNrs(edge) for edge in mesh.edges]
fblocks = [fes.GetDofNrs(face) for face in mesh.faces]

print (vblocks)
print (eblocks)
print (fblocks)

vinv = a.mat.CreateBlockSmoother(vblocks)
einv = a.mat.CreateBlockSmoother(eblocks)
finv = a.mat.CreateBlockSmoother(fblocks)

vinv.Smooth(gfu.vec, f.vec)
einv.Smooth(gfu.vec, f.vec)
finv.Smooth(gfu.vec, f.vec)
print (gfu.vec)

[(0,), (1,), (2,), (3,)]
[(4, 5), (6, 7), (8, 9), (10, 11), (12, 13)]
[(14,), (15,)]
       0
       1
       1
       0
     3.6
       2
     3.6
       2
       0
       0
 -6.93335e-32
 -1.66533e-15
     3.6
      -2
      -2
       2

Embedding Finite Element Spaces

This interpolation can be used to transform functions from one finite element space \(V_{src}\) to another one \(V_{dst}\). We use the dual space of the destination space:

\[\int_{node} u_{dst} v_{dual} = \int_{node} u_{src} v_{dual} \qquad \forall \, v_{dual} \; \forall \, \text{nodes}\]

The left hand side leads to a non-symmetric square matrix, the right hand side to a rectangular matrix.

As an example we implement the transformation from an vector valued \(H^1\) space into \(H(\operatorname{div})\):

from ngsolve import *
from ngsolve.webgui import Draw
from netgen.geom2d import unit_square
mesh = Mesh(unit_square.GenerateMesh(maxh=0.2))

fesh1 = VectorH1(mesh, order=2)
feshdiv = HDiv(mesh, order=2)

gfuh1 = GridFunction(fesh1)
gfuh1.Set ( (x*x,y*y) )

gfuhdiv = GridFunction(feshdiv, name="uhdiv")

Build the matrices, and use a direct solver:

amixed = BilinearForm(trialspace=fesh1, testspace=feshdiv)
ahdiv = BilinearForm(feshdiv)

u,v = feshdiv.TnT()
vdual = v.Operator("dual")
uh1 = fesh1.TrialFunction()

dS = dx(element_boundary=True)
ahdiv += u*vdual * dx + u*vdual * dS
ahdiv.Assemble()

amixed += uh1*vdual*dx + uh1*vdual*dS
amixed.Assemble()

# build transformation operator:
transform = ahdiv.mat.Inverse() @ amixed.mat
gfuhdiv.vec.data = transform * gfuh1.vec

Draw (gfuh1)
Draw (gfuhdiv)

BaseWebGuiScene

We implement a linear operator performing the fast conversion by Gauss-Seidel smoothing:

class MyBlockInverse(BaseMatrix):
    def __init__ (self, mat, eblocks, fblocks):
        super(MyBlockInverse, self).__init__()
        self.mat = mat
        self.einv = mat.CreateBlockSmoother(eblocks)
        self.finv = mat.CreateBlockSmoother(fblocks)
        self.res = self.mat.CreateColVector()

    def CreateRowVector(self):
        return self.mat.CreateColVector()
    def CreateColVector(self):
        return self.mat.CreateRowVector()

    def Mult(self, x, y):
        # y[:] = 0
        # self.einv.Smooth(y,x)    #   y = y +  A_E^-1  (x - A y)
        # self.finv.Smooth(y,x)    #   y = y +  A_E^-1  (x - A y)

        # the same, but we see how to transpose that
        y.data = self.einv * x
        self.res.data = x - self.mat * y
        y.data += finv * self.res

    def MultTrans(self, x, y):
        y.data = self.finv.T * x
        self.res.data = x - self.mat.T * y
        y.data += einv.T * self.res


eblocks = [feshdiv.GetDofNrs(edge) for edge in mesh.edges]
fblocks = [feshdiv.GetDofNrs(face) for face in mesh.faces]

transform = MyBlockInverse(ahdiv.mat, eblocks, fblocks) @ amixed.mat
gfuhdiv.vec.data = transform * gfuh1.vec

Auxiliary Space Preconditioning

Nepomnyaschikh 91, Hiptmair-Xu 07, ….

Assume we have a complicated problem with some complicated discretization, and we have good preconditioners for a nodal finite element discretization for the Laplace operator. By auxiliary space preconditioning we can reuse the simple preconditioners for the complicated problems. It is simple, and works well in many cases.

As a simple example, we precondition a DG discretization by an \(H^1\) conforming method.

mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))

The DG discretization:

order=4
fesDG = L2(mesh, order=order, dgjumps=True)
u,v = fesDG.TnT()
aDG = BilinearForm(fesDG)
jump_u = u-u.Other()
jump_v = v-v.Other()
n = specialcf.normal(2)
mean_dudn = 0.5*n * (grad(u)+grad(u.Other()))
mean_dvdn = 0.5*n * (grad(v)+grad(v.Other()))
alpha = 4
h = specialcf.mesh_size
aDG = BilinearForm(fesDG)
aDG += grad(u)*grad(v) * dx
aDG += alpha*order**2/h*jump_u*jump_v * dx(skeleton=True)
aDG += alpha*order**2/h*u*v * ds(skeleton=True)
aDG += (-mean_dudn*jump_v -mean_dvdn*jump_u) * dx(skeleton=True)
aDG += (-n*grad(u)*v-n*grad(v)*u)* ds(skeleton=True)
aDG.Assemble()

fDG = LinearForm(fesDG)
fDG += 1*v * dx
fDG.Assemble()
gfuDG = GridFunction(fesDG)

The auxiliary \(H^1\) discretization:

fesH1 = H1(mesh, order=2, dirichlet=".*")
u,v = fesH1.TnT()
aH1 = BilinearForm(fesH1)
aH1 += grad(u)*grad(v)*dx
preH1 = Preconditioner(aH1, "bddc")
aH1.Assemble()

<ngsolve.comp.BilinearForm at 0x7f9f05f84eb0>

transform = fesH1.ConvertL2Operator(fesDG)
pre = transform @ preH1.mat @ transform.T + aDG.mat.CreateSmoother()

solvers.CG(mat=aDG.mat, rhs=fDG.vec, sol=gfuDG.vec, pre=pre, \
           maxsteps=200)

Draw (gfuDG)

CG iteration 1, residual = 0.2532119733219237
CG iteration 2, residual = 0.03283242148303869
CG iteration 3, residual = 0.008130862525681015
CG iteration 4, residual = 0.0034342488110603643
CG iteration 5, residual = 0.0011783456811530293
CG iteration 6, residual = 0.0009046659695335783
CG iteration 7, residual = 0.0009161002088522016
CG iteration 8, residual = 0.0010697558880604528
CG iteration 9, residual = 0.0007971399933259626
CG iteration 10, residual = 0.000662789636667445
CG iteration 11, residual = 0.0005095344857659155
CG iteration 12, residual = 0.00042466444937136347
CG iteration 13, residual = 0.00032617023029280325
CG iteration 14, residual = 0.00034395529498088624
CG iteration 15, residual = 0.00027105251267776236
CG iteration 16, residual = 0.00023336881656180314
CG iteration 17, residual = 0.0001798377168164044
CG iteration 18, residual = 0.0001356676383204489
CG iteration 19, residual = 8.969965836308495e-05
CG iteration 20, residual = 7.588906438782797e-05
CG iteration 21, residual = 6.938043743084663e-05
CG iteration 22, residual = 6.381887596483748e-05
CG iteration 23, residual = 5.180208019108754e-05
CG iteration 24, residual = 4.3817631293678964e-05
CG iteration 25, residual = 3.640724987679596e-05
CG iteration 26, residual = 2.980836083285286e-05
CG iteration 27, residual = 2.9223985181308503e-05
CG iteration 28, residual = 2.9089520392432463e-05
CG iteration 29, residual = 2.6390001567259255e-05
CG iteration 30, residual = 2.0442446246668556e-05
CG iteration 31, residual = 1.7799471861373654e-05
CG iteration 32, residual = 1.428035397673256e-05
CG iteration 33, residual = 1.172408643444086e-05
CG iteration 34, residual = 1.03990820923482e-05
CG iteration 35, residual = 9.357450254587714e-06
CG iteration 36, residual = 8.04613249625652e-06
CG iteration 37, residual = 6.434972748061881e-06
CG iteration 38, residual = 4.977423074804273e-06
CG iteration 39, residual = 4.189516701918544e-06
CG iteration 40, residual = 3.879825859447426e-06
CG iteration 41, residual = 3.424994163355379e-06
CG iteration 42, residual = 3.1039505254082392e-06
CG iteration 43, residual = 2.6848114735438076e-06
CG iteration 44, residual = 2.1289592891352863e-06
CG iteration 45, residual = 1.7175875308400938e-06
CG iteration 46, residual = 1.5033515289259946e-06
CG iteration 47, residual = 1.4944954266968019e-06
CG iteration 48, residual = 1.4120180027628892e-06
CG iteration 49, residual = 1.1860668787092936e-06
CG iteration 50, residual = 9.38436483757009e-07
CG iteration 51, residual = 7.855316166902367e-07
CG iteration 52, residual = 6.514966222954281e-07
CG iteration 53, residual = 5.847843124222843e-07
CG iteration 54, residual = 5.642748143010218e-07
CG iteration 55, residual = 4.856169352654797e-07
CG iteration 56, residual = 3.6490360544844745e-07
CG iteration 57, residual = 2.8732276253075304e-07
CG iteration 58, residual = 2.5242710520942045e-07
CG iteration 59, residual = 2.1120183910235296e-07
CG iteration 60, residual = 1.8157945371825142e-07
CG iteration 61, residual = 1.6765792567394388e-07
CG iteration 62, residual = 1.353430987382415e-07
CG iteration 63, residual = 1.0740625460955306e-07
CG iteration 64, residual = 8.692499337660225e-08
CG iteration 65, residual = 6.99102464203828e-08
CG iteration 66, residual = 6.16682084855965e-08
CG iteration 67, residual = 5.426467057609486e-08
CG iteration 68, residual = 4.492616630467071e-08
CG iteration 69, residual = 3.427094557093642e-08
CG iteration 70, residual = 2.736795235873319e-08
CG iteration 71, residual = 2.2107937585507486e-08
CG iteration 72, residual = 1.868591295158383e-08
CG iteration 73, residual = 1.7316902148540638e-08
CG iteration 74, residual = 1.5508098318254025e-08
CG iteration 75, residual = 1.3238529365902683e-08
CG iteration 76, residual = 9.994763302499639e-09
CG iteration 77, residual = 7.8833133404764e-09
CG iteration 78, residual = 6.7873917080074224e-09
CG iteration 79, residual = 6.238631290597011e-09
CG iteration 80, residual = 5.761549944097009e-09
CG iteration 81, residual = 4.886605172386132e-09
CG iteration 82, residual = 3.7486015148701065e-09
CG iteration 83, residual = 2.7484507255154553e-09
CG iteration 84, residual = 2.280167928231959e-09
CG iteration 85, residual = 2.0226924556861182e-09
CG iteration 86, residual = 1.8702257796507456e-09
CG iteration 87, residual = 1.5609302585058095e-09
CG iteration 88, residual = 1.2269789259839587e-09
CG iteration 89, residual = 9.112910421711045e-10
CG iteration 90, residual = 7.127268960570408e-10
CG iteration 91, residual = 5.937551150270976e-10
CG iteration 92, residual = 5.51420345007548e-10
CG iteration 93, residual = 4.998776637894567e-10
CG iteration 94, residual = 4.1240181286042545e-10
CG iteration 95, residual = 3.2669176027200394e-10
CG iteration 96, residual = 2.599581160147041e-10
CG iteration 97, residual = 2.1815248859947405e-10
CG iteration 98, residual = 2.0149551888013363e-10
CG iteration 99, residual = 1.8834681060520743e-10
CG iteration 100, residual = 1.6606780443995193e-10
CG iteration 101, residual = 1.3400679424698756e-10
CG iteration 102, residual = 1.0597986884306437e-10
CG iteration 103, residual = 8.29949236055755e-11
CG iteration 104, residual = 7.177342176620354e-11
CG iteration 105, residual = 6.589102789283192e-11
CG iteration 106, residual = 5.88909580403242e-11
CG iteration 107, residual = 4.7475559498018766e-11
CG iteration 108, residual = 3.831939163542913e-11
CG iteration 109, residual = 2.981978886185597e-11
CG iteration 110, residual = 2.529741938897603e-11
CG iteration 111, residual = 2.4025898588910242e-11
CG iteration 112, residual = 2.3189757188373442e-11
CG iteration 113, residual = 2.0566658349221437e-11
CG iteration 114, residual = 1.6123922268873695e-11
CG iteration 115, residual = 1.2686047253736232e-11
CG iteration 116, residual = 1.0442133746666682e-11
CG iteration 117, residual = 1.0074235133755179e-11
CG iteration 118, residual = 9.313960121844095e-12
CG iteration 119, residual = 8.444022057968894e-12
CG iteration 120, residual = 7.023516296999095e-12
CG iteration 121, residual = 5.340599899207211e-12
CG iteration 122, residual = 4.201138029957767e-12
CG iteration 123, residual = 3.832348860032032e-12
CG iteration 124, residual = 3.6215301967497565e-12
CG iteration 125, residual = 3.0685172154425e-12
CG iteration 126, residual = 2.462462578966119e-12
CG iteration 127, residual = 1.8937381804356026e-12
CG iteration 128, residual = 1.5445175514430808e-12
CG iteration 129, residual = 1.3365865282386032e-12
CG iteration 130, residual = 1.2596961871684347e-12
CG iteration 131, residual = 1.124637848451851e-12
CG iteration 132, residual = 9.044195430190062e-13
CG iteration 133, residual = 6.755352781093621e-13
CG iteration 134, residual = 5.478309620646218e-13
CG iteration 135, residual = 4.587901445406802e-13
CG iteration 136, residual = 4.207573217764919e-13
CG iteration 137, residual = 3.8906948929441636e-13
CG iteration 138, residual = 3.4572985303551005e-13
CG iteration 139, residual = 2.659787680410384e-13
CG iteration 140, residual = 2.1117431987616406e-13

BaseWebGuiScene