From d0e4ca9bd67070e4e5c97299e9b1ecf2564385c1 Mon Sep 17 00:00:00 2001
From: Peter Munch <peterrmuench@gmail.com>
Date: Sat, 7 May 2022 23:25:13 +0200
Subject: [PATCH] Improve asserts in MF/FEEval regarding FENothing

---
 include/deal.II/matrix_free/fe_evaluation.h | 13 +++++
 include/deal.II/matrix_free/matrix_free.h   | 58 ++++++++++++++++++++-
 include/deal.II/matrix_free/util.h          |  7 +++
 3 files changed, 77 insertions(+), 1 deletion(-)

diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h
index 5c80c2f994..120e9dbb31 100644
--- a/include/deal.II/matrix_free/fe_evaluation.h
+++ b/include/deal.II/matrix_free/fe_evaluation.h
@@ -6599,6 +6599,19 @@ FEEvaluation<dim,
   (void)dof_no;
   (void)first_selected_component;
 
+  Assert(
+    this->data->dofs_per_component_on_cell > 0,
+    ExcMessage(
+      "There is nothing useful you can do with an FEEvaluation object with "
+      "FE_Nothing, i.e., without DoFs! If you have passed to "
+      "MatrixFree::reinit() a collection of finite elements also containing "
+      "FE_Nothing, please check - before creating FEEvaluation - the category "
+      "of the current range by calling either "
+      "MatrixFree::get_cell_range_category(range) or "
+      "MatrixFree::get_face_range_category(range). The returned category "
+      "is the index of the active FE, which you can use to exclude "
+      "FE_Nothing."));
+
 #  ifdef DEBUG
   // print error message when the dimensions do not match. Propose a possible
   // fix
diff --git a/include/deal.II/matrix_free/matrix_free.h b/include/deal.II/matrix_free/matrix_free.h
index f977f5d904..da75ab3aaa 100644
--- a/include/deal.II/matrix_free/matrix_free.h
+++ b/include/deal.II/matrix_free/matrix_free.h
@@ -1884,6 +1884,24 @@ public:
   get_face_quadrature(const unsigned int quad_index         = 0,
                       const unsigned int hp_active_fe_index = 0) const;
 
+  /**
+   * Return the category the current batch range of cells was assigned to.
+   * Categories run between the given values in the field
+   * AdditionalData::cell_vectorization_category for non-hp-DoFHandler types
+   * and return the active FE index in the hp-adaptive case.
+   */
+  unsigned int
+  get_cell_range_category(
+    const std::pair<unsigned int, unsigned int> cell_batch_range) const;
+
+  /**
+   * Return the category of the cells on the two sides of the current batch
+   * range of faces.
+   */
+  std::pair<unsigned int, unsigned int>
+  get_face_range_category(
+    const std::pair<unsigned int, unsigned int> face_batch_range) const;
+
   /**
    * Return the category the current batch of cells was assigned to. Categories
    * run between the given values in the field
@@ -1894,7 +1912,7 @@ public:
   get_cell_category(const unsigned int cell_batch_index) const;
 
   /**
-   * Return the category on the cells on the two sides of the current batch of
+   * Return the category of the cells on the two sides of the current batch of
    * faces.
    */
   std::pair<unsigned int, unsigned int>
@@ -2780,6 +2798,44 @@ MatrixFree<dim, Number, VectorizedArrayType>::get_face_quadrature(
 
 
 
+template <int dim, typename Number, typename VectorizedArrayType>
+inline unsigned int
+MatrixFree<dim, Number, VectorizedArrayType>::get_cell_range_category(
+  const std::pair<unsigned int, unsigned int> range) const
+{
+  const auto result = get_cell_category(range.first);
+
+#  ifdef DEBUG
+  for (unsigned int i = range.first; i < range.second; ++i)
+    Assert(result == get_cell_category(i),
+           ExcMessage(
+             "The cell batches of the range have different categories!"));
+#  endif
+
+  return result;
+}
+
+
+
+template <int dim, typename Number, typename VectorizedArrayType>
+inline std::pair<unsigned int, unsigned int>
+MatrixFree<dim, Number, VectorizedArrayType>::get_face_range_category(
+  const std::pair<unsigned int, unsigned int> range) const
+{
+  const auto result = get_face_category(range.first);
+
+#  ifdef DEBUG
+  for (unsigned int i = range.first; i < range.second; ++i)
+    Assert(result == get_face_category(i),
+           ExcMessage(
+             "The face batches of the range have different categories!"));
+#  endif
+
+  return result;
+}
+
+
+
 template <int dim, typename Number, typename VectorizedArrayType>
 inline unsigned int
 MatrixFree<dim, Number, VectorizedArrayType>::get_cell_category(
diff --git a/include/deal.II/matrix_free/util.h b/include/deal.II/matrix_free/util.h
index 323f2093e0..b0b01f6fe1 100644
--- a/include/deal.II/matrix_free/util.h
+++ b/include/deal.II/matrix_free/util.h
@@ -110,6 +110,13 @@ namespace internal
     inline std::pair<Quadrature<dim - 1>, Quadrature<dim - 1>>
     get_unique_face_quadratures(const Quadrature<dim> &quad)
     {
+      AssertThrow(
+        quad.size() > 0,
+        ExcMessage(
+          "There is nothing useful you can do with a MatrixFree/FEEvaluation "
+          "object when using a quadrature formula with zero "
+          "quadrature points!"));
+
       if (dim == 2 || dim == 3)
         {
           for (unsigned int i = 1; i <= 4; ++i)
-- 
2.39.5