FEEval::check_vector_compatibility: improve assert message

include/deal.II/matrix_free/fe_evaluation.h

@@ -2857,6 +2857,7 @@ inline FEEvaluationBase<dim,
                                                      active_quad_index,
                                                      face_type),
       is_interior_face,
+      dof_no,
       quad_no,
       first_selected_component)
   , scratch_data_array(matrix_free.acquire_scratch_data())
@@ -3279,11 +3280,15 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
                           "FEEvaluation. In that case, you must pass an "
                           "std::vector<VectorType> or a BlockVector to " +
                           "read_dof_values and distribute_local_to_global."));
-        internal::check_vector_compatibility(*src[comp], *this->dof_info);
+        internal::check_vector_compatibility(*src[comp],
+                                             *this->matrix_free,
+                                             this->get_dof_handler_index());
       }
   else
     {
-      internal::check_vector_compatibility(*src[0], *this->dof_info);
+      internal::check_vector_compatibility(*src[0],
+                                           *this->matrix_free,
+                                           this->get_dof_handler_index());
     }
 
   // Case 2: contiguous indices which use reduced storage of indices and can

include/deal.II/matrix_free/fe_evaluation_data.h

@@ -380,6 +380,12 @@ class FEEvaluationData
   const std::vector<unsigned int> &
   get_internal_dof_numbering() const;
 
+  /**
+   * Return the number of the DoFHandler of the present cell.
+   */
+  unsigned int
+  get_dof_handler_index() const;
+
   /**
    * Return the number of the quadrature formula of the present cell.
    */
@@ -641,6 +647,7 @@ class FEEvaluationData
    */
   FEEvaluationData(const InitializationData &initialization_data,
                    const bool                is_interior_face,
+                   const unsigned int        dof_no,
                    const unsigned int        quad_no,
                    const unsigned int        first_selected_component);
 
@@ -678,6 +685,13 @@ class FEEvaluationData
    */
   const MappingInfoStorageType *mapping_data;
 
