Multi-level Extension

86. Multi-level Extension#

The goal is to extend boundary data (e.g. Dirichlet values) onto the domain, such that the resulting energy norm is optimal (up to a constant factor). We do so by a multi-level decomposition of the boundary data.

86.1. Efficiently computable multi-level decomposition#

We are given a nested sequence of spaces V0⊂V1⊂…VL. For given u∈VL, we want to compute a decomposition

u=u0+∑l=1Lwl

with u0∈V0 and wl∈Vl such that

‖u0‖A2+∑hl−2‖wl‖L22⪯‖u‖H12

We know that there exist such decompositions, but do not yet have a cheap algorithm for computing one. This we define next.

Let Nl=dimVl, and φl,i be the hat-basis-functions of Vl. Nodal interpolation of u onto level l is not stable, but local averaging with the hat-basis is a good choice. We define

ul:=Πlu:=∑i=1Nl(u,φl,i)L2(1,φl,i)L2φl,i

If u is a constant function, than Πlu reproduces the same constant function.

For u∈VL we define

u0:=Π0uwl=Πlu−Πl−1u1≤l≤L−1wL=u−ΠL−1u,

which is a decomposition u=u0+∑wl.

Theorem There holds

‖u0‖H12+∑l=1Lhl−2‖wl‖L22⪯L2‖u‖H12

Proof: We use that there exists a stable decomposition u=u~0+∑k=1Lw~k. The computable choice is

wl=Πlu−Πl−1u=(Πl−Πl−1)u~0+∑k=1L(Πl−Πl−1)w~k

For k≤l there holds

‖(Πl−Πl−1)w~k‖L2⪯hl‖w~k‖H1⪯hlhk‖w~k‖L2,

for k≥l there holds

‖(Πl−Πl−1)w~k‖L2⪯‖w~k‖L2⪯hlhk‖w~k‖L2.

Using the triangle inequality, the proven estimates and Cauchy-Schwarz:

∑l=1Lhl−2‖wl‖L22≤∑l=1Lhl−2(∑k=1L‖(Πl−Πl−1)w~k‖L2)2≤∑l=1Lhl−2(∑k=1Lhlhk‖w~k‖L2)2=∑l=1L(∑k=1L1hk‖w~k‖L2)2=L(∑k=1L1⋅1hk‖w~k‖L2)2≤L∑k=1L12∑k=1L1hk2‖w~k‖L22=L2∑k=1L1hk2‖w~k‖L22⪯L2‖u‖H12

86.2. Algorithm#

The quasi-interpolant of u∈VL to ul∈Vl,

Πlu:=∑i=1Nl(u,φl,i)L2(1,φl,i)L2φl,i

can be implemented as follows. Let Ml be the mass-matrix on level l, and P~l=PLPL−1⋯Pl+1 be the cumulated prolongation operators from level l up to the final level L. Then

  • (u,φl,i)L2=(P~lTMLu)i

  • The values (1,φl,i)L2 are the entries of the vector Ml1, where 1=(1,…,1).

86.3. Extending boundary data#

We use this coarse-level interpolants to extend Dirichlet data onto the domain. On every level we have a trivial extension

El0:H1(∂Ω)→H1(Ω)

by just setting all interior nodal values to zero. Thus, the function decays very quickly from the boundary. If the mesh gets more and more refined, the decay gets steeper and steeper.

from ngsolve import *
from ngsolve.webgui import Draw

mesh = Mesh(unit_square.GenerateMesh(maxh=0.3))
fes = H1(mesh, order=1, dirichlet="left|bottom", autoupdate=True)
u,v = fes.TnT()
gfu = GridFunction(fes, autoupdate=True)

for i in range(4):
    mesh.Refine()

a = BilinearForm(grad(u)*grad(v)*dx).Assemble()

copy mesh complete

bnd = mesh.Boundaries("left|bottom")
gfu.Set (1-x-y+0.3*sin(10*pi*x), definedon=bnd)

Draw (gfu, deformation=True)
print ("Norm(u) = ", InnerProduct((a.mat*gfu.vec).Evaluate(), gfu.vec))

Norm(u) =  38.983324464566756

The method is now to decompose the function into different levels, and use the trivial extension on every level to obtain a better extension:

El=El−1Πl−1u+El0(u−Πl−1u)
class MLExtension:
    def __init__ (self, fes, level, bndmass, dofs):
        self.fes = fes
        self.level = level
        self.bndmass = bndmass
        self.dofs = dofs
        
        ones = bndmass.CreateRowVector()
        ones[:] = 1
        Mones = (bndmass*ones).Evaluate()
        self.inv = DiagonalMatrix(Mones).Inverse(dofs)
        
        if level > 0:
            self.prol = fes.Prolongation().CreateMatrix(level)
            self.rest = self.prol.CreateTranspose()
            coarsebndmass = self.rest @ bndmass @ self.prol # multiply matrices
            coarsedofs = BitArray(self.prol.width)
            coarsedofs[:] = False
            for i in range(len(coarsedofs)):
                coarsedofs[i] = dofs[i]
            self.coarseext = MLExtension(fes, level-1, coarsebndmass, coarsedofs)
        
    def Extend (self, x):
        Mx = (self.bndmass * x).Evaluate()
        ext = self.ExtendRec(Mx)
        x[~self.dofs] = ext
        
    def ExtendRec (self, Mx):
        sol = (self.inv * Mx).Evaluate()
        if self.level == 0:
            return sol
        
        if self.level > 0:
            self.fes.Prolongation().Restrict(self.level, Mx)
            xc = self.coarseext.ExtendRec(Mx)
            pxc = (self.fes.Prolongation().Operator(self.level) * xc).Evaluate()
            
        pxc[self.dofs] = sol
        return pxc

bndmass = BilinearForm(u*v*ds(bnd), check_unused=False).Assemble().mat
ext = MLExtension(fes, fes.mesh.levels-1, bndmass, fes.GetDofs(bnd))
ext.Extend(gfu.vec)

Draw (gfu, deformation=True)
print ("Norm(uext) = ", InnerProduct(a.mat*gfu.vec, gfu.vec))

Norm(uext) =  3.541248535042474
Contents