diff --git a/docs/src/PETScPC/poisson.py b/docs/src/PETScPC/poisson.py
deleted file mode 100644
index 2df07d8..0000000
--- a/docs/src/PETScPC/poisson.py
+++ /dev/null
@@ -1,158 +0,0 @@
-# Vertex Patch smoothing for p-Multigrid preconditioners for the Poisson problem
-# ===============================================================================
-#
-# In this tutorial, we explore using `PETSc PC` as a preconditioner inside NGSolve preconditioning infrastructure.
-# We will focus our attention on the Poisson problem, and we will consider different preconditioning strategies.
-# Not all the preconditioning strategies are equally effective for all the discretizations of the Poisson problem, and this demo is intended to provide a starting point for the exploration of the preconditioning strategies rather than providing a definitive answer.
-# We begin by creating a discretisation of the Poisson problem using H1 elements, in particular, we consider the usual variational formulation
-#
-# .. math::
-#
-#    \text{find } u\in H^1_0(\Omega) \text{ s.t. } a(u,v) := \int_{\Omega} \nabla u\cdot \nabla v \; d\vec{x} = L(v) := \int_{\Omega} fv\; d\vec{x}\qquad v\in H^1_0(\Omega).
-#
-# Such a discretisation can easily be constructed using NGSolve as follows: ::
-
-from ngsolve import *
-import netgen.gui
-import netgen.meshing as ngm
-from mpi4py.MPI import COMM_WORLD
-
-if COMM_WORLD.rank == 0:
-   ngmesh = unit_square.GenerateMesh(maxh=0.1)
-   for _ in range(4):
-      ngmesh.Refine()
-   mesh = Mesh(ngmesh.Distribute(COMM_WORLD))
-else:
-   mesh = Mesh(ngm.Mesh.Receive(COMM_WORLD))
-order = 2
-fes = H1(mesh, order=order, dirichlet="left|right|top|bottom")
-print("Number of degrees of freedom: ", fes.ndof)
-u,v = fes.TnT()
-a = BilinearForm(grad(u)*grad(v)*dx)
-f = LinearForm(fes)
-f += 32 * (y*(1-y)+x*(1-x)) * v * dx
-a.Assemble()
-f.Assemble()
-
-# We now construct an NGSolve preconditioner wrapping a `PETSc PC`, we begin considering an Algebraic MultiGrid preconditioner constructed using `HYPRE`.
-# In this tutorial, we will use the Krylov solver implemented inside NGSolve to solve the linear system. ::
-
-from ngsPETSc.pc import *
-from ngsolve.krylovspace import CG
-pre = Preconditioner(a, "PETScPC", pc_type="hypre")
-gfu = GridFunction(fes)
-print("-------------------|HYPRE p={}|-------------------".format(order))
-gfu.vec.data = CG(a.mat, rhs=f.vec, pre=pre.mat, printrates=True)
-Draw(gfu)
-
-# We see that the HYPRE preconditioner is quite effective for the Poisson problem discretised using linear elements, but it is not as effective for higher-order elements.
-#
-# .. list-table:: Preconditioners performance HYPRE
-#    :widths: auto
-#    :header-rows: 1
-#
-#    * - Preconditioner
-#      - p=1
-#      - p=2
-#      - p=3
-#      - p=4
-#    * - HYPRE
-#      - 10 (1.57e-12)
-#      - 70 (2.22e-13)
-#      - 42 (3.96e-13)
-#      - 70 (1.96e-13)
-#
-# To overcome this issue we will use a two-level additive Schwarz preconditioner.
-# In this case, we will use as fine space correction, the inverse of the local matrices associated with the patch of a vertex. ::
-
-def VertexPatchBlocks(mesh, fes):
- blocks = []
- freedofs = fes.FreeDofs()
- for v in mesh.vertices:
-     vdofs = set()
-     for el in mesh[v].elements:
-         vdofs |= set(d for d in fes.GetDofNrs(el)
-                      if freedofs[d])
-     blocks.append(vdofs)
- return blocks
-
-blocks = VertexPatchBlocks(mesh, fes)
-dofs = BitArray(fes.ndof); dofs[:] = True
-blockjac = ASMPreconditioner(a.mat, dofs, blocks=blocks,
-                               solverParameters={"pc_type": "asm",
-                                                 "sub_ksp_type": "preonly",
-                                                 "sub_pc_type": "lu"})  
-
-# We now isolate the degrees of freedom associated with the vertices and construct a two-level additive Schwarz preconditioner, where the coarse space correction is the inverse of the local matrices associated with the vertices. ::
-
-def VertexDofs(mesh, fes):
-   vertexdofs = BitArray(fes.ndof)
-   vertexdofs[:] = False
-   for v in mesh.vertices:
-      for d in fes.GetDofNrs(v):
-         vertexdofs[d] = True
-   vertexdofs &= fes.FreeDofs()
-   return vertexdofs
-
-vertexdofs = VertexDofs(mesh, fes)
-preCoarse = PETScPreconditioner(a.mat, vertexdofs, solverParameters={"pc_type": "hypre"})
-pretwo = preCoarse + blockjac
-print("-------------------|Additive Schwarz p={}|-------------------".format(order))
-gfu.vec.data = CG(a.mat, rhs=f.vec, pre=pretwo, printrates=True)
-
-# We can see that the two-level additive Schwarz preconditioner where the coarse space correction is performed using HYPRE is more effective than using just HYPRE preconditioner for higher-order elements.
-#
-# .. list-table:: Preconditioners performance HYPRE and Two Level Additive Schwarz
-#    :widths: auto
-#    :header-rows: 1
-#
-#    * - Preconditioner
-#      - p=1
-#      - p=2
-#      - p=3
-#      - p=4
-#    * - HYPRE
-#      - 10 (1.57e-12)
-#      - 70 (2.22e-13)
-#      - 42 (3.96e-13)
-#      - 70 (1.96e-13)
-#    * - Additive Schwarz
-#      - 44 (1.96e-12)
-#      - 45 (1.28e-12)
-#      - 45 (1.29e-12)
-#      - 45 (1.45e-12)
-#       
-#   
-# We can also use the PETSc preconditioner as an auxiliary space preconditioner.
-# Let us consider the discontinuous Galerkin discretisation of the Poisson problem. ::
-
-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*3**2/h*jump_u*jump_v * dx(skeleton=True)
-aDG += alpha*3**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)
-fDG = LinearForm(fesDG)
-fDG += 1*v * dx
-aDG.Assemble()
-fDG.Assemble()
-
-# We can now use the PETSc PC assembled for the conforming Poisson problem as an auxiliary space preconditioner for the DG discretisation. ::
-
-from ngsPETSc import pc
-smoother = Preconditioner(aDG, "PETScPC", pc_type="jacobi")
-transform = fes.ConvertL2Operator(fesDG)
-preDG = transform @ pre.mat @ transform.T + smoother.mat
-gfuDG = GridFunction(fesDG)
-print("-------------------|Auxiliary Space preconditioner p={}|-------------------".format(order))
-gfuDG.vec.data = CG(aDG.mat, rhs=fDG.vec, pre=preDG, printrates=True)
-Draw(gfuDG)
\ No newline at end of file
diff --git a/docs/src/PETScPC/stokes.py b/docs/src/PETScPC/stokes.py
deleted file mode 100644
index 56964a9..0000000
--- a/docs/src/PETScPC/stokes.py
+++ /dev/null
@@ -1,347 +0,0 @@
-# Preconditioning the Stokes problem 
-# ===================================
-#
-# In this tutorial, we explore constructing preconditioners for saddle point problems using `PETSc PC`.
-# In particular, we will consider a Bernardi-Raugel inf-sup stable discretization of the Stokes problem, i.e.
-#
-# .. math::       
-#   
-#    \text{find } (\vec{u},p) \in [H^1_{0}(\Omega)]^d\times L^2(\Omega) \text{ s.t. }
-#   
-#    \begin{cases} 
-#       (\nabla \vec{u},\nabla \vec{v})_{L^2(\Omega)} + (\nabla\cdot \vec{v}, p)_{L^2(\Omega)}  = (\vec{f},\vec{v})_{L^2(\Omega)} \qquad \forany v\in H^1_{0}(\Omega)\\
-#       (\nabla\cdot \vec{u},q)_{L^2(\Omega)} = 0 \qquad q\in L^2(\Omega)
-#    \end{cases}
-#
-# Such a discretization can easily be constructed using NGSolve as follows: ::
-
-from ngsolve import *
-from ngsolve import BilinearForm as BF
-from netgen.occ import *
-import netgen.gui
-import netgen.meshing as ngm
-from mpi4py.MPI import COMM_WORLD
-
-if COMM_WORLD.rank == 0:
-   shape = Rectangle(2,0.41).Circle(0.2,0.2,0.05).Reverse().Face()
-   shape.edges.name="wall"
-   shape.edges.Min(X).name="inlet"
-   shape.edges.Max(X).name="outlet"
-   geo = OCCGeometry(shape, dim=2)
-   ngmesh = geo.GenerateMesh(maxh=0.1)
-   mesh = Mesh(ngmesh.Distribute(COMM_WORLD))
-else:
-   mesh = Mesh(ngm.Mesh.Receive(COMM_WORLD))
-nu = Parameter(1.0)
-V = VectorH1(mesh, order=4, dirichlet="wall|inlet|cyl", autoupdate=True)
-Q = L2(mesh, order=3, autoupdate=True)
-u,v = V.TnT(); p,q = Q.TnT()
-a = BilinearForm(nu*InnerProduct(Grad(u),Grad(v))*dx)
-a.Assemble()
-b = BilinearForm(div(u)*q*dx)
-b.Assemble()
-gfu = GridFunction(V, name="u")
-gfp = GridFunction(Q, name="p")
-uin = CoefficientFunction( (1.5*4*y*(0.41-y)/(0.41*0.41), 0) )
-gfu.Set(uin, definedon=mesh.Boundaries("inlet"))
-f = LinearForm(V).Assemble()
-g = LinearForm(Q).Assemble();
-
-# We can now explore what happens if we use a Schour complement preconditioner for the saddle point problem.
-# We can construct the Schur complement preconditioner using the following code: ::
-
-K = BlockMatrix( [ [a.mat, b.mat.T], [b.mat, None] ] )
-from ngsPETSc.pc import *
-apre = Preconditioner(a, "PETScPC", pc_type="lu")
-S = (b.mat @ apre.mat @ b.mat.T).ToDense().NumPy()
-from numpy.linalg import pinv
-from scipy.sparse import coo_matrix
-from ngsolve.la import SparseMatrixd 
-Sinv = coo_matrix(pinv(S))
-Sinv = la.SparseMatrixd.CreateFromCOO(indi=Sinv.row, 
-                                      indj=Sinv.col,
-                                      values=Sinv.data,
-                                      h=S.shape[0],
-                                      w=S.shape[1])
-C = BlockMatrix( [ [a.mat.Inverse(V.FreeDofs()), None],
-                   [None, Sinv] ] )
-
-rhs = BlockVector (  [f.vec, g.vec] )
-sol = BlockVector( [gfu.vec, gfp.vec] )
-print("-----------|Schur|-----------")
-solvers.MinRes (mat=K, pre=C, rhs=rhs, sol=sol, tol=1e-8,
-                printrates=True, initialize=False)
-Draw(gfu)
-
-# The Schur complement preconditioner converges in a few iterations, but it is not very efficient since we need to invert the Schur complement.
-#
-# .. list-table:: Preconditioners performance
-#    :widths: auto
-#    :header-rows: 1
-#
-#    * - Preconditioner
-#      - Iterations
-#    * - Schur complement
-#      - 4 (9.15e-14)
-#
-# Notice that the Schur complement is dense hence inverting it is not a good idea. Not only that but to perform the inversion of the Schur complement had to write a lot of "boilerplate" code.
-# Since our discretization is inf-sup stable it is possible to prove that the mass matrix of the pressure space is spectrally equivalent to the Schur complement.
-# This means that we can use the mass matrix of the pressure space as a preconditioner for the Schur complement.
-# Notice that we still need to invert the mass matrix and we will do so using a `PETSc PC` of type Jacobi, which is the exact inverse since we are using `P0` elements.
-# We will also invert the Laplacian block using a `PETSc PC` of type `LU`. ::
-
-m = BilinearForm((1/nu)*p*q*dx).Assemble()
-mpre = Preconditioner(m, "PETScPC", pc_type="lu")
-apre = a.mat.Inverse(V.FreeDofs())
-C = BlockMatrix( [ [apre, None], [None, mpre] ] )
-
-gfu.vec.data[:] = 0; gfp.vec.data[:] = 0;
-gfu.Set(uin, definedon=mesh.Boundaries("inlet"))
-sol = BlockVector( [gfu.vec, gfp.vec] )
-
-print("-----------|Mass & LU|-----------")
-solvers.MinRes (mat=K, pre=C, rhs=rhs, sol=sol, tol=1e-8,
-                maxsteps=100, printrates=True, initialize=False)
-Draw(gfu)
-
-# .. list-table:: Preconditioners performance
-#    :widths: auto
-#    :header-rows: 1
-#
-#    * - Preconditioner
-#      - Iterations
-#    * - Schur complement
-#      - 4 (9.15e-14)
-#    * - Mass & LU
-#      - 66 (2.45e-08)
-#   
-# We can also construct a multi-grid preconditioner for the top left block of the saddle point problem, as we have seen in :doc:`poisson.py`. ::
-
-def DoFInfo(mesh, fes):
-   blocks = []
-   freedofs = fes.FreeDofs()
-   vertexdofs = BitArray(fes.ndof)
-   vertexdofs[:] = False
-   for v in mesh.vertices:
-      vdofs = set()
-      vdofs |= set(d for d in fes.GetDofNrs(v) if freedofs[d])
-      for ed in mesh[v].edges:
-         vdofs |= set(d for d in fes.GetDofNrs(ed) if freedofs[d])
-      for fc in mesh[v].faces:
-         vdofs |= set(d for d in fes.GetDofNrs(fc) if freedofs[d])
-      blocks.append(vdofs)
-      for d in fes.GetDofNrs(v):
-         vertexdofs[d] = True
-   vertexdofs &= fes.FreeDofs()
-   return vertexdofs, blocks 
-
-vertexdofs, blocks = DoFInfo(mesh, V)
-blockjac = a.mat.CreateBlockSmoother(blocks)
-preH = PETScPreconditioner(a.mat, vertexdofs, solverParameters={"pc_type":"hypre"})
-twolvpre = preH + blockjac
-C = BlockMatrix( [ [twolvpre, None], [None, mpre] ] )
-gfu.vec.data[:] = 0; gfp.vec.data[:] = 0;
-gfu.Set(uin, definedon=mesh.Boundaries("inlet"))
-print("-----------|Mass & Two Level Additive Schwarz|-----------")
-solvers.MinRes (mat=K, pre=C, rhs=rhs, sol=sol, tol=1e-8,
-                maxsteps=100, printrates=True, initialize=False)
- 
-# .. list-table:: Preconditioners performance
-#    :widths: auto
-#    :header-rows: 1
-#
-#    * - Preconditioner
-#      - Iterations
-#    * - Schur complement
-#      - 4 (9.15e-14)
-#    * - Mass & LU
-#      - 66 (2.45e-08)
-#    * - Mass & Two Level Additive Schwarz
-#      - 100 (4.68e-06)
-#   
-# The mass matrix as a preconditioner doesn't seem to be ideal, in fact, our Krylov solver took many iterations to converge with a direct LU factorization of the velocity block and did not converge at all with `HYPRE`.
-# To resolve this issue we resort to an augmented Lagrangian formulation, i.e.
-#
-# .. math::
-#    \begin{cases} 
-#       (\nabla \vec{u},\nabla \vec{v})_{L^2(\Omega)} + (\nabla\cdot \vec{v}, p)_{L^2(\Omega)} + \gamma (\nabla\cdot \vec{u},\nabla\cdot\vec{v})_{L^2(\Omega)} = (\vec{f},\vec{v})_{L^2(\Omega)} \qquad \forany v\in H^1_{0}(\Omega)\\
-#       (\nabla\cdot \vec{u},q)_{L^2(\Omega)} = 0 \qquad q\in L^2(\Omega)
-#    \end{cases}
-#
-# This formulation can easily be constructed by adding a new velocity block in the `BlockMatrix`, as follows: ::
-
-gamma = Parameter(1e8)
-aG = BilinearForm(nu*InnerProduct(Grad(u),Grad(v))*dx+gamma*div(u)*div(v)*dx)
-aG.Assemble()
-aGpre = Preconditioner(aG, "PETScPC", pc_type="lu")
-mG = BilinearForm((1/nu+gamma)*p*q*dx).Assemble()
-mGpre = Preconditioner(mG, "PETScPC", pc_type="jacobi")
-
-K = BlockMatrix( [ [aG.mat, b.mat.T], [b.mat, None] ] )
-C = BlockMatrix( [ [aGpre.mat, None], [None, mGpre.mat] ] )
-
-gfu.vec.data[:] = 0; gfp.vec.data[:] = 0;
-gfu.Set(uin, definedon=mesh.Boundaries("inlet"))
-sol = BlockVector( [gfu.vec, gfp.vec] )
-
-print("-----------|Augmented LU|-----------")
-solvers.MinRes (mat=K, pre=C, rhs=rhs, sol=sol, tol=1e-10,
-                printrates=True, initialize=False)
-Draw(gfu)
-
-# Using an augmented Lagrangian formulation, we were able to converge in only two iterations.
-# This is because the augmented Lagrangian formulation improves the spectral equivalence of the mass matrix of the pressure space and the Schur complement.
-# 
-# .. list-table:: Preconditioners performance
-#    :widths: auto
-#    :header-rows: 1
-#
-#    * - Preconditioner
-#      - Iterations
-#    * - Schur complement
-#      - 4 (9.15e-14)
-#    * - Mass & LU
-#      - 66 (2.45e-08)
-#    * - Mass & Two Level Additive Schwarz
-#      - 100 (4.68e-06)
-#    * - Augmented Lagrangian LU
-#      - 2 (9.24e-8)
-#
-# Notice that so far we have been inverting the matrix corresponding to the Laplacian block using a direct LU factorization.
-# This is not ideal for large problems, and we can use a `Hypre` preconditioner for the Laplacian block. ::
-
-smoother = aG.mat.CreateBlockSmoother(blocks)
-preHG = PETScPreconditioner(aG.mat, vertexdofs, solverParameters={"pc_type":"hypre"})
-twolvpre = preHG + smoother
-C = BlockMatrix( [ [twolvpre, None], [None, mGpre] ] )
-gfu.vec.data[:] = 0; gfp.vec.data[:] = 0;
-gfu.Set(uin, definedon=mesh.Boundaries("inlet"))
-print("-----------|Augmented Two Level Additive Schwarz|-----------")
-solvers.MinRes (mat=K, pre=C, rhs=rhs, sol=sol, tol=1e-10,
-                printrates=True, initialize=False)
-Draw(gfu)
-
-# Our first attempt at using a `HYPRE` preconditioner for the Laplacian block did not converge.
-#
-# .. list-table:: Preconditioners performance
-#    :widths: auto
-#    :header-rows: 1
-#
-#    * - Preconditi¬oner
-#      - Iterations
-#    * - Schur complement
-#      - 4 (9.15e-14)
-#    * - Mass & LU
-#      - 66 (2.45e-08)
-#    * - Mass & Two Level Additive Schwarz
-#      - 100 (4.68e-06)
-#    * - Augmented Lagrangian LU
-#      - 2 (9.24e-08)
-#    * - Augmented Two Level Additive Schwarz
-#      - 100 (1.06e-03)
-#
-# This is because the top left block of the saddle point problem now contains the augmentation term, which has a very large kernel.
-# It is well known that algebraic multi-grid methods do not work well with indefinite problems, and this is what we are observing here.
-# We begin by constructing the augmented Lagrangian formulation in more numerical linear algebra terms, i.e. 
-#
-# .. math::
-#    \begin{bmatrix}
-#       A + B^T (\gamma M^{-1}) B & B^T \\
-#       B & 0
-#    \end{bmatrix}
-#    \begin{bmatrix}
-#       u \\
-#       p
-#    \end{bmatrix}
-#    =
-#    \begin{bmatrix}
-#       f \\
-#       0
-#    \end{bmatrix}
-#
-# We can construct this linear algebra problem inside NGSolve as follows ::
-
-d = BilinearForm((1/gamma)*p*q*dx)
-d.Assemble()
-dpre = PETScPreconditioner(d.mat, Q.FreeDofs(), solverParameters={"pc_type":"lu"})
-aG = a.mat + b.mat.T@dpre@b.mat
-aG = coo_matrix(aG.ToDense().NumPy())
-aG = la.SparseMatrixd.CreateFromCOO(indi=aG.row, 
-                                      indj=aG.col,
-                                      values=aG.data,
-                                      h=aG.shape[0],
-                                      w=aG.shape[1])
-K = BlockMatrix( [ [aG, b.mat.T], [b.mat, None] ] )
-pre = PETScPreconditioner(aG, V.FreeDofs(), solverParameters={"pc_type":"lu"})
-C = BlockMatrix( [ [pre, None], [None, mGpre.mat] ] )
-
-gfu.vec.data[:] = 0; gfp.vec.data[:] = 0;
-gfu.Set(uin, definedon=mesh.Boundaries("inlet"))
-sol = BlockVector( [gfu.vec, gfp.vec] )
-
-print("-----------|Boffi--Lovadina Augmentation LU|-----------")
-solvers.MinRes (mat=K, pre=C, rhs=rhs, sol=sol, tol=1e-10,
-                printrates=True, initialize=False)
-Draw(gfu)
-
-# We can now think of a more efficient way to invert the matrix corresponding to the augmentation term.
-# In fact, since we know that the augmentation block has a lower rank than the Laplacian block, we can use the Sherman-Morrisson-Woodbory formula to invert the augmentation block.
-#
-# .. math::
-#    (A + B^T(\gamma M^{-1})B)^{-1} = A^{-1} - A^{-1}B^T(\frac{1}{\gamma}M^{-1} + BA^{-1}B^T)^{-1}BA^{-1}
-#
-# We will do this in two different ways first we will invert the :math:`(\frac{1}{\gamma}M^{-1} + BA^{-1}B^T)` block using a direct LU factorization.
-# Then we will notice that since the penalisation parameter is large we can ignore the :math:`\frac{1}{\gamma}M^{-1}` term and use the mass matrix since it is spectrally to the Schur complement. ::
-
-SM = (d.mat + b.mat@apre@b.mat.T).ToDense().NumPy()
-SM = coo_matrix(SM)
-SM = la.SparseMatrixd.CreateFromCOO(indi=SM.row, 
-                                      indj=SM.col,
-                                      values=SM.data,
-                                      h=SM.shape[0],
-                                      w=SM.shape[1])
-
-SMinv = PETScPreconditioner(SM, Q.FreeDofs(), solverParameters={"pc_type":"lu"})
-
-C = BlockMatrix( [ [apre - apre@b.mat.T@SMinv@b.mat@apre, None], [None, mGpre.mat] ] )
-
-gfu.vec.data[:] = 0; gfp.vec.data[:] = 0;
-gfu.Set(uin, definedon=mesh.Boundaries("inlet"))
-sol = BlockVector( [gfu.vec, gfp.vec] )
-
-print("-----------|Boffi--Lovadina Augmentation Sherman-Morrisson-Woodbory|-----------")
-solvers.MinRes (mat=K, pre=C, rhs=rhs, sol=sol, tol=1e-10,
-                printrates=True, initialize=False)
-Draw(gfu)
-
-C = BlockMatrix( [ [apre + apre@(b.mat.T@mpre.mat@b.mat)@apre, None], [None, mGpre.mat] ] )
-
-gfu.vec.data[:] = 0; gfp.vec.data[:] = 0;
-gfu.Set(uin, definedon=mesh.Boundaries("inlet"))
-sol = BlockVector( [gfu.vec, gfp.vec] )
-
-print("-----------|Boffi--Lovadina Augmentation Sherman-Morrisson-Woodbory|-----------")
-solvers.MinRes (mat=K, pre=C, rhs=rhs, sol=sol, tol=1e-13,
-                printrates=True, initialize=False)
-Draw(gfu)
-
-# We see that a purely algebraic approach based on the Sherman-Morrisson-Woodbory formula is more efficient for the augmented Lagrangian formulation, then a naive two-level additive Schwarz approach.
-#
-# .. list-table:: Preconditioners performance
-#    :widths: auto
-#    :header-rows: 1
-#
-#    * - Preconditioner
-#      - Iterations
-#    * - Schur complement
-#      - 4 (9.15e-14)
-#    * - Mass & LU
-#      - 66 (2.45e-08)
-#    * - Mass & Two Level Additive Schwarz
-#      - 100 (4.68e-06)
-#    * - Augmented Lagrangian LU
-#      - 2 (9.24e-08)
-#    * - Augmented Two Level Additive Schwarz
-#      - 100 (1.06e-03)
-#    * - Augmentation Schermon-Morrisson-Woodbory
-#      - 84 (1.16e-07)