From a6ec42c8b2db23f0d198e1cbae5eb0ea6e447e63 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Fri, 14 Jan 2022 22:27:30 +0100 Subject: [PATCH] FEEval: do not cast const away during read_dof_values() --- include/deal.II/matrix_free/fe_evaluation.h | 93 ++++++++++++++----- include/deal.II/matrix_free/type_traits.h | 3 +- .../matrix_free/vector_access_internal.h | 20 ++++ 3 files changed, 91 insertions(+), 25 deletions(-) diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index ac12008a35..0bdcbe6890 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -3270,17 +3270,35 @@ FEEvaluationBase:: namespace internal { + // given a block vector return the underlying vector type + // including constness (specified by bool) + template + struct ConstBlockVectorSelector; + + template + struct ConstBlockVectorSelector + { + using BaseVectorType = const typename VectorType::BlockType; + }; + + template + struct ConstBlockVectorSelector + { + using BaseVectorType = typename VectorType::BlockType; + }; + // allows to select between block vectors and non-block vectors, which // allows to use a unified interface for extracting blocks on block vectors // and doing nothing on usual vectors template - struct BlockVectorSelector - {}; + struct BlockVectorSelector; template struct BlockVectorSelector { - using BaseVectorType = typename VectorType::BlockType; + using BaseVectorType = typename ConstBlockVectorSelector< + VectorType, + std::is_const::value>::BaseVectorType; static BaseVectorType * get_vector_component(VectorType &vec, const unsigned int component) @@ -3329,6 +3347,20 @@ namespace internal } }; + template + struct BlockVectorSelector, false> + { + using BaseVectorType = const VectorType; + + static const BaseVectorType * + get_vector_component(const std::vector &vec, + const unsigned int component) + { + AssertIndexRange(component, vec.size()); + return &vec[component]; + } + }; + template struct BlockVectorSelector, false> { @@ -3342,6 +3374,20 @@ namespace internal return vec[component]; } }; + + template + struct BlockVectorSelector, false> + { + using BaseVectorType = const VectorType; + + static const BaseVectorType * + get_vector_component(const std::vector &vec, + const unsigned int component) + { + AssertIndexRange(component, vec.size()); + return vec[component]; + } + }; } // namespace internal @@ -4183,15 +4229,13 @@ namespace internal template std::pair< std::array::type, - IsBlockVector::type>:: - value>::BaseVectorType *, + VectorType, + IsBlockVector::value>::BaseVectorType *, n_components>, std::array< const std::vector::type, - IsBlockVector::type>::value>:: - BaseVectorType::value_type>> *, + VectorType, + IsBlockVector::value>::BaseVectorType::value_type>> *, n_components>> get_vector_data(VectorType & src, const unsigned int first_index, @@ -4203,32 +4247,33 @@ namespace internal // of components is checked in the internal data std::pair< std::array::type, - IsBlockVector::type>:: - value>::BaseVectorType *, + VectorType, + IsBlockVector::value>::BaseVectorType *, n_components>, std::array< const std::vector< ArrayView::type, - IsBlockVector::type>:: - value>::BaseVectorType::value_type>> *, + VectorType, + IsBlockVector::value>::BaseVectorType::value_type>> *, n_components>> src_data; for (unsigned int d = 0; d < n_components; ++d) src_data.first[d] = internal::BlockVectorSelector< - typename std::remove_const::type, - IsBlockVector::type>::value>:: - get_vector_component( - const_cast::type &>(src), - d + first_index); + VectorType, + IsBlockVector::value>::get_vector_component(src, + d + + first_index); for (unsigned int d = 0; d < n_components; ++d) - src_data.second[d] = get_shared_vector_data(*src_data.first[d], - is_valid_mode_for_sm, - active_fe_index, - dof_info); + src_data.second[d] = get_shared_vector_data( + const_cast::type, + IsBlockVector::type>::value>:: + BaseVectorType &>(*src_data.first[d]), + is_valid_mode_for_sm, + active_fe_index, + dof_info); return src_data; } diff --git a/include/deal.II/matrix_free/type_traits.h b/include/deal.II/matrix_free/type_traits.h index b57c05049a..bce4761ae0 100644 --- a/include/deal.II/matrix_free/type_traits.h +++ b/include/deal.II/matrix_free/type_traits.h @@ -199,7 +199,8 @@ namespace internal { static const bool value = has_begin::value && - (has_local_element::value || is_serial_vector::value) && + (has_local_element::value || + is_serial_vector::type>::value) && std::is_same::value; }; diff --git a/include/deal.II/matrix_free/vector_access_internal.h b/include/deal.II/matrix_free/vector_access_internal.h index 3164675497..7222b3321b 100644 --- a/include/deal.II/matrix_free/vector_access_internal.h +++ b/include/deal.II/matrix_free/vector_access_internal.h @@ -247,10 +247,20 @@ namespace internal VectorizedArrayType *dof_values, std::integral_constant) const { +#ifdef DEBUG + // in debug mode, run non-vectorized version because this path + // has additional checks (e.g., regarding ghosting) + process_dofs_vectorized(dofs_per_cell, + dof_index, + vec, + dof_values, + std::integral_constant()); +#else const Number *vec_ptr = vec.begin() + dof_index; for (unsigned int i = 0; i < dofs_per_cell; ++i, vec_ptr += VectorizedArrayType::size()) dof_values[i].load(vec_ptr); +#endif } @@ -340,7 +350,17 @@ namespace internal VectorizedArrayType &res, std::integral_constant) const { +#ifdef DEBUG + // in debug mode, run non-vectorized version because this path + // has additional checks (e.g., regarding ghosting) + process_dof_gather(indices, + vec, + constant_offset, + res, + std::integral_constant()); +#else res.gather(vec.begin() + constant_offset, indices); +#endif } -- 2.39.5