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 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
 3.06756e-15
 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()

../../_images/i-tutorials_unit-2.10-dualbasis_dualbasis_8_0.png

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()
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:

[11]:

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()

[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.25284243668046297
iteration 1 error = 0.032978067384423895
iteration 2 error = 0.008565533056584225
iteration 3 error = 0.003350791073889276
iteration 4 error = 0.001235225038618856
iteration 5 error = 0.000895302737947024
iteration 6 error = 0.0009408102631568238
iteration 7 error = 0.0010428541678774179
iteration 8 error = 0.0008241218795830562
iteration 9 error = 0.0006517325370372178
iteration 10 error = 0.0005341681600736738
iteration 11 error = 0.0004077364183234147
iteration 12 error = 0.00035057733041789403
iteration 13 error = 0.00033153548156673505
iteration 14 error = 0.00028803580573035615
iteration 15 error = 0.0002312154589429102
iteration 16 error = 0.00018655018570322998
iteration 17 error = 0.00013157695014631517
iteration 18 error = 9.230174022085532e-05
iteration 19 error = 7.485395694726418e-05
iteration 20 error = 7.207046355389521e-05
iteration 21 error = 5.9713438908785194e-05
iteration 22 error = 4.9829695090184656e-05
iteration 23 error = 4.120697126467163e-05
iteration 24 error = 3.41764899033569e-05
iteration 25 error = 2.9765047174175662e-05
iteration 26 error = 3.0227068897520353e-05
iteration 27 error = 2.9571437689583936e-05
iteration 28 error = 2.5690373149739806e-05
iteration 29 error = 2.020293864722904e-05
iteration 30 error = 1.7555484171470644e-05
iteration 31 error = 1.3996042684636183e-05
iteration 32 error = 1.1567524127370245e-05
iteration 33 error = 1.0573917718227942e-05
iteration 34 error = 9.869130802509562e-06
iteration 35 error = 7.968202828295557e-06
iteration 36 error = 6.502761855274084e-06
iteration 37 error = 4.734565244528577e-06
iteration 38 error = 4.169406582627728e-06
iteration 39 error = 3.740963548260426e-06
iteration 40 error = 3.5148808051781854e-06
iteration 41 error = 3.050386409022662e-06
iteration 42 error = 2.7094137652647153e-06
iteration 43 error = 2.07655204056756e-06
iteration 44 error = 1.6480108115371886e-06
iteration 45 error = 1.47467817278642e-06
iteration 46 error = 1.4114503319115066e-06
iteration 47 error = 1.4006035122630387e-06
iteration 48 error = 1.1481091661340297e-06
iteration 49 error = 9.279751292280902e-07
iteration 50 error = 7.684324049143166e-07
iteration 51 error = 6.343737423517505e-07
iteration 52 error = 5.868300135075061e-07
iteration 53 error = 5.483251077318047e-07
iteration 54 error = 4.601522921980794e-07
iteration 55 error = 3.5654272425279007e-07
iteration 56 error = 2.793730241945098e-07
iteration 57 error = 2.466033071081388e-07
iteration 58 error = 2.0116576888680717e-07
iteration 59 error = 1.7441528258855818e-07
iteration 60 error = 1.552070033063872e-07
iteration 61 error = 1.2841858136453518e-07
iteration 62 error = 9.985763939316765e-08
iteration 63 error = 8.08856726025781e-08
iteration 64 error = 6.482197424796254e-08
iteration 65 error = 5.788278775485292e-08
iteration 66 error = 5.155782466351719e-08
iteration 67 error = 4.207825530042541e-08
iteration 68 error = 3.349444609875453e-08
iteration 69 error = 2.68486711463308e-08
iteration 70 error = 2.142356403620811e-08
iteration 71 error = 1.8330570047070182e-08
iteration 72 error = 1.7353354273815782e-08
iteration 73 error = 1.518534436997985e-08
iteration 74 error = 1.3273851059294744e-08
iteration 75 error = 1.0106789943533343e-08
iteration 76 error = 7.987735063386033e-09
iteration 77 error = 6.969750611404846e-09
iteration 78 error = 6.5144315988776994e-09
iteration 79 error = 6.059090244119175e-09
iteration 80 error = 4.875319110184188e-09
iteration 81 error = 3.861739632308879e-09
iteration 82 error = 2.9665506559001218e-09
iteration 83 error = 2.3995608341140846e-09
iteration 84 error = 2.1812926073426197e-09
iteration 85 error = 2.036305947459066e-09
iteration 86 error = 1.6651955286634425e-09
iteration 87 error = 1.283258478782406e-09
iteration 88 error = 9.670470794800201e-10
iteration 89 error = 7.828922434338418e-10
iteration 90 error = 6.674761986226291e-10
iteration 91 error = 6.03412481370858e-10
iteration 92 error = 5.644285777584519e-10
iteration 93 error = 4.4825552419084787e-10
iteration 94 error = 3.4440585189783814e-10
iteration 95 error = 2.822365078350753e-10
iteration 96 error = 2.488561253830576e-10
iteration 97 error = 2.3078260347537464e-10
iteration 98 error = 2.1689142905582958e-10
iteration 99 error = 1.86402783155493e-10
iteration 100 error = 1.4801015185106908e-10
iteration 101 error = 1.2201651732794116e-10
iteration 102 error = 1.0348297931535949e-10
iteration 103 error = 9.479983019879548e-11
iteration 104 error = 8.721469338934802e-11
iteration 105 error = 7.443117940682301e-11
iteration 106 error = 6.280343850514982e-11
iteration 107 error = 5.017380130592685e-11
iteration 108 error = 3.921995089176934e-11
iteration 109 error = 3.4729006072585804e-11
iteration 110 error = 3.211177863537385e-11
iteration 111 error = 2.922151599403755e-11
iteration 112 error = 2.3831039071107446e-11
iteration 113 error = 1.8027585299927186e-11
iteration 114 error = 1.4199876699663712e-11
iteration 115 error = 1.2486551276030442e-11
iteration 116 error = 1.1945924860564022e-11
iteration 117 error = 1.1170399066666007e-11
iteration 118 error = 9.282960246017481e-12
iteration 119 error = 7.271561012586448e-12
iteration 120 error = 5.4419894690663774e-12
iteration 121 error = 4.599084570639826e-12
iteration 122 error = 4.26034528452037e-12
iteration 123 error = 3.81806500318616e-12
iteration 124 error = 3.2222286423260747e-12
iteration 125 error = 2.6082196447179107e-12
iteration 126 error = 2.0409992913604886e-12
iteration 127 error = 1.6285185376177415e-12
iteration 128 error = 1.4603556589926205e-12
iteration 129 error = 1.3226999283113427e-12
iteration 130 error = 1.1433191322070923e-12
iteration 131 error = 9.044682421525567e-13
iteration 132 error = 6.630037139724958e-13
iteration 133 error = 5.362345450065622e-13
iteration 134 error = 4.733717009042433e-13
iteration 135 error = 4.281852587748644e-13
iteration 136 error = 3.802497903539686e-13
iteration 137 error = 3.193450616643783e-13
iteration 138 error = 2.3783196413483784e-13

[ ]: