L2 projection from corse mesh FESpace to refined mesh FESpace

Wenzheng_Kuang · June 28, 2022, 4:04am

Dear developers,

Assume Q0 is a lowest order L2 FESpace on the current coarse mesh. After hierarchical refinement of the mesh, I construct a new L2 FESpace Q on the refined mesh, without updating the “old” Q0 on the coarser mesh. I am wondering if NGSolve supports calculating the projection from L2 space Q0 on the coarser mesh to L2 space Q on the finer mesh, or vice versa? I try to do this by using BilinearForm(trialspace=Q0, testspace=Q), but this throws an error as expected. Is there any way to do this?

mneunteufel · June 28, 2022, 7:03am

Hi,

it is possible to interpolate a GridFunction to another one living on a different mesh.

gf1 = GridFunction(L2(mesh1,order=0))
gf1.Set(…)

gf2 = GridFunction(L2(mesh2,order=0))
gf2.Set(gf1)

This works in both directions. The Set method does an L2-projection.

Best,
Michael

hvwahl · June 28, 2022, 7:18am

Hi,

you can also do the L2 projection of a FE function from one mesh to the other explicitly

from ngsolve import *
from ngsolve.meshes import MakeStructured2DMesh
from ngsolve.webgui import Draw

mesh1 = MakeStructured2DMesh(nx=10, ny=10)
Q0 = L2(mesh1, order=0)
gfu_q0 = GridFunction(Q0)
gfu_q0.Set(y)

Draw(gfu_q0, mesh1, ‘q0’)

mesh2 = MakeStructured2DMesh(nx=40, ny=40)
Q = L2(mesh2, order=0)

f = LinearForm(Q)
f += gfu_q0 * Q.TestFunction() * dx
f.Assemble()

gfu_q = GridFunction(Q)
gfu_q.vec.data = Q.InvM() * f.vec

Draw(gfu_q, mesh2, ‘q’)

Best wishes,
Henry

joachim · June 28, 2022, 12:51pm

Hi,

I guess you want to utilize the hierarchy of spaces, and not considering coarse and fine spaces as unrelated finite element spaces.

You can assemble matrices only on the current (fine) mesh. However, some spaces provide grid-transfer operators (=prolongation matrix).

An examples is the (P1-iso-P2)-P1 element pairing for Stokes:
https://docu.ngsolve.org/latest/i-tutorials/unit-2.6-stokes/stokes.html#(P1-iso-P2)-P1

or multigrid-algorithms in the i-FEM tutorials:

github.com

JSchoeberl/iFEM/blob/master/multigrid/algorithms.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "mediterranean-workplace",
   "metadata": {},
   "source": [
    "Multigrid and Multilevel Methods\n",
    "===\n",
    "\n",
    "Multigrid (MG) and Multilevel (ML) algorithms provide preconditioners with optimal condition numbers $\\kappa (C^{-1} A) = O(1)$, and optimal computational complexity $O(N)$.\n",
    "\n",
    "They can be seen as extension of the two level overlapping domain decomposition method to more levels.\n",
    "\n",
    "* Ulrich Trottenberg, Cornelius W. Oosterlee, Anton Schuller: Multigrid, Academic Press, 2001\n",
    "* Wolfgang Hackbusch: Multi-Grid Methods and Applications, Springer, 1985\n",
    "\n",
    "In short, both methods are sub-space correction methods where the space splitting is defined by all basis functions together from a hierarchy of grids. Multilevel methods are additive Schwarz methods, while Multigrid methods are the multiplicative Schwarz methods. Both can be implemented efficiently by recursive algorithms. \n",
    "\n",
    "We will present different theories for their analysis, one is based on sub-space decomposition using the ASM lemma, the other one uses the smoothing-and-approximation properties."

This file has been truncated. show original

Should work for the embedding of the coarse-L2(order=0) into a finer one.

Joachim

Wenzheng_Kuang · December 29, 2022, 3:47am

Hi Dr. Schoberl,

Thanks for your previous reply.
However, when I follow the link github.com/JSchoeberl/iFEM/blob/master/m…rid/algorithms.ipynb and just change the H1 space to L2 (order=0) space, there is error in “fes.Prolongtaion().LevelDofs(level-1)” saying “Illegal level 0 num levels = 0”. It seems that for now there is no LevelDofs() updating for the L2 space when the mesh is refined, right?

str · December 30, 2022, 8:53pm

Dear All,
I’ve been also curious about the L2 projection operator from the fine to the coarse levels - I guess , this is what is meant in the topic title). I’ve been testing an implementation of it a little (see other recent posts of mine) and found that MinRes, for example, throws an error when the prolongation matrix is computed between L2 spaces:
netgen.libngpy._meshing.NgException: Prolongation::GetNDofLevel not overloaded
This behaviour doesn’t manifest for prolongations between H1 spaces. I attach the file where the only thing that happens is the multiplication of a course vector with something that involves prolongation matrix. In fact, it should be the coarse mass matrix computed from the fine level mass matrix. Any help is appreciated!

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 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.141592
coef=CoefficientFunction(((3 * pi * pi + 1)*sin(pi * x)*sin(pi * y)*sin(pi * z)))

#FE Spaces
#changing all the L2 spaces to H1 resolves the bug
Vh = L2(mesh, order=1) #fine level space
u, v = Vh.TnT()

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

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

#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
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#=coarse mass matrix

#computing coarse mass_matrix * coarse vector
diff = gfz.vec.CreateVector()
diff += coarse_to_fine_to_coarse_mass * gfz.vec
#!!!error: netgen.libngpy._meshing.NgException: Prolongation::GetNDofLevel not overloaded
exit(0);