Prolongation and Restriction operators

Main Forum Old-Forum - do not post here
1

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=PM_HR. I construct a fine-level vector u_h and compare M_hu_h and PM_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)

2

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.matgfz.vec
diff += (-1) coarse_to_fine_to_coarse_mass * gfz.vec
print("Diff: ", sqrt(InnerProduct(diff, diff))) #zero!
exit(0)