FEEval: do not use cell as indication if reinit() was called

diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h
index 86e54949e2e91ce574055b40c1162c3a7a95a86e..13aa4761c94392da544eb5cb55fb863deff1de4f 100644 (file)
--- a/include/deal.II/matrix_free/fe_evaluation.h
+++ b/include/deal.II/matrix_free/fe_evaluation.h
@@ -5008,7 +5008,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
   submit_value(const Tensor<1, n_components_, VectorizedArrayType> val_in,
                const unsigned int                                  q_point)
 {
-  Assert(this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
+  Assert(this->is_reinitialized, ExcNotInitialized());
   AssertIndexRange(q_point, this->n_quadrature_points);
   Assert(this->J_value != nullptr,
          internal::ExcMatrixFreeAccessToUninitializedMappingField(
@@ -5046,7 +5046,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
     const Tensor<1, n_components_, Tensor<1, dim, VectorizedArrayType>> grad_in,
     const unsigned int                                                  q_point)
 {
-  Assert(this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
+  Assert(this->is_reinitialized, ExcNotInitialized());
   AssertIndexRange(q_point, this->n_quadrature_points);
   Assert(this->J_value != nullptr,
          internal::ExcMatrixFreeAccessToUninitializedMappingField(
@@ -5158,7 +5158,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
                        hessian_in,
     const unsigned int q_point)
 {
-  Assert(this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
+  Assert(this->is_reinitialized, ExcNotInitialized());
   AssertIndexRange(q_point, this->n_quadrature_points);
   Assert(this->J_value != nullptr,
          internal::ExcMatrixFreeAccessToUninitializedMappingField(
@@ -5309,7 +5309,7 @@ inline Tensor<1, n_components_, VectorizedArrayType>
 FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
   integrate_value() const
 {
-  Assert(this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
+  Assert(this->is_reinitialized, ExcNotInitialized());
 #  ifdef DEBUG
   Assert(this->values_quad_submitted == true,
          internal::ExcAccessToUninitializedField());
@@ -5646,7 +5646,7 @@ FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>::submit_value(
   const VectorizedArrayType val_in,
   const unsigned int        q_point)
 {
-  Assert(this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
+  Assert(this->is_reinitialized, ExcNotInitialized());
   AssertIndexRange(q_point, this->n_quadrature_points);
   Assert(this->J_value != nullptr,
          internal::ExcMatrixFreeAccessToUninitializedMappingField(
@@ -5699,7 +5699,7 @@ FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>::
   submit_gradient(const Tensor<1, dim, VectorizedArrayType> grad_in,
                   const unsigned int                        q_point)
 {
-  Assert(this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
+  Assert(this->is_reinitialized, ExcNotInitialized());
   AssertIndexRange(q_point, this->n_quadrature_points);
   Assert(this->J_value != nullptr,
          internal::ExcMatrixFreeAccessToUninitializedMappingField(
@@ -6016,7 +6016,7 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
   submit_divergence(const VectorizedArrayType div_in,
                     const unsigned int        q_point)
 {
-  Assert(this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
+  Assert(this->is_reinitialized, ExcNotInitialized());
   AssertIndexRange(q_point, this->n_quadrature_points);
   Assert(this->J_value != nullptr,
          internal::ExcMatrixFreeAccessToUninitializedMappingField(
@@ -6078,7 +6078,7 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
   // could have used base class operator, but that involves some overhead
   // which is inefficient. it is nice to have the symmetric tensor because
   // that saves some operations
-  Assert(this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
+  Assert(this->is_reinitialized, ExcNotInitialized());
   AssertIndexRange(q_point, this->n_quadrature_points);
   Assert(this->J_value != nullptr,
          internal::ExcMatrixFreeAccessToUninitializedMappingField(
@@ -6372,7 +6372,7 @@ FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>::submit_value(
   const VectorizedArrayType val_in,
   const unsigned int        q_point)
 {
-  Assert(this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
+  Assert(this->is_reinitialized, ExcNotInitialized());
   AssertIndexRange(q_point, this->n_quadrature_points);
 #  ifdef DEBUG
   this->values_quad_submitted = true;
@@ -6421,7 +6421,7 @@ FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>::submit_gradient(
   const VectorizedArrayType grad_in,
   const unsigned int        q_point)
 {
-  Assert(this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
+  Assert(this->is_reinitialized, ExcNotInitialized());
   AssertIndexRange(q_point, this->n_quadrature_points);
 #  ifdef DEBUG
   this->gradients_quad_submitted = true;
@@ -6907,6 +6907,7 @@ FEEvaluation<dim,
     this->cell_ids[i] = cell_index * VectorizedArrayType::size() + i;
 
 #  ifdef DEBUG
+  this->is_reinitialized           = true;
   this->dof_values_initialized     = false;
   this->values_quad_initialized    = false;
   this->gradients_quad_initialized = false;
@@ -6945,6 +6946,10 @@ FEEvaluation<dim,
     cell->get_mg_dof_indices(this->local_dof_indices);
   else
     cell->get_dof_indices(this->local_dof_indices);
+
+#  ifdef DEBUG
+  this->is_reinitialized = true;
+#  endif
 }
 
 
@@ -6971,6 +6976,10 @@ FEEvaluation<dim,
                     "instead"));
   Assert(this->mapped_geometry.get() != 0, ExcNotInitialized());
   this->mapped_geometry->reinit(cell);
+
+#  ifdef DEBUG
+  this->is_reinitialized = true;
+#  endif
 }
 
 
@@ -7640,6 +7649,7 @@ FEFaceEvaluation<dim,
        ->normals_times_jacobians[!this->is_interior_face][offsets];
 
 #  ifdef DEBUG
+  this->is_reinitialized           = true;
   this->dof_values_initialized     = false;
   this->values_quad_initialized    = false;
   this->gradients_quad_initialized = false;
@@ -7776,6 +7786,7 @@ FEFaceEvaluation<dim,
        .normals_times_jacobians[!this->is_interior_face][offsets];
 
 #  ifdef DEBUG
+  this->is_reinitialized           = true;
   this->dof_values_initialized     = false;
   this->values_quad_initialized    = false;
   this->gradients_quad_initialized = false;

diff --git a/include/deal.II/matrix_free/fe_evaluation_data.h b/include/deal.II/matrix_free/fe_evaluation_data.h
index de3cb3753535a152fa04ddc3e7f504d6bd589c5e..ee815a8eba61dfa6db5acf057f318dbbbabe40cc 100644 (file)
--- a/include/deal.II/matrix_free/fe_evaluation_data.h
+++ b/include/deal.II/matrix_free/fe_evaluation_data.h
@@ -472,7 +472,7 @@ public:
   get_cell_ids() const
   {
     // implemented inline to avoid compilation problems on Windows
-    Assert(cell != numbers::invalid_unsigned_int, ExcNotInitialized());
+    Assert(is_reinitialized, ExcNotInitialized());
     return cell_ids;
   }
 
@@ -484,7 +484,7 @@ public:
   get_cell_or_face_batch_id() const
   {
     // implemented inline to avoid compilation problems on Windows
-    Assert(cell != numbers::invalid_unsigned_int, ExcNotInitialized());
+    Assert(is_reinitialized, ExcNotInitialized());
     return cell;
   }
 
@@ -496,7 +496,7 @@ public:
   get_cell_or_face_ids() const
   {
     // implemented inline to avoid compilation problems on Windows
-    Assert(cell != numbers::invalid_unsigned_int, ExcNotInitialized());
+    Assert(is_reinitialized, ExcNotInitialized());
     if (!is_face || dof_access_index ==
                       internal::MatrixFreeFunctions::DoFInfo::dof_access_cell)
       return cell_ids;
@@ -705,6 +705,12 @@ protected:
    */
   Number *hessians_quad;
 
+  /**
+   * Debug information to track whether the reinit() function of derived classes
+   * was called.
+   */
+  bool is_reinitialized;
+
   /**
    * Debug information to track whether dof values have been initialized
    * before accessed. Used to control exceptions when uninitialized data is
@@ -895,6 +901,7 @@ inline FEEvaluationData<dim, Number, is_face>::FEEvaluationData(
   , quadrature_weights(
       descriptor != nullptr ? descriptor->quadrature_weights.begin() : nullptr)
 #  ifdef DEBUG
+  , is_reinitialized(false)
   , dof_values_initialized(false)
   , values_quad_initialized(false)
   , gradients_quad_initialized(false)
@@ -944,6 +951,7 @@ inline FEEvaluationData<dim, Number, is_face>::FEEvaluationData(
   , is_interior_face(true)
   , dof_access_index(internal::MatrixFreeFunctions::DoFInfo::dof_access_cell)
   , mapped_geometry(mapped_geometry)
+  , is_reinitialized(false)
 {
   mapping_data = &mapped_geometry->get_data_storage();
   jacobian     = mapped_geometry->get_data_storage().jacobians[0].begin();
@@ -974,6 +982,7 @@ FEEvaluationData<dim, Number, is_face>::operator=(const FEEvaluationData &other)
   quadrature_weights = other.quadrature_weights;
 
 #  ifdef DEBUG
+  is_reinitialized           = false;
   dof_values_initialized     = false;
   values_quad_initialized    = false;
   gradients_quad_initialized = false;
@@ -1242,7 +1251,7 @@ template <int dim, typename Number, bool is_face>
 inline internal::MatrixFreeFunctions::GeometryType
 FEEvaluationData<dim, Number, is_face>::get_cell_type() const
 {
-  Assert(cell != numbers::invalid_unsigned_int, ExcNotInitialized());
+  Assert(is_reinitialized, ExcNotInitialized());
   return cell_type;
 }