From 45cf863d2d1dacca033ad79b313d8cec2460dd98 Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Tue, 5 Mar 2019 19:24:58 +0100 Subject: [PATCH] add more vector_access_XYZ internal function in FEEvaluation and improve const-correctness of reader class --- include/deal.II/matrix_free/fe_evaluation.h | 209 ++++++++++++++++-- .../fe_evaluation_type_traits_02.cc | 9 + .../fe_evaluation_type_traits_02.output | 6 + 3 files changed, 205 insertions(+), 19 deletions(-) diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index 5d3b0f4128..5ae7b1905c 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -36,6 +36,7 @@ #include #include + DEAL_II_NAMESPACE_OPEN @@ -3592,6 +3593,96 @@ namespace internal + // same for const access + template ::value, + VectorType>::type * = nullptr> + inline typename VectorType::value_type + vector_access(const VectorType &vec, const unsigned int entry) + { + return vec.local_element(entry); + } + + + + template ::value, + VectorType>::type * = nullptr> + inline void + vector_access_add(VectorType & vec, + const unsigned int entry, + const typename VectorType::value_type &val) + { + vec.add_local_element(entry, val); + } + + + + template ::value, + VectorType>::type * = nullptr> + inline void + vector_access_add(VectorType & vec, + const unsigned int entry, + const typename VectorType::value_type &val) + { + vector_access(vec, entry) += val; + } + + + + template ::value, + VectorType>::type * = nullptr> + inline void + vector_access_add_global(VectorType & vec, + const types::global_dof_index entry, + const typename VectorType::value_type &val) + { + vec.add(entry, val); + } + + + + template ::value, + VectorType>::type * = nullptr> + inline void + vector_access_add_global(VectorType & vec, + const types::global_dof_index entry, + const typename VectorType::value_type &val) + { + vec(entry) += val; + } + + + + template ::value, + VectorType>::type * = nullptr> + inline void + vector_access_set(VectorType & vec, + const unsigned int entry, + const typename VectorType::value_type &val) + { + vec.set_local_element(entry, val); + } + + + + template ::value, + VectorType>::type * = nullptr> + inline void + vector_access_set(VectorType & vec, + const unsigned int entry, + const typename VectorType::value_type &val) + { + vector_access(vec, entry) = val; + } + + + // this is to make sure that the parallel partitioning in VectorType // is really the same as stored in MatrixFree. // version below is when has_partitioners_are_compatible == false @@ -3612,6 +3703,7 @@ namespace internal } + // same as above for has_partitioners_are_compatible == true template < typename VectorType, @@ -3632,17 +3724,28 @@ namespace internal "compatible vector.")); } - // A class to use the same code to read from and write to vector + + + // Below, three classes (VectorReader, VectorSetter, + // VectorDistributorLocalToGlobal) implement the same interface and can be + // used to to read from vector, set elements of a vector and add to elements + // of the vector. + + // 1. A class to read data from vector template struct VectorReader { template void - process_dof(const unsigned int index, VectorType &vec, Number &res) const + process_dof(const unsigned int index, + const VectorType & vec, + Number & res) const { res = vector_access(vec, index); } + + template void process_dofs_vectorized(const unsigned int dofs_per_cell, @@ -3658,11 +3761,12 @@ namespace internal } + template void process_dofs_vectorized(const unsigned int dofs_per_cell, const unsigned int dof_index, - VectorType & vec, + const VectorType & vec, VectorizedArray *dof_values, std::integral_constant) const { @@ -3673,6 +3777,8 @@ namespace internal vec, dof_index + v + i * VectorizedArray::n_array_elements); } + + template void process_dofs_vectorized_transpose(const unsigned int dofs_per_cell, @@ -3688,11 +3794,12 @@ namespace internal } + template void process_dofs_vectorized_transpose(const unsigned int dofs_per_cell, const unsigned int * dof_indices, - VectorType & vec, + const VectorType & vec, VectorizedArray *dof_values, std::integral_constant) const { @@ -3702,6 +3809,8 @@ namespace internal dof_values[d][v] = vector_access(vec, dof_indices[v] + d); } + + // variant where VectorType::value_type is the same as Number -> can call // gather template @@ -3715,12 +3824,14 @@ namespace internal res.gather(vec.begin() + constant_offset, indices); } + + // variant where VectorType::value_type is not the same as Number -> must // manually load the data template void process_dof_gather(const unsigned int * indices, - VectorType & vec, + const VectorType & vec, const unsigned int constant_offset, VectorizedArray &res, std::integral_constant) const @@ -3730,37 +3841,47 @@ namespace internal res[v] = vector_access(vec, indices[v] + constant_offset); } + + template void process_dof_global(const types::global_dof_index index, - VectorType & vec, + const VectorType & vec, Number & res) const { - res = const_cast(vec)(index); + res = vec(index); } + + void pre_constraints(const Number &, Number &res) const { res = Number(); } + + template void process_constraint(const unsigned int index, const Number weight, - VectorType & vec, + const VectorType & vec, Number & res) const { res += weight * vector_access(vec, index); } + + void post_constraints(const Number &sum, Number &write_pos) const { write_pos = sum; } + + void process_empty(VectorizedArray &res) const { @@ -3768,7 +3889,10 @@ namespace internal } }; - // A class to use the same code to read from and write to vector + + + // 2. A class to add values to the vector during + // FEEvaluation::distribute_local_to_global() call template struct VectorDistributorLocalToGlobal { @@ -3776,9 +3900,11 @@ namespace internal void process_dof(const unsigned int index, VectorType &vec, Number &res) const { - vector_access(vec, index) += res; + vector_access_add(vec, index, res); } + + template void process_dofs_vectorized(const unsigned int dofs_per_cell, @@ -3798,6 +3924,8 @@ namespace internal } } + + template void process_dofs_vectorized(const unsigned int dofs_per_cell, @@ -3809,12 +3937,14 @@ namespace internal for (unsigned int i = 0; i < dofs_per_cell; ++i) for (unsigned int v = 0; v < VectorizedArray::n_array_elements; ++v) - vector_access(vec, - dof_index + v + - i * VectorizedArray::n_array_elements) += - dof_values[i][v]; + vector_access_add(vec, + dof_index + v + + i * VectorizedArray::n_array_elements, + dof_values[i][v]); } + + template void process_dofs_vectorized_transpose(const unsigned int dofs_per_cell, @@ -3827,6 +3957,8 @@ namespace internal true, dofs_per_cell, dof_values, dof_indices, vec.begin()); } + + template void process_dofs_vectorized_transpose(const unsigned int dofs_per_cell, @@ -3838,9 +3970,11 @@ namespace internal for (unsigned int d = 0; d < dofs_per_cell; ++d) for (unsigned int v = 0; v < VectorizedArray::n_array_elements; ++v) - vector_access(vec, dof_indices[v] + d) += dof_values[d][v]; + vector_access_add(vec, dof_indices[v] + d, dof_values[d][v]); } + + // variant where VectorType::value_type is the same as Number -> can call // scatter template @@ -3864,6 +3998,8 @@ namespace internal # endif } + + // variant where VectorType::value_type is not the same as Number -> must // manually append all data template @@ -3876,24 +4012,30 @@ namespace internal { for (unsigned int v = 0; v < VectorizedArray::n_array_elements; ++v) - vector_access(vec, indices[v] + constant_offset) += res[v]; + vector_access_add(vec, indices[v] + constant_offset, res[v]); } + + template void process_dof_global(const types::global_dof_index index, VectorType & vec, Number & res) const { - vec(index) += res; + vector_access_add_global(vec, index, res); } + + void pre_constraints(const Number &input, Number &res) const { res = input; } + + template void process_constraint(const unsigned int index, @@ -3901,20 +4043,25 @@ namespace internal VectorType & vec, Number & res) const { - vector_access(vec, index) += weight * res; + vector_access_add(vec, index, weight * res); } + + void post_constraints(const Number &, Number &) const {} + + void process_empty(VectorizedArray &) const {} }; - // A class to use the same code to read from and write to vector + + // 3. A class to set elements of the vector template struct VectorSetter { @@ -3925,6 +4072,8 @@ namespace internal vector_access(vec, index) = res; } + + template void process_dofs_vectorized(const unsigned int dofs_per_cell, @@ -3939,6 +4088,8 @@ namespace internal dof_values[i].store(vec_ptr); } + + template void process_dofs_vectorized(const unsigned int dofs_per_cell, @@ -3956,6 +4107,8 @@ namespace internal dof_values[i][v]; } + + template void process_dofs_vectorized_transpose(const unsigned int dofs_per_cell, @@ -3968,6 +4121,8 @@ namespace internal false, dofs_per_cell, dof_values, dof_indices, vec.begin()); } + + template void process_dofs_vectorized_transpose(const unsigned int dofs_per_cell, @@ -3982,6 +4137,8 @@ namespace internal vector_access(vec, dof_indices[v] + i) = dof_values[i][v]; } + + template void process_dof_gather(const unsigned int * indices, @@ -3993,6 +4150,8 @@ namespace internal res.scatter(indices, vec.begin() + constant_offset); } + + template void process_dof_gather(const unsigned int * indices, @@ -4006,6 +4165,8 @@ namespace internal vector_access(vec, indices[v] + constant_offset) = res[v]; } + + template void process_dof_global(const types::global_dof_index index, @@ -4015,10 +4176,14 @@ namespace internal vec(index) = res; } + + void pre_constraints(const Number &, Number &) const {} + + template void process_constraint(const unsigned int, @@ -4027,15 +4192,21 @@ namespace internal Number &) const {} + + void post_constraints(const Number &, Number &) const {} + + void process_empty(VectorizedArray &) const {} }; + + // 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 diff --git a/tests/matrix_free/fe_evaluation_type_traits_02.cc b/tests/matrix_free/fe_evaluation_type_traits_02.cc index f87bf95c1b..6f57b2aab2 100644 --- a/tests/matrix_free/fe_evaluation_type_traits_02.cc +++ b/tests/matrix_free/fe_evaluation_type_traits_02.cc @@ -116,5 +116,14 @@ main() << internal::has_set_local_element>::value << std::endl; + // now check internal::vector_access_[set|add] wrappers + deallog << "internal::vector_access_set:" << std::endl; + internal::vector_access_set(dummy, 0, 0); + internal::vector_access_set(dummy2, 0, 0); + + deallog << "internal::vector_access_add:" << std::endl; + internal::vector_access_add(dummy, 0, 0); + internal::vector_access_add(dummy2, 0, 0); + deallog << "OK" << std::endl; } diff --git a/tests/matrix_free/fe_evaluation_type_traits_02.output b/tests/matrix_free/fe_evaluation_type_traits_02.output index d6fcec8248..8413d1fad5 100644 --- a/tests/matrix_free/fe_evaluation_type_traits_02.output +++ b/tests/matrix_free/fe_evaluation_type_traits_02.output @@ -5,4 +5,10 @@ DEAL::Dummy2 = 0 DEAL::has_set_local_element: DEAL::Dummy = 1 DEAL::Dummy2 = 0 +DEAL::internal::vector_access_set: +DEAL::Dummy::set_local_element() +DEAL::Dummy2::local_element() +DEAL::internal::vector_access_add: +DEAL::Dummy::add_local_element() +DEAL::Dummy2::local_element() DEAL::OK -- 2.39.5