From 128f0e24ecd2f6b470d81f518f78e7befbff7f54 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Sat, 18 Feb 2017 15:07:04 +0100 Subject: [PATCH] Simplify vector access in FEEvaluation. --- include/deal.II/matrix_free/fe_evaluation.h | 383 +++++--------------- 1 file changed, 81 insertions(+), 302 deletions(-) diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index a30c0c7fd1..df343f87c9 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -184,29 +184,17 @@ public: * MatrixFree object and lead to a structure that does not effectively use * vectorization in the evaluate routines based on these values (instead, * VectorizedArray::n_array_elements same copies are worked on). + * + * If the given vector template class is a block vector (determined through + * the template function 'IsBlockVector::value', which checks + * for vectors derived from dealii::BlockVectorBase) or an + * std::vector or std::vector, this function reads + * @p n_components blocks from the block vector starting at the index + * @p first_index. For non-block vectors, @p first_index is ignored. */ template - void read_dof_values (const VectorType &src); - - /** - * For a collection of several vector @p src, read out the values on the - * degrees of freedom of the current cell for @p n_components (template - * argument), starting at @p first_index, and store them internally. Similar - * functionality as the function ConstraintMatrix::read_dof_values. Note - * that if vectorization is enabled, the DoF values for several cells are - * set. - */ - template - void read_dof_values (const std::vector &src, - const unsigned int first_index=0); - - /** - * Reads data from several vectors. Same as other function with std::vector, - * but accepts a vector of pointers to vectors. - */ - template - void read_dof_values (const std::vector &src, - const unsigned int first_index=0); + void read_dof_values (const VectorType &src, + const unsigned int first_index = 0); /** * For the vector @p src, read out the values on the degrees of freedom of @@ -228,29 +216,17 @@ public: * MatrixFree object and lead to a structure that does not effectively use * vectorization in the evaluate routines based on these values (instead, * VectorizedArray::n_array_elements same copies are worked on). + * + * If the given vector template class is a block vector (determined through + * the template function 'IsBlockVector::value', which checks + * for vectors derived from dealii::BlockVectorBase) or an + * std::vector or std::vector, this function reads + * @p n_components blocks from the block vector starting at the index + * @p first_index. For non-block vectors, @p first_index is ignored. */ template - void read_dof_values_plain (const VectorType &src); - - /** - * For a collection of several vector @p src, read out the values on the - * degrees of freedom of the current cell for @p n_components (template - * argument), starting at @p first_index, and store them internally. Similar - * functionality as the function DoFAccessor::read_dof_values. Note that if - * vectorization is enabled, the DoF values for several cells are set. - */ - template - void read_dof_values_plain (const std::vector &src, - const unsigned int first_index=0); - - /** - * Reads data from several vectors without resolving constraints. Same as - * other function with std::vector, but accepts a vector of pointers to - * vectors. - */ - template - void read_dof_values_plain (const std::vector &src, - const unsigned int first_index=0); + void read_dof_values_plain (const VectorType &src, + const unsigned int first_index = 0); /** * Takes the values stored internally on dof values of the current cell and @@ -267,37 +243,27 @@ public: * MatrixFree object and lead to a structure that does not effectively use * vectorization in the evaluate routines based on these values (instead, * VectorizedArray::n_array_elements same copies are worked on). + * + * If the given vector template class is a block vector (determined through + * the template function 'IsBlockVector::value', which checks + * for vectors derived from dealii::BlockVectorBase) or an + * std::vector or std::vector, this function + * writes to @p n_components blocks of the block vector starting at the + * index @p first_index. For non-block vectors, @p first_index is ignored. */ template - void distribute_local_to_global (VectorType &dst) const; - - /** - * Takes the values stored internally on dof values of the current cell for - * a vector-valued problem consisting of @p n_components (template argument) - * and sums them into the collection of vectors vector @p dst, starting at - * index @p first_index. The function also applies constraints during the - * write operation. The functionality is hence similar to the function - * ConstraintMatrix::distribute_local_to_global. If vectorization is - * enabled, the DoF values for several cells are used. - */ - template - void distribute_local_to_global (std::vector &dst, - const unsigned int first_index=0) const; - - /** - * Write data to several vectors. Same as other function with std::vector, - * but accepts a vector of pointers to vectors. - */ - template - void distribute_local_to_global (std::vector &dst, - const unsigned int first_index=0) const; + void distribute_local_to_global (VectorType &dst, + const unsigned int first_index = 0) const; /** * Takes the values stored internally on dof values of the current cell and - * sums them into the vector @p dst. The function also applies constraints - * during the write operation. The functionality is hence similar to the - * function ConstraintMatrix::distribute_local_to_global. Note that if - * vectorization is enabled, the DoF values for several cells are used. + * writes them into the vector @p dst. The function skips the degrees of + * freedom which are constrained. As opposed to the + * distribute_local_to_global method, the old values at the position given + * by the current cell are overwritten. Thus, if a degree of freedom is + * associated to more than one cell (as usual in continuous finite + * elements), the values will be overwritten and only the value written last + * is retained. * * If this class was constructed without a MatrixFree object and the * information is acquired on the fly through a @@ -307,30 +273,17 @@ public: * MatrixFree object and lead to a structure that does not effectively use * vectorization in the evaluate routines based on these values (instead, * VectorizedArray::n_array_elements same copies are worked on). + * + * If the given vector template class is a block vector (determined through + * the template function 'IsBlockVector::value', which checks + * for vectors derived from dealii::BlockVectorBase) or an + * std::vector or std::vector, this function + * writes to @p n_components blocks of the block vector starting at the + * index @p first_index. For non-block vectors, @p first_index is ignored. */ template - void set_dof_values (VectorType &dst) const; - - /** - * Takes the values stored internally on dof values of the current cell for - * a vector-valued problem consisting of @p n_components (template argument) - * and sums them into the collection of vectors vector @p dst, starting at - * index @p first_index. The function also applies constraints during the - * write operation. The functionality is hence similar to the function - * ConstraintMatrix::distribute_local_to_global. Note that if vectorization - * is enabled, the DoF values for several cells are used. - */ - template - void set_dof_values (std::vector &dst, - const unsigned int first_index=0) const; - - /** - * Write data to several vectors. Same as other function with std::vector, - * but accepts a vector of pointers to vectors. - */ - template - void set_dof_values (std::vector &dst, - const unsigned int first_index=0) const; + void set_dof_values (VectorType &dst, + const unsigned int first_index = 0) const; //@} @@ -2772,6 +2725,32 @@ namespace internal return &vec; } }; + + template + struct BlockVectorSelector,false> + { + typedef VectorType BaseVectorType; + + static BaseVectorType *get_vector_component (std::vector &vec, + const unsigned int component) + { + AssertIndexRange (component, vec.size()); + return &vec[component]; + } + }; + + template + struct BlockVectorSelector,false> + { + typedef VectorType BaseVectorType; + + static BaseVectorType *get_vector_component (std::vector &vec, + const unsigned int component) + { + AssertIndexRange (component, vec.size()); + return vec[component]; + } + }; } @@ -3125,70 +3104,15 @@ template inline void FEEvaluationBase -::read_dof_values (const VectorType &src) +::read_dof_values (const VectorType &src, + const unsigned int first_index) { // select between block vectors and non-block vectors. Note that the number // of components is checked in the internal data typename internal::BlockVectorSelector::value>::BaseVectorType *src_data[n_components]; for (unsigned int d=0; d::value>::get_vector_component(const_cast(src), d); - - internal::VectorReader reader; - read_write_operation (reader, src_data); - -#ifdef DEBUG - dof_values_initialized = true; -#endif -} - - - -template -template -inline -void -FEEvaluationBase -::read_dof_values (const std::vector &src, - const unsigned int first_index) -{ - AssertIndexRange (first_index, src.size()); - Assert (n_fe_components == 1, ExcNotImplemented()); - Assert ((n_fe_components == 1 ? - (first_index+n_components <= src.size()) : true), - ExcIndexRange (first_index + n_components_, 0, src.size())); - - VectorType *src_data [n_components]; - for (unsigned int comp=0; comp(&src[comp+first_index]); - - internal::VectorReader reader; - read_write_operation (reader, src_data); - -#ifdef DEBUG - dof_values_initialized = true; -#endif -} - - - -template -template -inline -void -FEEvaluationBase -::read_dof_values (const std::vector &src, - const unsigned int first_index) -{ - AssertIndexRange (first_index, src.size()); - Assert (n_fe_components == 1, ExcNotImplemented()); - Assert ((n_fe_components == 1 ? - (first_index+n_components <= src.size()) : true), - ExcIndexRange (first_index + n_components_, 0, src.size())); - - const VectorType *src_data [n_components]; - for (unsigned int comp=0; comp(src[comp+first_index]); + src_data[d] = internal::BlockVectorSelector::value>::get_vector_component(const_cast(src), d+first_index); internal::VectorReader reader; read_write_operation (reader, src_data); @@ -3205,14 +3129,15 @@ template inline void FEEvaluationBase -::read_dof_values_plain (const VectorType &src) +::read_dof_values_plain (const VectorType &src, + const unsigned int first_index) { // select between block vectors and non-block vectors. Note that the number // of components is checked in the internal data const typename internal::BlockVectorSelector::value>::BaseVectorType *src_data[n_components]; for (unsigned int d=0; d::value>::get_vector_component(const_cast(src), d); + src_data[d] = internal::BlockVectorSelector::value>::get_vector_component(const_cast(src), d+first_index); read_dof_values_plain (src_data); } @@ -3224,49 +3149,8 @@ template inline void FEEvaluationBase -::read_dof_values_plain (const std::vector &src, - const unsigned int first_index) -{ - AssertIndexRange (first_index, src.size()); - Assert (n_fe_components == 1, ExcNotImplemented()); - Assert ((n_fe_components == 1 ? - (first_index+n_components <= src.size()) : true), - ExcIndexRange (first_index + n_components_, 0, src.size())); - const VectorType *src_data [n_components]; - for (unsigned int comp=0; comp -template -inline -void -FEEvaluationBase -::read_dof_values_plain (const std::vector &src, - const unsigned int first_index) -{ - AssertIndexRange (first_index, src.size()); - Assert (n_fe_components == 1, ExcNotImplemented()); - Assert ((n_fe_components == 1 ? - (first_index+n_components <= src.size()) : true), - ExcIndexRange (first_index + n_components_, 0, src.size())); - const VectorType *src_data [n_components]; - for (unsigned int comp=0; comp -template -inline -void -FEEvaluationBase -::distribute_local_to_global (VectorType &dst) const +::distribute_local_to_global (VectorType &dst, + const unsigned int first_index) const { Assert (dof_values_initialized==true, internal::ExcAccessToUninitializedField()); @@ -3276,33 +3160,7 @@ FEEvaluationBase typename internal::BlockVectorSelector::value>::BaseVectorType *dst_data[n_components]; for (unsigned int d=0; d::value>::get_vector_component(dst, d); - - internal::VectorDistributorLocalToGlobal distributor; - read_write_operation (distributor, dst_data); -} - - - -template -template -inline -void -FEEvaluationBase -::distribute_local_to_global (std::vector &dst, - const unsigned int first_index) const -{ - AssertIndexRange (first_index, dst.size()); - Assert (n_fe_components == 1, ExcNotImplemented()); - Assert ((n_fe_components == 1 ? - (first_index+n_components <= dst.size()) : true), - ExcIndexRange (first_index + n_components_, 0, dst.size())); - Assert (dof_values_initialized==true, - internal::ExcAccessToUninitializedField()); - - VectorType *dst_data [n_components]; - for (unsigned int comp=0; comp::value>::get_vector_component(dst, d+first_index); internal::VectorDistributorLocalToGlobal distributor; read_write_operation (distributor, dst_data); @@ -3315,33 +3173,8 @@ template inline void FEEvaluationBase -::distribute_local_to_global (std::vector &dst, - const unsigned int first_index) const -{ - AssertIndexRange (first_index, dst.size()); - Assert (n_fe_components == 1, ExcNotImplemented()); - Assert ((n_fe_components == 1 ? - (first_index+n_components <= dst.size()) : true), - ExcIndexRange (first_index + n_components_, 0, dst.size())); - Assert (dof_values_initialized==true, - internal::ExcAccessToUninitializedField()); - - VectorType *dst_data [n_components]; - for (unsigned int comp=0; comp distributor; - read_write_operation (distributor, dst_data); -} - - - -template -template -inline -void -FEEvaluationBase -::set_dof_values (VectorType &dst) const +::set_dof_values (VectorType &dst, + const unsigned int first_index) const { Assert (dof_values_initialized==true, internal::ExcAccessToUninitializedField()); @@ -3351,61 +3184,7 @@ FEEvaluationBase typename internal::BlockVectorSelector::value>::BaseVectorType *dst_data[n_components]; for (unsigned int d=0; d::value>::get_vector_component(dst, d); - - internal::VectorSetter setter; - read_write_operation (setter, dst_data); -} - - - -template -template -inline -void -FEEvaluationBase -::set_dof_values (std::vector &dst, - const unsigned int first_index) const -{ - AssertIndexRange (first_index, dst.size()); - Assert (n_fe_components == 1, ExcNotImplemented()); - Assert ((n_fe_components == 1 ? - (first_index+n_components <= dst.size()) : true), - ExcIndexRange (first_index + n_components_, 0, dst.size())); - - Assert (dof_values_initialized==true, - internal::ExcAccessToUninitializedField()); - - VectorType *dst_data [n_components]; - for (unsigned int comp=0; comp setter; - read_write_operation (setter, dst_data); -} - - - -template -template -inline -void -FEEvaluationBase -::set_dof_values (std::vector &dst, - const unsigned int first_index) const -{ - AssertIndexRange (first_index, dst.size()); - Assert (n_fe_components == 1, ExcNotImplemented()); - Assert ((n_fe_components == 1 ? - (first_index+n_components <= dst.size()) : true), - ExcIndexRange (first_index + n_components_, 0, dst.size())); - - Assert (dof_values_initialized==true, - internal::ExcAccessToUninitializedField()); - - VectorType *dst_data [n_components]; - for (unsigned int comp=0; comp::value>::get_vector_component(dst, d+first_index); internal::VectorSetter setter; read_write_operation (setter, dst_data); -- 2.39.5