FEEval::check_vector_compatibility: improve assert message

[peterrum]

The assert message:

Additional information: 
    The parallel layout of the given vector is not compatible with the 
    parallel partitioning in MatrixFree. A potential reason is that 
    you did not use MatrixFree::initialize_dof_vector() to get a
    compatible vector. If you did, the dof_handler_index used there
    and the one passed to the constructor of FEEvaluation do not match.  

is a bit too general. In particular, the hints at the end could be improved by adding actual indices.

This PR tries to improve the assert message. For the following program:

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

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

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

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

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

using namespace dealii;

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

  AssertDimension(argc, 2);

  const unsigned int dim = 2;

  parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
  GridGenerator::hyper_cube(tria);

  FE_Q<dim>            fe_1(1);
  FE_Q<dim>            fe_2(2);
  QGauss<1>            quad(1);
  MappingQGeneric<dim> mapping(1);

  DoFHandler<dim> dof_handler_1(tria);
  dof_handler_1.distribute_dofs(fe_1);

  DoFHandler<dim> dof_handler_2(tria);
  dof_handler_2.distribute_dofs(fe_2);

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

  typename MatrixFree<dim, double>::AdditionalData additional_data;
  additional_data.mapping_update_flags = update_gradients | update_values;

  MatrixFree<dim, double> matrix_free;
  matrix_free.reinit(mapping,
                     std::vector<const DoFHandler<dim> *>{&dof_handler_1,
                                                          &dof_handler_2,
                                                          &dof_handler_1},
                     std::vector<const AffineConstraints<double> *>{
                       &constraints, &constraints, &constraints},
                     quad,
                     additional_data);

  using VectorType = LinearAlgebra::distributed::Vector<double>;

  switch (atoi(argv[1]))
    {
      case 0:
        {
          FEEvaluation<dim, -1, 0, 1, double> phi(matrix_free, 0);
          phi.reinit(0);

          VectorType vec;

          phi.read_dof_values(vec);
          break;
        }
      case 1:
        {
          FEEvaluation<dim, -1, 0, 1, double> phi(matrix_free, 1);
          phi.reinit(0);

          VectorType vec;

          phi.read_dof_values(vec);
          break;
        }
      case 2:
        {
          FEEvaluation<dim, -1, 0, 1, double> phi(matrix_free, 0);
          phi.reinit(0);

          VectorType vec;
          matrix_free.initialize_dof_vector(vec, 1);

          phi.read_dof_values(vec);
          break;
        }
      case 3:
        {
          FEEvaluation<dim, -1, 0, 1, double> phi(matrix_free, 1);
          phi.reinit(0);

          VectorType vec;
          matrix_free.initialize_dof_vector(vec, 0);

          phi.read_dof_values(vec);
          break;
        }
    }
}

the following assert messages are printed:

Additional information: 
    The parallel layout of the given vector is compatible neither with the
    Partitioner of the current FEEvaluation with dof_handler_index=0 nor
    with any Partioner in MatrixFree. A potential reason is that you did
    not use MatrixFree::initialize_dof_vector() to get a compatible
    vector.
    
Additional information: 
    The parallel layout of the given vector is compatible neither with the
    Partitioner of the current FEEvaluation with dof_handler_index=1 nor
    with any Partioner in MatrixFree. A potential reason is that you did
    not use MatrixFree::initialize_dof_vector() to get a compatible
    vector.
    
Additional information: 
    The parallel layout of the given vector is not compatible with the
    Partitioner of the current FEEvaluation with dof_handler_index=0.
    However, the underlying MatrixFree contais Partitioner objects that
    are compatible. They have the following dof_handler_index values: 1.
    Did you want to pass any of these values to the constructor of the
    current FEEvaluation object or did you not use
    MatrixFree::initialize_dof_vector() with dof_handler_index=0 to get a
    compatible vector.        
    
Additional information: 
    The parallel layout of the given vector is not compatible with the
    Partitioner of the current FEEvaluation with dof_handler_index=1.
    However, the underlying MatrixFree contais Partitioner objects that
    are compatible. They have the following dof_handler_index values: 0,
    2. Did you want to pass any of these values to the constructor of the
    current FEEvaluation object or did you not use
    MatrixFree::initialize_dof_vector() with dof_handler_index=1 to get a
    compatible vector.  

Written jointly with @mschreter.

[peterrum]

@kronbichler I had to introduce the dof index in the FEEvaluationData class. Is that fine? Or are there other ways to get the dof index of FEEvaluation?

[kronbichler]

I am not happy to add additional variables to the FEEvaluation base class for the sole purpose of an assertion. I had rather hoped to work towards reducing the footprint of these classes given the performance-critical context in which they are used. But I think it is pretty easy to achieve what you want by letting the get_dof_handler_index() sit in the FEEvaluationBase class instead and implement the function there: Simply run through the matrix_free.n_components() and let it compare pointers of dof_info and matrix_free.get_dof_info(i) until we find the right one (or the pointer is zero).

[kronbichler]

How much code does this include file pull in?

[kronbichler]

It's also a borderline use case that can easily be done differently, in case there is a huge amount of code getting in.

[peterrum]

https://www.boost.org/doc/libs/1_50_0/boost/algorithm/string/join.hpp. I have added an ifdef so that it is only included in release mode. Would that work for you?

[peterrum]

Simply run through the matrix_free.n_components() and let it compare pointers of dof_info and matrix_free.get_dof_info(i) until we find the right one (or the pointer is zero).

I made the modification.

[kronbichler]

I checked and it seems we include boost/algorithm/string/join.hpp already from before via boost/algorithm/string.hpp, which comes in via boost/io/wtk/read.hpp, which in turn comes in via the geometry functionality of boost, which eventually gets its way into our code via grid_tools.h. (Boost pulls in an incredible amount of header files, for example compiling step-37 I see 625,558 lines of code in the final compilation unit after the preprocessor has finished. Probably 65% of that is boost, 20% is Trilinos, the rest is STL or our own code.)

[peterrum]

/rebuild