# 2.10 Dual basis functions¶

We use dual basis functions to define interpolation operators, define transfer operators between different finite element spaces, and auxiliar 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*}

[1]:

import netgen.gui
from ngsolve import *
from netgen.geom2d import unit_square
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.

[2]:

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
8.91822e-18
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 interpolaion:

[3]:

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 vise versa. Thus, the obtained matrix A is block-triangular:

[4]:

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.

[5]:

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})$$:

[6]:

import netgen.gui
from ngsolve import *
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:

[7]:

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)


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

[8]:

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.

[9]:

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


The DG discretization:

[10]:

order=4
fesDG = L2(mesh, order=order, dgjumps=True)
u,v = fesDG.TnT()
jump_u = u-u.Other()
jump_v = v-v.Other()
n = specialcf.normal(2)
alpha = 4
h = specialcf.mesh_size
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)

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


The auxiliary $$H^1$$ discretization:

[11]:

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

[11]:

<ngsolve.comp.BilinearForm at 0x7f18277209b0>

[12]:

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)

iteration 0 error = 0.25330526037462214
iteration 1 error = 0.033096246037648165
iteration 2 error = 0.008611208192688617
iteration 3 error = 0.004161316791286654
iteration 4 error = 0.0014153167752984758
iteration 5 error = 0.0010294082919071873
iteration 6 error = 0.0009409975524653755
iteration 7 error = 0.001115058802510059
iteration 8 error = 0.0009072553532028597
iteration 9 error = 0.0007005597437154765
iteration 10 error = 0.0005825091563469165
iteration 11 error = 0.00045806716363981236
iteration 12 error = 0.0003501358759117548
iteration 13 error = 0.0003389579426000714
iteration 14 error = 0.0003009288036336665
iteration 15 error = 0.00024917896162954224
iteration 16 error = 0.00020393153941205143
iteration 17 error = 0.00015750033322040478
iteration 18 error = 0.00011665450896384623
iteration 19 error = 8.69386790084416e-05
iteration 20 error = 7.785593301103237e-05
iteration 21 error = 6.884157007377474e-05
iteration 22 error = 5.580418759769967e-05
iteration 23 error = 4.6587301525775e-05
iteration 24 error = 3.993871204144851e-05
iteration 25 error = 3.268865081952109e-05
iteration 26 error = 2.9155421446685284e-05
iteration 27 error = 2.895026601283125e-05
iteration 28 error = 2.715175136737116e-05
iteration 29 error = 2.2894435672077907e-05
iteration 30 error = 1.929409730319834e-05
iteration 31 error = 1.714430095995899e-05
iteration 32 error = 1.2932552421849758e-05
iteration 33 error = 1.1175658200042554e-05
iteration 34 error = 9.991739671222712e-06
iteration 35 error = 9.63849135821357e-06
iteration 36 error = 7.788256388277832e-06
iteration 37 error = 5.8846383890814665e-06
iteration 38 error = 4.82025641595475e-06
iteration 39 error = 4.147719419989545e-06
iteration 40 error = 3.6572807096711626e-06
iteration 41 error = 3.528176840424625e-06
iteration 42 error = 3.361345329993807e-06
iteration 43 error = 2.639834599746446e-06
iteration 44 error = 2.052806670695037e-06
iteration 45 error = 1.6586528592920206e-06
iteration 46 error = 1.555796207356698e-06
iteration 47 error = 1.5199816456631496e-06
iteration 48 error = 1.35861054399637e-06
iteration 49 error = 1.1159673333244119e-06
iteration 50 error = 8.810515092816279e-07
iteration 51 error = 7.223412178531473e-07
iteration 52 error = 6.222439483826799e-07
iteration 53 error = 5.869599046990032e-07
iteration 54 error = 5.271249628074154e-07
iteration 55 error = 4.276420883873427e-07
iteration 56 error = 3.466766571161947e-07
iteration 57 error = 2.7731147533246614e-07
iteration 58 error = 2.3651665434280193e-07
iteration 59 error = 1.9307560915483528e-07
iteration 60 error = 1.7433235630800746e-07
iteration 61 error = 1.4824475272995236e-07
iteration 62 error = 1.2078504304624718e-07
iteration 63 error = 9.28728408355508e-08
iteration 64 error = 7.362507051471966e-08
iteration 65 error = 6.104834176589066e-08
iteration 66 error = 5.802300802344468e-08
iteration 67 error = 4.7686280370093355e-08
iteration 68 error = 3.972120957369619e-08
iteration 69 error = 3.242850428456667e-08
iteration 70 error = 2.672991014523919e-08
iteration 71 error = 2.1444242344748922e-08
iteration 72 error = 1.7927730338148645e-08
iteration 73 error = 1.7401401742927986e-08
iteration 74 error = 1.4431907094196668e-08
iteration 75 error = 1.2301070308081875e-08
iteration 76 error = 9.139113364497722e-09
iteration 77 error = 7.377883337218838e-09
iteration 78 error = 6.198580256173309e-09
iteration 79 error = 5.734130490194182e-09
iteration 80 error = 5.481851034254635e-09
iteration 81 error = 4.494447993738023e-09
iteration 82 error = 3.416685458671028e-09
iteration 83 error = 2.651203429481029e-09
iteration 84 error = 2.3212227045637192e-09
iteration 85 error = 2.0741908219103763e-09
iteration 86 error = 1.973305602160263e-09
iteration 87 error = 1.657600194379589e-09
iteration 88 error = 1.227955190605025e-09
iteration 89 error = 9.174882763750405e-10
iteration 90 error = 7.771196280810747e-10
iteration 91 error = 6.922557345198841e-10
iteration 92 error = 6.512710135319149e-10
iteration 93 error = 5.651786195219746e-10
iteration 94 error = 4.5071510745810753e-10
iteration 95 error = 3.251276414483687e-10
iteration 96 error = 2.6747000363132605e-10
iteration 97 error = 2.556402841988816e-10
iteration 98 error = 2.379549114420168e-10
iteration 99 error = 2.05919292660053e-10
iteration 100 error = 1.6036235970562304e-10
iteration 101 error = 1.2318690391055095e-10
iteration 102 error = 1.0109606102873154e-10
iteration 103 error = 9.240005489030931e-11
iteration 104 error = 8.836556544410632e-11
iteration 105 error = 7.445129972839662e-11
iteration 106 error = 5.71138323708523e-11
iteration 107 error = 4.152798594663944e-11
iteration 108 error = 3.482125912662574e-11
iteration 109 error = 3.0543563042063194e-11
iteration 110 error = 2.7745697492572214e-11
iteration 111 error = 2.4981697174724386e-11
iteration 112 error = 1.981582959444435e-11
iteration 113 error = 1.5340965266378727e-11
iteration 114 error = 1.2343606348534228e-11
iteration 115 error = 1.0419630326857591e-11
iteration 116 error = 9.967037319381216e-12
iteration 117 error = 8.973603105462785e-12
iteration 118 error = 7.266034945591235e-12
iteration 119 error = 6.209644880962063e-12
iteration 120 error = 5.0656059023065446e-12
iteration 121 error = 4.044567068713039e-12
iteration 122 error = 3.5708504228421606e-12
iteration 123 error = 3.1516566586023165e-12
iteration 124 error = 2.534621059646172e-12
iteration 125 error = 2.135860425297587e-12
iteration 126 error = 1.7391017505428728e-12
iteration 127 error = 1.4431246503619532e-12
iteration 128 error = 1.2198496303883594e-12
iteration 129 error = 1.088262137300684e-12
iteration 130 error = 9.928308280412695e-13
iteration 131 error = 8.722434017263604e-13
iteration 132 error = 7.291963083719915e-13
iteration 133 error = 5.818089656333564e-13
iteration 134 error = 4.434039832248252e-13
iteration 135 error = 3.7433704729608704e-13
iteration 136 error = 3.426320221014901e-13
iteration 137 error = 3.179083168870281e-13
iteration 138 error = 2.880465065459007e-13
iteration 139 error = 2.311347703446632e-13

[ ]: