FEEval: do not cast const away during read_dof_values()

dealii · Jan 14, 2022 · 4a44987 · 4a44987
1 parent 3447d0a
commit 4a44987
Show file tree

Hide file tree

Showing 3 changed files with 80 additions and 23 deletions.
diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h
@@ -3270,6 +3270,22 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
 
 namespace internal
 {
+  template <typename VectorType, bool>
+  struct ConstBlockVectorSelector
+  {};
+
+  template <typename VectorType>
+  struct ConstBlockVectorSelector<VectorType, true>
+  {
+    using BaseVectorType = const typename VectorType::BlockType;
+  };
+
+  template <typename VectorType>
+  struct ConstBlockVectorSelector<VectorType, false>
+  {
+    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
@@ -3280,7 +3296,9 @@ namespace internal
   template <typename VectorType>
   struct BlockVectorSelector<VectorType, true>
   {
-    using BaseVectorType = typename VectorType::BlockType;
+    using BaseVectorType = typename ConstBlockVectorSelector<
+      VectorType,
+      std::is_const<VectorType>::value>::BaseVectorType;
 
     static BaseVectorType *
     get_vector_component(VectorType &vec, const unsigned int component)
@@ -3329,6 +3347,20 @@ namespace internal
     }
   };
 
+  template <typename VectorType>
+  struct BlockVectorSelector<const std::vector<VectorType>, false>
+  {
+    using BaseVectorType = const VectorType;
+
+    static const BaseVectorType *
+    get_vector_component(const std::vector<VectorType> &vec,
+                         const unsigned int             component)
+    {
+      AssertIndexRange(component, vec.size());
+      return &vec[component];
+    }
+  };
+
   template <typename VectorType>
   struct BlockVectorSelector<std::vector<VectorType *>, false>
   {
@@ -3342,6 +3374,20 @@ namespace internal
       return vec[component];
     }
   };
+
+  template <typename VectorType>
+  struct BlockVectorSelector<const std::vector<VectorType *>, false>
+  {
+    using BaseVectorType = const VectorType;
+
+    static const BaseVectorType *
+    get_vector_component(const std::vector<VectorType *> &vec,
+                         const unsigned int               component)
+    {
+      AssertIndexRange(component, vec.size());
+      return vec[component];
+    }
+  };
 } // namespace internal
 
 
@@ -4183,15 +4229,13 @@ namespace internal
   template <int n_components, typename VectorType>
   std::pair<
     std::array<typename internal::BlockVectorSelector<
-                 typename std::remove_const<VectorType>::type,
-                 IsBlockVector<typename std::remove_const<VectorType>::type>::
-                   value>::BaseVectorType *,
+                 VectorType,
+                 IsBlockVector<VectorType>::value>::BaseVectorType *,
                n_components>,
     std::array<
       const std::vector<ArrayView<const typename internal::BlockVectorSelector<
-        typename std::remove_const<VectorType>::type,
-        IsBlockVector<typename std::remove_const<VectorType>::type>::value>::
-                                    BaseVectorType::value_type>> *,
+        VectorType,
+        IsBlockVector<VectorType>::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<typename internal::BlockVectorSelector<
-                   typename std::remove_const<VectorType>::type,
-                   IsBlockVector<typename std::remove_const<VectorType>::type>::
-                     value>::BaseVectorType *,
+                   VectorType,
+                   IsBlockVector<VectorType>::value>::BaseVectorType *,
                  n_components>,
       std::array<
         const std::vector<
           ArrayView<const typename internal::BlockVectorSelector<
-            typename std::remove_const<VectorType>::type,
-            IsBlockVector<typename std::remove_const<VectorType>::type>::
-              value>::BaseVectorType::value_type>> *,
+            VectorType,
+            IsBlockVector<VectorType>::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<VectorType>::type,
-        IsBlockVector<typename std::remove_const<VectorType>::type>::value>::
-        get_vector_component(
-          const_cast<typename std::remove_const<VectorType>::type &>(src),
-          d + first_index);
+        VectorType,
+        IsBlockVector<VectorType>::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<typename internal::BlockVectorSelector<
+          typename std::remove_const<VectorType>::type,
+          IsBlockVector<typename std::remove_const<VectorType>::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
@@ -199,7 +199,8 @@ namespace internal
   {
     static const bool value =
       has_begin<T>::value &&
-      (has_local_element<T>::value || is_serial_vector<T>::value) &&
+      (has_local_element<T>::value ||
+       is_serial_vector<typename std::remove_const<T>::type>::value) &&
       std::is_same<typename T::value_type, Number>::value;
   };
 

diff --git a/include/deal.II/matrix_free/vector_access_internal.h b/include/deal.II/matrix_free/vector_access_internal.h
@@ -247,6 +247,12 @@ namespace internal
                             VectorizedArrayType *dof_values,
                             std::integral_constant<bool, true>) const
     {
+#ifdef DEBUG
+      for (unsigned int i = 0; i < dofs_per_cell; ++i)
+        for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
+          vector_access(vec, dof_index + v + i * VectorizedArrayType::size());
+#endif
+
       const Number *vec_ptr = vec.begin() + dof_index;
       for (unsigned int i = 0; i < dofs_per_cell;
            ++i, vec_ptr += VectorizedArrayType::size())
@@ -340,6 +346,11 @@ namespace internal
                        VectorizedArrayType &res,
                        std::integral_constant<bool, true>) const
     {
+#ifdef DEBUG
+      for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
+        vector_access(vec, indices[v] + constant_offset);
+#endif
+
       res.gather(vec.begin() + constant_offset, indices);
     }