From 4739d73d258e80ff788c5ad01eb1e0d7e20b5014 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Fri, 23 Jun 2023 11:28:58 +0200 Subject: [PATCH] Fix invalid data access in MatrixFreeTools::compute_diagonal --- include/deal.II/matrix_free/tools.h | 27 +++++++++++++++------------ 1 file changed, 15 insertions(+), 12 deletions(-) diff --git a/include/deal.II/matrix_free/tools.h b/include/deal.II/matrix_free/tools.h index a3b7ac947e..0afcd0992b 100644 --- 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, n_lanes> c_pools; -- 2.39.5