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 \(V_0 \subset V_1 \subset \ldots V_L\). For given \(u \in V_L\), we want to compute a decomposition

\[ u = u_0 + \sum_{l=1}^L w_l \]

with \(u_0 \in V_0\) and \(w_l \in V_l\) such that

\[ \| u_0 \|_A^2 + \sum h_l^{-2} \| w_l \|_{L_2}^2 \preceq \| u \|_{H^1}^2 \]

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

Let \(N_l = \operatorname{dim}{V_l}\), and \(\varphi_{l,i}\) be the hat-basis-functions of \(V_l\). Nodal interpolation of \(u\) onto level \(l\) is not stable, but local averaging with the hat-basis is a good choice. We define

\[ u_l := \Pi_l u := \sum_{i=1}^{N_l} \frac{ (u, \varphi_{l,i})_{L_2} } { (1, \varphi_{l,i} )_{L_2} } \varphi_{l,i} \]

If \(u\) is a constant function, than \(\Pi_l u\) reproduces the same constant function.

For \(u \in V_L\) we define

\[\begin{eqnarray*} u_0 & := & \Pi_0 u \\ w_l & = & \Pi_l u - \Pi_{l-1} u \qquad 1 \leq l \leq L-1 \\ w_L & = & u - \Pi_{L-1} u, \end{eqnarray*}\]

which is a decomposition \(u = u_0 + \sum w_l\).

Theorem There holds

\[ \| u_0 \|_{H^1}^2 + \sum_{l=1}^L h_l^{-2} \| w_l \|_{L_2}^2 \preceq L^2 \, \| u \|_{H^1}^2 \]

Proof: We use that there exists a stable decomposition \(u = \tilde u_0 + \sum_{k=1}^L \tilde w_k\). The computable choice is

\[ w_l = \Pi_l u - \Pi_{l-1} u = (\Pi_l - \Pi_{l-1}) \tilde u_0 + \sum_{k=1}^L (\Pi_l - \Pi_{l-1}) \tilde w_k \]

For \(k \leq l\) there holds

\[ \| (\Pi_l - \Pi_{l-1}) \tilde w_k \|_{L_2} \preceq h_l \, \| \tilde w_k \|_{H^1} \preceq\frac{h_l}{h_k} \| \tilde w_k \|_{L_2}, \]

for \(k \geq l\) there holds

\[ \| (\Pi_l - \Pi_{l-1}) \tilde w_k \|_{L_2} \preceq \| \tilde w_k \|_{L_2} \preceq \frac{h_l}{h_k} \| \tilde w_k \|_{L_2}. \]

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

\[\begin{eqnarray*} \sum_{l=1}^L h_l^{-2} \| w_l \|_{L_2}^2 & \leq & \sum_{l=1}^L h_l^{-2} \Big( \sum_{k=1}^L \| (\Pi_l - \Pi_{l-1}) \tilde w_k \|_{L_2} \Big)^2 \\ & \leq & \sum_{l=1}^L h_l^{-2} \Big( \sum_{k=1}^L \frac{h_l}{h_k} \| \tilde w_k \|_{L_2} \Big)^2 \\ & = & \sum_{l=1}^L \Big( \sum_{k=1}^L \frac{1}{h_k} \| \tilde w_k \|_{L_2} \Big)^2 \\ & = & L \Big( \sum_{k=1}^L 1 \cdot \frac{1}{h_k} \| \tilde w_k \|_{L_2} \Big)^2 \\ & \leq & L \sum_{k=1}^L 1^2 \, \sum_{k=1}^L \frac{1}{h_k^2} \| \tilde w_k \|_{L_2}^2 \\ & = & L^2 \sum_{k=1}^L \frac{1}{h_k^2} \| \tilde w_k \|_{L_2}^2 \\ & \preceq & L^2 \, \| u \|_{H^1}^2 \end{eqnarray*}\]

86.2. Algorithm

The quasi-interpolant of \(u \in V_L\) to \(u_l \in V_l\),

\[ \Pi_l u := \sum_{i=1}^{N_l} \frac{ (u, \varphi_{l,i})_{L_2} } { (1, \varphi_{l,i} )_{L_2} } \varphi_{l,i} \]

can be implemented as follows. Let \(M_l\) be the mass-matrix on level \(l\), and \(\tilde P_l = P_L P_{L-1} \cdots P_{l+1}\) be the cumulated prolongation operators from level \(l\) up to the final level \(L\). Then

• \( (u, \varphi_{l,i})_{L_2} = (\tilde P_l^T M_L u)_i \)

• The values \((1, \varphi_{l,i})_{L_2}\) are the entries of the vector \(M_l {\mathbb 1}\), where \({\mathbb 1} = (1,\ldots, 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

\[ E_l^0 : H^1(\partial \Omega) \rightarrow H^1(\Omega) \]

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:

\[ E_l = E_{l-1} \Pi_{l-1} u + E_l^0 (u - \Pi_{l-1} u) \]

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