Tetrahedral refinement and FiniteElement::prolongation/restriction.

[bangerth]

As mentioned in #16755 (comment), now that we refine tetrahedra in different ways depending on the cell's geometry, a bunch of things are likely wrong. Specifically (copied from there):

I'm coming a bit late to this, but if the subdivision of a cell into children now depends on the specific cell, how does that affect what we do in FiniteElement::get_prolongation_matrix()? This function provides a matrix that facilitates the prolongation of values from the mother to a child cell. I don't know whether we implement these for tet cells already, but wouldn't we have to adapt what the function does?

Specifically, we currently have

template <int dim, int spacedim>
const FullMatrix<double> &
FE_SimplexPoly<dim, spacedim>::get_prolongation_matrix(
  const unsigned int         child,
  const RefinementCase<dim> &refinement_case) const
{
  Assert(refinement_case == RefinementCase<dim>::isotropic_refinement,
         ExcNotImplemented());
  AssertDimension(dim, spacedim);

  // initialization upon first request
  if (this->prolongation[refinement_case - 1][child].n() == 0)
    {
      std::lock_guard<std::mutex> lock(prolongation_matrix_mutex);

      // if matrix got updated while waiting for the lock
      if (this->prolongation[refinement_case - 1][child].n() ==
          this->n_dofs_per_cell())
        return this->prolongation[refinement_case - 1][child];

      // now do the work. need to get a non-const version of data in order to
      // be able to modify them inside a const function
      auto &this_nonconst = const_cast<FE_SimplexPoly<dim, spacedim> &>(*this);

      std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
        RefinementCase<dim>::isotropic_refinement);
      isotropic_matrices.back().resize(
        GeometryInfo<dim>::n_children(RefinementCase<dim>(refinement_case)),
        FullMatrix<double>(this->n_dofs_per_cell(), this->n_dofs_per_cell()));

      FETools::compute_embedding_matrices(*this, isotropic_matrices, true);

      this_nonconst.prolongation[refinement_case - 1] =
        std::move(isotropic_matrices.back());
    }

  // finally return the matrix
  return this->prolongation[refinement_case - 1][child];
}

The use of GeometryInfo<dim> here is suspicious to begin with, but in any case, FETools::compute_embedding_matrices() only works by considering the geometry of the reference cell (=reference tetrahedron here).

[bangerth]

The test tests/fe/up_and_down.cc checks that this kind of thing actually works by ensuring that prolongation is an embedding operation that preserves the function. We should run this test, or a variant, also on tetrahedra on a number of different meshes.

[bangerth]

@dominiktassilostill FYI

[peterrum]

Picking up the discussion in #16755

When running multigrid_a_01 in 3D with the following modifications:

