Fix invalid data access in MatrixFreeTools::compute_diagonal

diff --git a/include/deal.II/matrix_free/tools.h b/include/deal.II/matrix_free/tools.h
index a3b7ac947e8cb7a207f346528e07e35a35aa5f94..0afcd0992be07996f2a3d1b38145d76b2b950a36 100644 (file)
--- a/include/deal.II/matrix_free/tools.h
+++ b/include/deal.II/matrix_free/tools.h
@@ -444,10 +444,12 @@ namespace MatrixFreeTools
 
       ComputeDiagonalHelper()
         : phi(nullptr)
+        , dofs_per_component(0)
       {}
 
       ComputeDiagonalHelper(const ComputeDiagonalHelper &)
         : phi(nullptr)
+        , dofs_per_component(0)
       {}
 
       void
@@ -460,9 +462,11 @@ namespace MatrixFreeTools
       {
         // if we are in hp mode and the number of unknowns changed, we must
         // clear the map of entries
-        if (this->phi != nullptr &&
-            this->phi->dofs_per_component != phi.dofs_per_component)
-          locally_relevant_constraints_hn_map.clear();
+        if (dofs_per_component != phi.dofs_per_component)
+          {
+            locally_relevant_constraints_hn_map.clear();
+            dofs_per_component = phi.dofs_per_component;
+          }
         this->phi = φ
       }
 
@@ -474,8 +478,7 @@ namespace MatrixFreeTools
         // STEP 1: get relevant information from FEEvaluation
         const auto &       dof_info        = phi->get_dof_info();
         const unsigned int n_fe_components = dof_info.start_components.back();
-        const unsigned int dofs_per_component = phi->dofs_per_component;
-        const auto &       matrix_free        = phi->get_matrix_free();
+        const auto &       matrix_free     = phi->get_matrix_free();
 
         // if we have a block vector with components with the same DoFHandler,
         // each component is described with same set of constraints and
@@ -799,7 +802,6 @@ namespace MatrixFreeTools
         // first guess
         const unsigned int degree =
           phi->get_shape_info().data.front().fe_degree;
-        const unsigned int dofs_per_component = phi->dofs_per_component;
         constraints_hn.reserve(Utilities::pow(degree + 1, dim - 1));
 
         // 1) collect hanging-node constraints for cell assuming
@@ -873,9 +875,9 @@ namespace MatrixFreeTools
         const unsigned int n_fe_components =
           phi->get_dof_info().start_components.back();
         const unsigned int comp =
-          n_fe_components == 1 ? i / phi->dofs_per_component : 0;
+          n_fe_components == 1 ? i / dofs_per_component : 0;
         const unsigned int i_comp =
-          n_fe_components == 1 ? (i % phi->dofs_per_component) : i;
+          n_fe_components == 1 ? (i % dofs_per_component) : i;
 
         // apply local constraint matrix from left and from right:
         // loop over all rows of transposed constrained matrix
@@ -894,10 +896,9 @@ namespace MatrixFreeTools
                 // apply constraint matrix from the left
                 Number temp = 0.0;
                 for (unsigned int k = c_pool.row[j]; k < c_pool.row[j + 1]; ++k)
-                  temp +=
-                    c_pool.val[k] *
-                    phi->begin_dof_values()[comp * phi->dofs_per_component +
-                                            c_pool.col[k]][v];
+                  temp += c_pool.val[k] *
+                          phi->begin_dof_values()[comp * dofs_per_component +
+                                                  c_pool.col[k]][v];
 
                 // apply constraint matrix from the right
                 diagonals_local_constrained
@@ -944,6 +945,8 @@ namespace MatrixFreeTools
                    Number,
                    VectorizedArrayType> *phi;
 
+      unsigned int dofs_per_component;
+
       unsigned int i;
 
       std::array<internal::LocalCSR<Number>, n_lanes> c_pools;