Prolongation and Restriction operators

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)

So this is what was my mistake. It is not M_h=P M_H R which should be true; only M_H=R M_h P holds. I post the corrected code that checks the last statement.

from netgen.csg import CSGeometry, OrthoBrick, Pnt
from ngsolve import *
from scipy.integrate import Radau

# Geometry
cube = CSGeometry()
cube.Add(OrthoBrick(Pnt(-2.0, -2.0, -2.0), Pnt(2.0, 2.0, 2.0)))

# Mesh sizes
N=4
maxh = 4.0/N #to be refined later
maxH = 4.0/N

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

#forcing
pi=3.14
coef=CoefficientFunction(((3 * pi * pi + 1)*sin(pi * x)*sin(pi * y)*sin(pi * z)))

#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

#Grid function
gfz=GridFunction(Rh)
gfz.Set(coef);

# Prolongation operator from mesh-level 0 to level 1:
prol = Rh_avatar.Prolongation().Operator(1)
coarse_to_fine_to_coarse_mass   = prol.T @ mass.mat @ prol

#Check if M_H=Restriction * M_h * Prolongation holds
diff  = gfz.vec.CreateVector()
diff += mass_c.mat*gfz.vec
diff += (-1)* coarse_to_fine_to_coarse_mass * gfz.vec
print("Diff: ", sqrt(InnerProduct(diff, diff))) #zero!
exit(0)