+  /**
+   * The number of the DoFHandler of the present cell among all
+   * DoFHandler objects available in the MatrixFree objects pass to derived
+   * classes.
+   */
+  const unsigned int dof_no;
+
   /**
    * The number of the quadrature formula of the present cell among all
    * quadrature formulas available in the MatrixFree objects pass to derived
@@ -997,6 +1011,7 @@ inline FEEvaluationData<dim, Number, is_face>::FEEvaluationData(
       InitializationData{&shape_info, nullptr, nullptr, 0, 0, nullptr},
       is_interior_face,
       0,
+      0,
       0)
 {}
 
@@ -1005,11 +1020,13 @@ template <int dim, typename Number, bool is_face>
 inline FEEvaluationData<dim, Number, is_face>::FEEvaluationData(
   const InitializationData &initialization_data,
   const bool                is_interior_face,
+  const unsigned int        dof_no,
   const unsigned int        quad_no,
   const unsigned int        first_selected_component)
   : data(initialization_data.shape_info)
   , dof_info(initialization_data.dof_info)
   , mapping_data(initialization_data.mapping_data)
+  , dof_no(dof_no)
   , quad_no(quad_no)
   , n_fe_components(dof_info != nullptr ? dof_info->start_components.back() : 0)
   , first_selected_component(first_selected_component)
@@ -1060,6 +1077,7 @@ inline FEEvaluationData<dim, Number, is_face>::FEEvaluationData(
   : data(nullptr)
   , dof_info(nullptr)
   , mapping_data(nullptr)
+  , dof_no(numbers::invalid_unsigned_int)
   , quad_no(numbers::invalid_unsigned_int)
   , n_fe_components(n_fe_components)
   , first_selected_component(first_selected_component)
@@ -1097,6 +1115,7 @@ template <int dim, typename Number, bool is_face>
 inline FEEvaluationData<dim, Number, is_face> &
 FEEvaluationData<dim, Number, is_face>::operator=(const FEEvaluationData &other)
 {
+  AssertDimension(dof_no, other.dof_no);
   AssertDimension(quad_no, other.quad_no);
   AssertDimension(n_fe_components, other.n_fe_components);
   AssertDimension(first_selected_component, other.first_selected_component);
@@ -1462,6 +1481,15 @@ FEEvaluationData<dim, Number, is_face>::get_internal_dof_numbering() const
 
 
 
+template <int dim, typename Number, bool is_face>
+inline unsigned int
+FEEvaluationData<dim, Number, is_face>::get_dof_handler_index() const
+{
+  return dof_no;
+}
+
+
+
 template <int dim, typename Number, bool is_face>
 inline unsigned int
 FEEvaluationData<dim, Number, is_face>::get_quadrature_index() const

include/deal.II/matrix_free/vector_access_internal.h

@@ -24,6 +24,8 @@
 #include <deal.II/matrix_free/dof_info.h>
 #include <deal.II/matrix_free/type_traits.h>
 
+#include <boost/algorithm/string/join.hpp>
+
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -169,42 +171,94 @@ namespace internal
   // version below is when has_partitioners_are_compatible == false
   // FIXME: this is incorrect for PETSc/Trilinos MPI vectors
   template <
+    int dim,
+    typename Number,
+    typename VectorizedArrayType,
     typename VectorType,
     typename std::enable_if<!has_partitioners_are_compatible<VectorType>,
                             VectorType>::type * = nullptr>
   inline void
   check_vector_compatibility(
-    const VectorType &                            vec,
-    const internal::MatrixFreeFunctions::DoFInfo &dof_info)
+    const VectorType &                                  vec,
+    const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+    const unsigned int                                  dof_index)
   {
     (void)vec;
-    (void)dof_info;
+    (void)matrix_free;
+    (void)dof_index;
 
-    AssertDimension(vec.size(), dof_info.vector_partitioner->size());
+    AssertDimension(
+      vec.size(),
+      matrix_free.get_dof_info(dof_index).vector_partitioner->size());
   }
 
 
 
   // same as above for has_partitioners_are_compatible == true
-  template <typename VectorType,
+  template <int dim,
+            typename Number,
+            typename VectorizedArrayType,
+            typename VectorType,
             typename std::enable_if<has_partitioners_are_compatible<VectorType>,
                                     VectorType>::type * = nullptr>
   inline void
   check_vector_compatibility(
-    const VectorType &                            vec,
-    const internal::MatrixFreeFunctions::DoFInfo &dof_info)
+    const VectorType &                                  vec,
+    const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+    const unsigned int                                  dof_index)
   {
     (void)vec;
-    (void)dof_info;
-    Assert(vec.partitioners_are_compatible(*dof_info.vector_partitioner),
-           ExcMessage(
-             "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."));
+    (void)matrix_free;
+    (void)dof_index;
+
+#if DEBUG
+    if (vec.partitioners_are_compatible(
+          *matrix_free.get_dof_info(dof_index).vector_partitioner) == false)
+      {
+        std::vector<std::string> dof_indices_with_compatible_partitioners;
+
+        for (unsigned int i = 0; i < matrix_free.n_components(); ++i)
+          if (vec.partitioners_are_compatible(
+                *matrix_free.get_dof_info(i).vector_partitioner))
+            dof_indices_with_compatible_partitioners.push_back(
+              std::to_string(i));
+
+        if (dof_indices_with_compatible_partitioners.empty())
+          {
+            Assert(false,
+                   ExcMessage("The parallel layout of the given vector is "
+                              "compatible neither with the Partitioner of the "
+                              "current FEEvaluation with dof_handler_index=" +
+                              std::to_string(dof_index) +
+                              " nor with any Partitioner in MatrixFree. A "
+                              "potential reason is that you did not use "
+                              "MatrixFree::initialize_dof_vector() to get a "
+                              "compatible vector."));
+          }
+        else
+          {
+            Assert(
+              false,
+              ExcMessage(
+                "The parallel layout of the given vector is "
+                "not compatible with the Partitioner of the "
+                "current FEEvaluation with dof_handler_index=" +
+                std::to_string(dof_index) +
+                ". However, the underlying "
+                "MatrixFree contains Partitioner objects that are compatible. "
+                "They have the following dof_handler_index values: " +
+                boost::algorithm::join(dof_indices_with_compatible_partitioners,
+                                       ", ") +
+                ". 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=" +
+                std::to_string(dof_index) +
+                " to get a "
+                "compatible vector?"));
+          }
+      }
+#endif
   }