MG two-level restrict_and_add functionality

[dominiktassilostill]

While working on #16856 (especially on @peterrum's comment #16856 (comment)) I noticed, that restriction function MGTwoLevelTransfer::restrict_and_add_internal() seems to be faulty, not only for simplex elements but also hypercube elements.

When running the test multigrid-global-coarsening/mg-transfer_a_05 with a constant function the resulting restricted vector is non-constant. As the refinement of Tet elements was changed this was expected for these kinds of elements but it is also happening for hypercube elements in 3D with polynomial degrees of 3 and above (test<3, double>(3, Functions::ConstantFunction<3, double>(1.), false); in main()).

Changing in MGTwoLevelTransfer::restrict_and_add_internal() prolongation_matrix to restriction_matrix reduced the error for hypercube elements but it does not vanish.

CellRestrictor<dim, double, VectorizedArrayType>
                  cell_restrictor(scheme.restriction_matrix, // was prolongation_matrix before,
                                  scheme.restriction_matrix_1d, // was prolongation_matrix_1d before,
                                  evaluation_data_fine.begin() +
                                    c * n_scalar_dofs_fine,
                                  evaluation_data_coarse.begin() +
                                    c * n_scalar_dofs_coarse);

With @kronbichler we discussed to introduce new transfer.schemes in MGTwoLevelTransfer for the different refinement choices of the Tetrahedrons to fix the simplex transfer, yet I'm not sure where the error is for the hypercube elements.

[kronbichler]

When running the test multigrid-global-coarsening/mg-transfer_a_05 with a constant function the resulting restricted vector is non-constant.

Could you please elaborate on this a bit more? I would not expect the function to remain constant, because the restriction transfers residuals (i.e., quantities multiplied by test functions and integrated over the domain, or functionals in mathematical terms), so the question is what a constant input vector to the restriction would mean. In any case, I would not change the prolongation matrix without further insight, because we should implement the transpose of the prolongation (which one can interpret as applying to dual functionals).

[bangerth]

Along the same lines as @kronbichler , you typically choose the restriction operation as the transpose of the interpolation from coarse to fine grid. It is true that if you start with a constant function, interpolating it onto a finer grid should result in a function that is constant and has the same value. The important point is that the restriction is simply the transpose of this operation, and is -- despite its name -- not the interpolation from fine to coarse mesh. In other words, the multigrid "restriction" operation is not what the FiniteElement::get_restriction_matrix() function returns.

[bangerth]

See also #16906.

[dominiktassilostill]

Thank you both for the clarification. If the restriction operator is the transpose of the prolongation operator, it makes sense that the residual does not stay constant, opposite to injection with the matrix returned by FiniteElement::get_restriction_matrix(). As far as I can see the transposed prolongation matrix is implemented in the MGTwoLevelTransfer::restrict_and_add_internal() function.
The point that was unclear to me is that 'restriction' refers to different matrices/operations in the different contexts, your new note @bangerth is really of benefit to avoid confusion, as I just thought the restriction should be independent of what kind of vector gets restricted (e.g. solution or residual).