Hello All! I would like to use the operator of the L2 projection onto a coarser level but applied to a fine-level, trial function, not just to a grid function. So I employ a mixed formulation using prolongation and restriction matrices.

In order to test my understanding of the code I check if the fine level mass matrix M_h is exactly the restriction matrix R times the coarse level mass matrix M_H times the prolongation matrix P, or M_h=P*M_H*R. I construct a fine-level vector u_h and compare M_h*u_h and P*M_H*R u_h . The difference between two is not zero, though. Am I missing something? Many thanks!

from netgen.csg import CSGeometry, OrthoBrick, Pnt

from ngsolve import *

# Geometry

cube = CSGeometry()

cube.Add(OrthoBrick(Pnt(-2.0, -2.0, -2.0), Pnt(2.0, 2.0, 2.0)))

#Mesh generation

N=4

maxh = 4.0/N #to be refined later

maxH = 4.0/N

mesh = Mesh(cube.GenerateMesh(maxh=maxh, quad_dominated=False))

mesh_c = Mesh(cube.GenerateMesh(maxh=maxH, quad_dominated=False))

#FE Spaces

Vh = H1(mesh, order=1) #fine level space

u, v = Vh.TnT()

Rh = H1(mesh_c, order=1)#coarse level space

z, r = Rh.TnT()

Rh_avatar= H1(mesh, order=1)#used only to compute prolongation operator from coarse to fine

# Bilinear forms:

mass_c = BilinearForm(Rh, symmetric=True)#coarse level mass matrix M_H

mass_c += (u * v) * dx

mass_c.Assemble()

#Mesh refinement from level=0 to level=1

mesh.ngmesh.Refine()

Vh.Update()

Rh_avatar.Update() #Rh remains to be coarse of level=0

# Bilinear forms:

mass = BilinearForm(Vh, symmetric=True) #fine level mass matrix M_h

mass += (u * v) * dx

mass.Assemble()

#Grid function

coef=CoefficientFunction(((3 * pi * pi + 1)*sin(pi * x)*sin(pi * y)*sin(pi * z)))#an example function

gfu=GridFunction(Vh)

gfu.Set(coef);

# Prolongation operator from mesh-level 0 to level 1:

prol = Rh_avatar.Prolongation().Operator(1)

fine_to_coarse_to_fine_mass = prol @ mass_c.mat @ prol.T # with mass.mat here diff=0 trivially

#Check if M_h=Restriction * M_H * Prolongation holds

diff = gfu.vec.CreateVector()

diff += mass.mat*gfu.vec

diff -= fine_to_coarse_to_fine_mass * gfu.vec

print("Diff: ", sqrt(InnerProduct(diff, diff))) #non-zero!

exit(0)