Unresolved matrix-free simplex projection problems

[drwells]

Part of #11166. I tried to do a matrix-free projection on master (after writing #11203) and things don't quite work. Here's the code:

#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/function_lib.h>

#include <deal.II/distributed/shared_tria.h>

#include <deal.II/dofs/dof_handler.h>

#include <deal.II/fe/mapping_fe.h>

#include <deal.II/grid/grid_in.h>

#include <deal.II/lac/la_parallel_vector.h>

#include <deal.II/matrix_free/fe_evaluation.h>
#include <deal.II/matrix_free/matrix_free.h>

#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools_project.h>

#include <deal.II/simplex/fe_lib.h>
#include <deal.II/simplex/grid_generator.h>
#include <deal.II/simplex/quadrature_lib.h>

using namespace dealii;

template <int dim, int spacedim = dim>
void
test(const FiniteElement<dim, spacedim> &fe)
{
  ConditionalOStream pcout(std::cout,
                           Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0);
  parallel::shared::Triangulation<dim, spacedim> tria(MPI_COMM_WORLD);
  if (dim == 2)
  {
      GridGenerator::subdivided_hyper_cube_with_simplices(tria, 10);
  }
  if (dim == 3)
  {
    GridGenerator::subdivided_hyper_cube_with_simplices(tria, 100);
#if 0
    GridIn<dim> gi;
    gi.attach_triangulation(tria);
    pcout << "Reading mesh:\n";
    gi.read_exodusii("whole_heart_with_fibers_v5.e");
    pcout << "Reading mesh: Done.\n";
#endif
  }
  pcout << "Number of global active cells: " << tria.n_global_active_cells() << '\n';
  pcout << "Number of local active cells:  " << tria.n_active_cells() << '\n';

  pcout << "Distributing DoFs:\n";
  DoFHandler<dim> dof_handler(tria);
  dof_handler.distribute_dofs(fe);
  pcout << "Distributing DoFs: Done.\n";

  pcout << "Number of global DoFs: " << dof_handler.n_dofs() << '\n';
  pcout << "Number of local DoFs:  " << dof_handler.n_locally_owned_dofs() << '\n';

  MappingFE<dim> mapping(Simplex::FE_P<dim>(1));

  MatrixFree<dim, double> matrix_free;

  AffineConstraints<double> constraints;
  constraints.close();

  Simplex::QGauss<dim> quadrature(3);

  pcout << "Initializing MatrixFree:\n";
  matrix_free.reinit(mapping, dof_handler, constraints, quadrature);
  pcout << "Initializing MatrixFree: Done.\n";

  LinearAlgebra::distributed::Vector<double> dst, src;

  matrix_free.initialize_dof_vector(dst);
  matrix_free.initialize_dof_vector(src);

  pcout << "Projecting:\n";
  VectorTools::project(mapping, dof_handler, constraints, quadrature,
                       Functions::ConstantFunction<dim>(42.0), dst);
  dst.update_ghost_values();
  pcout << "Projecting: Done.\n";
  for (const auto &value : dst)
  {
      if (std::abs(value - 42.0) > 1e-6)
          pcout << value << '\n';
  }

  pcout << "Building Patches:\n";
  DataOut<dim> data_out;
  data_out.attach_dof_handler(dof_handler);
  data_out.add_data_vector(dst, "u");
  data_out.build_patches(mapping);
  pcout << "Building Patches: Done.\n";

  pcout << "Writing PVTU files:\n";
  data_out.write_vtu_with_pvtu_record("./", "solution", 0, MPI_COMM_WORLD, 2, 8);
  pcout << "Writing PVTU files: Done.\n";
}

int
main(int argc, char **argv)
{
  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);

  const unsigned int dim = 2;
  test<dim>(Simplex::FE_P<dim>(2)); // 1 doesn't work (nans at step 1)
}

If we switch everything over to using FE_Q and quadrilaterals things are fine but the values printed out here are all off: most are around 30-40.

[peterrum]

I had a quick look at the implementation of VectorTools::project and one line striked me in particular - and is definitely not correct in the context of simplex:

dealii/include/deal.II/numerics/vector_tools_project.templates.h

Line 193 in 203108e

QGauss<1>(dof.get_fe().degree + 2),

@kronbichler Do you have an idea how we could "fix" this implicit assumption? Let me verify tomorrow if this is the only issue.

[drwells]

One possibility would be to let a finite element (or a finite element collection) return the correct quadrature rules for evaluating mass matrices themselves.

[drwells]

Either way - hard-coding the correct quadrature rules into vector_tools_project.templates.h fixes the issue.

[peterrum]

But only for fe_degree=1? I don't think i have implemented Simplex::QGauss<1>(4).

[drwells]

I was able to project a constant correctly for TET10 (I set everything to Simplex::QGauss<dim>(3))

[peterrum]

Excellent 👍