diff --git a/tests/multigrid-global-coarsening/multigrid_a_01.cc b/tests/multigrid-global-coarsening/multigrid_a_01.cc
index b8b8991feb..424d3cd9f6 100644
--- a/tests/multigrid-global-coarsening/multigrid_a_01.cc
+++ b/tests/multigrid-global-coarsening/multigrid_a_01.cc
@@ -130,7 +130,7 @@ test(const unsigned int n_refinements,
 
   deallog << dim << ' ' << fe_degree_fine << ' ' << n_refinements << ' '
           << (do_simplex_mesh ? "tri " : "quad") << ' '
-          << solver_control.last_step() << std::endl;
+          << solver_control.last_step() << ' ' << dof_handlers[max_level].n_dofs() << std::endl;
 
   static unsigned int counter = 0;
 
@@ -168,6 +168,7 @@ main(int argc, char **argv)
 
   for (const auto ones_on_diagonal : {false, true})
     {
+      if(false)
       for (unsigned int n_refinements = 2; n_refinements <= 4; ++n_refinements)
         for (unsigned int degree = 1; degree <= 4; ++degree)
           test<2>(n_refinements,
@@ -177,6 +178,6 @@ main(int argc, char **argv)
 
       for (unsigned int n_refinements = 2; n_refinements <= 4; ++n_refinements)
         for (unsigned int degree = 1; degree <= 2; ++degree)
-          test<2>(n_refinements, degree, true /*triangle*/, ones_on_diagonal);
+          test<3>(n_refinements, degree, true /*triangle*/, ones_on_diagonal);
     }
 }

I get the following output (the second last numbers are the number of iterations):

OLD:
DEAL:0::3 1 2 tri  3 649
DEAL:0::3 2 2 tri  3 4241
DEAL:0::3 1 3 tri  3 4241
DEAL:0::3 2 3 tri  4 30497
DEAL:0::3 1 4 tri  4 30497
DEAL:0::3 2 4 tri  4 230977
DEAL:0::3 1 2 tri  3 649
DEAL:0::3 2 2 tri  3 4241
DEAL:0::3 1 3 tri  3 4241
DEAL:0::3 2 3 tri  4 30497
DEAL:0::3 1 4 tri  4 30497
DEAL:0::3 2 4 tri  4 230977

NEW:
DEAL:0::3 1 2 tri  3 649
DEAL:0::3 2 2 tri  3 4241
DEAL:0::3 1 3 tri  3 4241
DEAL:0::3 2 3 tri  4 30497
DEAL:0::3 1 4 tri  3 30497
DEAL:0::3 2 4 tri  4 230977
DEAL:0::3 1 2 tri  3 649
DEAL:0::3 2 2 tri  3 4241
DEAL:0::3 1 3 tri  3 4241
DEAL:0::3 2 3 tri  4 30497
DEAL:0::3 1 4 tri  3 30497
DEAL:0::3 2 4 tri  4 230977

The number of iterations look very similar before and after the patch. Is this related to the fact that this is FE_SimplexP and not FE_SimplexDGP?

[kronbichler]

If we run this test one more iteration and start with a subdivided hyper cube with 1 subdivision and an additional global refinement instead, you start to see better iteration counts:

OLD:
DEAL:0::3 1 2 tri  3 645
DEAL:0::3 2 2 tri  4 4233
DEAL:0::3 1 3 tri  4 4233
DEAL:0::3 2 3 tri  4 30481
DEAL:0::3 1 4 tri  5 30481
DEAL:0::3 2 4 tri  5 230945
DEAL:0::3 1 5 tri  6 230945
DEAL:0::3 2 5 tri  6 1797185

NEW:
DEAL:0::3 1 2 tri  3 645
DEAL:0::3 2 2 tri  3 4233
DEAL:0::3 1 3 tri  3 4233
DEAL:0::3 2 3 tri  3 30481
DEAL:0::3 1 4 tri  3 30481
DEAL:0::3 2 4 tri  4 230945
DEAL:0::3 1 5 tri  4 230945
DEAL:0::3 2 5 tri  4 1797185

This shows that either the embedding (for continuous shape functions) is correct also when the geometric configuration changes, or that the errors made in this process are insensitive within multigrid (because it is only a high-frequency error that gets damped on the finer level). I also checked with polynomial degree 3 (where I observed a bug, see #16775) and got 4-5 iterations.
That said, this experiment is either indicating a correct interpolation or it is not affecting multigrid enough to expose the issue with the way the solution gets embedded (I believe it is the latter). @dominiktassilostill do you want to have a look at the embedding indicated above?

[peterrum]

I also checked with polynomial degree 3 (where I observed a bug, see #16775) and got 4-5 iterations.

Right, we now can also test degree 3 👍

[dominiktassilostill]

When testing the new refinement method compared to the old one, I only looked at DG. What I did observe was when going to a polynomial degree of 3 there was an segfault for pure geometric multigrid (so only h-refinement). This also happened for a ch-strategy in ExaDG. For a polynomial degree of 2 or less this didn't occur, also this error didn't show if p-refinement was used first in the multigrid preconditioner. Yet, the error occurs also if chosen_line_tetrahedron is set to 0, i.e. doing the identical refinement as before the patch.
@kronbichler you probably observed the same bug, I will have a look at the embedding.

[peterrum]

I am not sure what the best approach is. But let me document what I am currently thinking:

• use the (potentially new) anisotropic refinement flags in CellAccessor::set_refine_flag()
• if the default is chosen (isotropic_refinement) run the algorithm added in the last PR, i.e., search shortest distance
• store the refinement type in internal::TriangulationImplementation::TriaObjects::refinement_cases
• use this information within multigrid

Manually setting the type in CellAccessor::set_refine_flag() seems to be somewhat useless but this is probably what we will need here:

dealii/include/deal.II/fe/fe_tools.templates.h

Lines 1736 to 1739 in 53bb2ba

	Triangulation<dim, spacedim> tria;
	GridGenerator::reference_cell(tria, reference_cell);
	tria.begin_active()->set_refine_flag(RefinementCase<dim>(ref_case));
	tria.execute_coarsening_and_refinement();

.

[peterrum]

@dominiktassilostill In a first step, what you could do is to prolongate a linear function from a coarse grid to a fine grid and to interpolate a linear function from a fine grid to a coarse grid. I guess the resulting solution is not linear anymore?

[dominiktassilostill]

I checked the prolongation matrices for the 3 different refinement cases and some of them are different. As the first 4 children are the same in each case, they also have the same prolongation matrices, for the rest there are only a few different entries (as most of the entries are 0 anyways) so as @kronbichler suggested the interpolation error is probably small and the tests do not show the effect.

@dominiktassilostill In a first step, what you could do is to prolongate a linear function from a coarse grid to a fine grid and to interpolate a linear function from a fine grid to a coarse grid. I guess the resulting solution is not linear anymore?

@peterrum does this cover the outcome of your suggested experiment or was it intended to check something different?

[kronbichler]

What @peterrum was suggesting is the approximation quality of the transfer, which means considering not the entries in the matrix itself, but the result: To this, you interpolate a linear function on the coarse grid, say a single Simplex. You can visualize it there with ParaView. Then prolongate the solution to the finer grid and visualize the result again. Since the linear function should be in the function space for both cases, it should look exactly the same. I believe that the two cases we introduced in the other PR do reproduce this behavior, but cause some artifacts for p=2; p=3 because the vertices and lines will be aligned differently, so different sections of the function will get stitched together. I believe the problem cannot happen for linear shape functions because there only the vertices matter, which are the same for any of the refinement cases.

Either way, this problem should eventually be fixed. The interesting question would be whether the combination of restriction and prologation "almost" cures the problem in the special case of multigrid, apart from a fine-scale error.

use the (potentially new) anisotropic refinement flags in CellAccessor::set_refine_flag()

Are you referring to a new name here? It is not really an anisotropic refinement, so I would not use the existing data structures for that, but set up a new path that stores these as std::uint8_t for each mother cell.

[bangerth]

The procedure @peterrum describes is in essence what tests/fe/up_and_down.cc and tests/multigrid-global-coarsening/interpolate_up_down_01.cc test. Take a look there -- we already have these kinds of tests, we just don't run them on tetrahedral meshes so far.