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]:

from ngsolve import *
from ngsolve.webgui import Draw
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

[3]:

BaseWebGuiScene

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]:

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:

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

[7]:

BaseWebGuiScene

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

[11]:

<ngsolve.comp.BilinearForm at 0x7ffba91d26f0>

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

WARNING: maxsteps is deprecated, use maxiter instead!
CG iteration 1, residual = 0.25330526037462214
CG iteration 2, residual = 0.03309624603764793
CG iteration 3, residual = 0.008611208192688572
CG iteration 4, residual = 0.004161316791286645
CG iteration 5, residual = 0.0014153167752984812
CG iteration 6, residual = 0.0010294082919071756
CG iteration 7, residual = 0.0009409975524653673
CG iteration 8, residual = 0.0011150588025100386
CG iteration 9, residual = 0.0009072553532028417
CG iteration 10, residual = 0.0007005597437154712
CG iteration 11, residual = 0.0005825091563469057
CG iteration 12, residual = 0.00045806716363980797
CG iteration 13, residual = 0.000350135875911749
CG iteration 14, residual = 0.00033895794260006715
CG iteration 15, residual = 0.0003009288036336635
CG iteration 16, residual = 0.00024917896162954045
CG iteration 17, residual = 0.00020393153941204994
CG iteration 18, residual = 0.00015750033322040432
CG iteration 19, residual = 0.00011665450896384501
CG iteration 20, residual = 8.693867900844134e-05
CG iteration 21, residual = 7.785593301103119e-05
CG iteration 22, residual = 6.88415700737739e-05
CG iteration 23, residual = 5.580418759769907e-05
CG iteration 24, residual = 4.658730152577426e-05
CG iteration 25, residual = 3.9938712041447935e-05
CG iteration 26, residual = 3.268865081952005e-05
CG iteration 27, residual = 2.9155421446684847e-05
CG iteration 28, residual = 2.8950266012830194e-05
CG iteration 29, residual = 2.7151751367370894e-05
CG iteration 30, residual = 2.28944356720773e-05
CG iteration 31, residual = 1.929409730319816e-05
CG iteration 32, residual = 1.71443009599587e-05
CG iteration 33, residual = 1.2932552421849839e-05
CG iteration 34, residual = 1.1175658200042496e-05
CG iteration 35, residual = 9.991739671222797e-06
CG iteration 36, residual = 9.638491358213626e-06
CG iteration 37, residual = 7.788256388277974e-06
CG iteration 38, residual = 5.8846383890815114e-06
CG iteration 39, residual = 4.82025641595477e-06
CG iteration 40, residual = 4.1477194199895195e-06
CG iteration 41, residual = 3.6572807096711096e-06
CG iteration 42, residual = 3.5281768404246238e-06
CG iteration 43, residual = 3.361345329993782e-06
CG iteration 44, residual = 2.639834599746442e-06
CG iteration 45, residual = 2.0528066706949604e-06
CG iteration 46, residual = 1.6586528592918764e-06
CG iteration 47, residual = 1.5557962073562954e-06
CG iteration 48, residual = 1.51998164566208e-06
CG iteration 49, residual = 1.3586105439931027e-06
CG iteration 50, residual = 1.1159673333146486e-06
CG iteration 51, residual = 8.810515092619467e-07
CG iteration 52, residual = 7.223412178036316e-07
CG iteration 53, residual = 6.222439482912041e-07
CG iteration 54, residual = 5.869599045065223e-07
CG iteration 55, residual = 5.271249619744767e-07
CG iteration 56, residual = 4.2764208617633895e-07
CG iteration 57, residual = 3.466766520287763e-07
CG iteration 58, residual = 2.773114643857253e-07
CG iteration 59, residual = 2.3651662651386224e-07
CG iteration 60, residual = 1.930755493486055e-07
CG iteration 61, residual = 1.743321874966022e-07
CG iteration 62, residual = 1.482442349947516e-07
CG iteration 63, residual = 1.207839093379714e-07
CG iteration 64, residual = 9.286997960336695e-08
CG iteration 65, residual = 7.361889760148165e-08
CG iteration 66, residual = 6.103213218390078e-08
CG iteration 67, residual = 5.798413753305417e-08
CG iteration 68, residual = 4.754288830696697e-08
CG iteration 69, residual = 3.934273045630499e-08
CG iteration 70, residual = 3.1532701808006754e-08
CG iteration 71, residual = 2.537986158590304e-08
CG iteration 72, residual = 1.9777168390646605e-08
CG iteration 73, residual = 1.7233670623280566e-08
CG iteration 74, residual = 1.6524970268472126e-08
CG iteration 75, residual = 1.39422842167317e-08
CG iteration 76, residual = 1.2119311914285624e-08
CG iteration 77, residual = 9.146462581493334e-09
CG iteration 78, residual = 7.523297824882635e-09
CG iteration 79, residual = 6.263353050416141e-09
CG iteration 80, residual = 5.765856431747356e-09
CG iteration 81, residual = 5.490605095372075e-09
CG iteration 82, residual = 4.496038421429952e-09
CG iteration 83, residual = 3.4169535462511374e-09
CG iteration 84, residual = 2.6512985782540083e-09
CG iteration 85, residual = 2.3212494299601334e-09
CG iteration 86, residual = 2.0741985402427964e-09
CG iteration 87, residual = 1.9733072468390902e-09
CG iteration 88, residual = 1.6576005213733278e-09
CG iteration 89, residual = 1.2279552732362014e-09
CG iteration 90, residual = 9.174882947575928e-10
CG iteration 91, residual = 7.771196346521387e-10
CG iteration 92, residual = 6.922557367473674e-10
CG iteration 93, residual = 6.512710143145366e-10
CG iteration 94, residual = 5.65178619653505e-10
CG iteration 95, residual = 4.507151074735273e-10
CG iteration 96, residual = 3.251276414432282e-10
CG iteration 97, residual = 2.674700036219682e-10
CG iteration 98, residual = 2.5564028417710706e-10
CG iteration 99, residual = 2.3795491139126934e-10
CG iteration 100, residual = 2.0591929252014045e-10
CG iteration 101, residual = 1.6036235947545246e-10
CG iteration 102, residual = 1.2318690339505628e-10
CG iteration 103, residual = 1.010960603013828e-10
CG iteration 104, residual = 9.240005362137239e-11
CG iteration 105, residual = 8.836556279209545e-11
CG iteration 106, residual = 7.445129267060336e-11
CG iteration 107, residual = 5.7113818976168324e-11
CG iteration 108, residual = 4.1527963128373796e-11
CG iteration 109, residual = 3.4821212507518504e-11
CG iteration 110, residual = 3.0543468656049e-11
CG iteration 111, residual = 2.7745482615326168e-11
CG iteration 112, residual = 2.498109459544577e-11
CG iteration 113, residual = 1.9814458146642636e-11
CG iteration 114, residual = 1.5338311006705303e-11
CG iteration 115, residual = 1.2338556475367854e-11
CG iteration 116, residual = 1.0409825400434616e-11
CG iteration 117, residual = 9.949464056624519e-12
CG iteration 118, residual = 8.925398346163866e-12
CG iteration 119, residual = 7.191566320001598e-12
CG iteration 120, residual = 6.114138905544884e-12
CG iteration 121, residual = 5.000595058424829e-12
CG iteration 122, residual = 4.050949882776923e-12
CG iteration 123, residual = 3.592894127502257e-12
CG iteration 124, residual = 3.15693023737965e-12
CG iteration 125, residual = 2.53737750012793e-12
CG iteration 126, residual = 2.1358046429079365e-12
CG iteration 127, residual = 1.7366171992889772e-12
CG iteration 128, residual = 1.436173546785084e-12
CG iteration 129, residual = 1.2028856050504542e-12
CG iteration 130, residual = 1.0452116697780117e-12
CG iteration 131, residual = 9.291844203513888e-13
CG iteration 132, residual = 8.162568277638706e-13
CG iteration 133, residual = 6.785851180002148e-13
CG iteration 134, residual = 5.382062590437613e-13
CG iteration 135, residual = 4.181747786049532e-13
CG iteration 136, residual = 3.606191990614787e-13
CG iteration 137, residual = 3.370837050888976e-13
CG iteration 138, residual = 3.087711118808599e-13
CG iteration 139, residual = 2.766120786902507e-13
CG iteration 140, residual = 2.1759148818729274e-13

[12]:

BaseWebGuiScene

[ ]: