MatrixFree: pass SM vectors to functions

diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h
index a27a0e3a9d615c1f1ec5d0c0c93340bb4c01ae6a..ce85c837bab7c312340e725831d1fa589ef3ea84 100644 (file)
--- a/include/deal.II/matrix_free/evaluation_kernels.h
+++ b/include/deal.II/matrix_free/evaluation_kernels.h
@@ -3186,9 +3186,10 @@ namespace internal
   {
     template <int fe_degree, int n_q_points_1d>
     static bool
-    run(const unsigned int n_components,
-        const unsigned int n_face_orientations,
-        const Number2 *    src_ptr,
+    run(const unsigned int                          n_components,
+        const unsigned int                          n_face_orientations,
+        const Number2 *                             src_ptr,
+        const std::vector<ArrayView<const Number>> *sm_ptr,
         const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &data,
         const MatrixFreeFunctions::DoFInfo &                       dof_info,
         VectorizedArrayType *                                      values_quad,
@@ -3211,9 +3212,12 @@ namespace internal
           return false;
         }
 
+      (void)sm_ptr;
+
       Processor<fe_degree, n_q_points_1d> p(n_components,
                                             false,
                                             src_ptr,
+                                            sm_ptr,
                                             data,
                                             dof_info,
                                             values_quad,
@@ -3248,9 +3252,10 @@ namespace internal
       using Number2_                  = const Number2;
 
       Processor(
-        const unsigned int n_components,
-        const bool         integrate,
-        const Number2 *    global_vector_ptr,
+        const unsigned int                          n_components,
+        const bool                                  integrate,
+        const Number2 *                             global_vector_ptr,
+        const std::vector<ArrayView<const Number>> *sm_ptr,
         const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &data,
         const MatrixFreeFunctions::DoFInfo &                       dof_info,
         VectorizedArrayType *                                      values_quad,
@@ -3270,6 +3275,7 @@ namespace internal
         : n_components(n_components)
         , integrate(integrate)
         , global_vector_ptr(global_vector_ptr)
+        , sm_ptr(sm_ptr)
         , data(data)
         , dof_info(dof_info)
         , values_quad(values_quad)
@@ -3402,9 +3408,10 @@ namespace internal
                              subface_index);
       }
 
-      const unsigned int n_components;
-      const bool         integrate;
-      const Number2 *    global_vector_ptr;
+      const unsigned int                          n_components;
+      const bool                                  integrate;
+      const Number2 *                             global_vector_ptr;
+      const std::vector<ArrayView<const Number>> *sm_ptr;
       const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &data;
       const MatrixFreeFunctions::DoFInfo &                       dof_info;
       VectorizedArrayType *                                      values_quad;
@@ -3432,9 +3439,10 @@ namespace internal
   {
     template <int fe_degree, int n_q_points_1d>
     static bool
-    run(const unsigned int n_components,
-        const unsigned int n_face_orientations,
-        Number2 *          dst_ptr,
+    run(const unsigned int                           n_components,
+        const unsigned int                           n_face_orientations,
+        Number2 *                                    dst_ptr,
+        const std::vector<ArrayView<const Number2>> *sm_ptr,
         const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &data,
         const MatrixFreeFunctions::DoFInfo &                       dof_info,
         VectorizedArrayType *                                      values_array,
@@ -3453,6 +3461,8 @@ namespace internal
                                       face_orientations,
         const Table<2, unsigned int> &orientation_map)
     {
+      (void)sm_ptr;
+
       if (dst_ptr == nullptr)
         {
           AssertDimension(n_face_orientations, 1);
@@ -3482,6 +3492,7 @@ namespace internal
                                             n_components,
                                             true,
                                             dst_ptr,
+                                            sm_ptr,
                                             data,
                                             dof_info,
                                             values_quad,
@@ -3517,10 +3528,11 @@ namespace internal
 
 
       Processor(
-        VectorizedArrayType *values_array,
-        const unsigned int   n_components,
-        const bool           integrate,
-        Number2 *            global_vector_ptr,
+        VectorizedArrayType *                       values_array,
+        const unsigned int                          n_components,
+        const bool                                  integrate,
+        Number2 *                                   global_vector_ptr,
+        const std::vector<ArrayView<const Number>> *sm_ptr,
         const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &data,
         const MatrixFreeFunctions::DoFInfo &                       dof_info,
         VectorizedArrayType *                                      values_quad,
@@ -3541,6 +3553,7 @@ namespace internal
         , n_components(n_components)
         , integrate(integrate)
         , global_vector_ptr(global_vector_ptr)
+        , sm_ptr(sm_ptr)
         , data(data)
         , dof_info(dof_info)
         , values_quad(values_quad)
@@ -3695,9 +3708,10 @@ namespace internal
       VectorizedArrayType *values_array;
 
 
-      const unsigned int n_components;
-      const bool         integrate;
-      Number2 *          global_vector_ptr;
+      const unsigned int                          n_components;
+      const bool                                  integrate;
+      Number2 *                                   global_vector_ptr;
+      const std::vector<ArrayView<const Number>> *sm_ptr;
       const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &data;
       const MatrixFreeFunctions::DoFInfo &                       dof_info;
       VectorizedArrayType *                                      values_quad;

diff --git a/include/deal.II/matrix_free/evaluation_template_factory.h b/include/deal.II/matrix_free/evaluation_template_factory.h
index dbb4e57a424ca7fc3856b1c8e72e5b49f86ac558..3052af2dffd40d673753dbcc538558fdc699a795 100644 (file)
--- a/include/deal.II/matrix_free/evaluation_template_factory.h
+++ b/include/deal.II/matrix_free/evaluation_template_factory.h
@@ -99,9 +99,10 @@ namespace internal
 
     static bool
     gather_evaluate(
-      const unsigned int n_components,
-      const std::size_t  n_face_orientations,
-      const Number *     src_ptr,
+      const unsigned int                          n_components,
+      const std::size_t                           n_face_orientations,
+      const Number *                              src_ptr,
+      const std::vector<ArrayView<const Number>> *sm_ptr,
       const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &data,
       const MatrixFreeFunctions::DoFInfo &                       dof_info,
       VectorizedArrayType *                                      values_quad,
@@ -121,9 +122,10 @@ namespace internal
 
     static bool
     integrate_scatter(
-      const unsigned int n_components,
-      const std::size_t  n_face_orientations,
-      Number *           dst_ptr,
+      const unsigned int                          n_components,
+      const std::size_t                           n_face_orientations,
+      Number *                                    dst_ptr,
+      const std::vector<ArrayView<const Number>> *sm_ptr,
       const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &data,
       const MatrixFreeFunctions::DoFInfo &                       dof_info,
       VectorizedArrayType *                                      values_array,

diff --git a/include/deal.II/matrix_free/evaluation_template_factory.templates.h b/include/deal.II/matrix_free/evaluation_template_factory.templates.h
index 168cf1c82c8067c1041601323224b87bf8fab13c..1a5c2548999eaa2c584e043631eda2a0faf6aafe 100644 (file)
--- a/include/deal.II/matrix_free/evaluation_template_factory.templates.h
+++ b/include/deal.II/matrix_free/evaluation_template_factory.templates.h
@@ -200,9 +200,10 @@ namespace internal
   template <int dim, typename Number, typename VectorizedArrayType>
   bool
   FEFaceEvaluationFactory<dim, Number, VectorizedArrayType>::gather_evaluate(
-    const unsigned int n_components,
-    const std::size_t  n_face_orientations,
-    const Number *     src_ptr,
+    const unsigned int                          n_components,
+    const std::size_t                           n_face_orientations,
+    const Number *                              src_ptr,
+    const std::vector<ArrayView<const Number>> *sm_ptr,
     const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &data,
     const MatrixFreeFunctions::DoFInfo &                       dof_info,
     VectorizedArrayType *                                      values_quad,
@@ -230,6 +231,7 @@ namespace internal
       n_components,
       n_face_orientations,
       src_ptr,
+      sm_ptr,
       data,
       dof_info,
       values_quad,
@@ -252,9 +254,10 @@ namespace internal
   template <int dim, typename Number, typename VectorizedArrayType>
   bool
   FEFaceEvaluationFactory<dim, Number, VectorizedArrayType>::integrate_scatter(
-    const unsigned int n_components,
-    const std::size_t  n_face_orientations,
-    Number *           dst_ptr,
+    const unsigned int                          n_components,
+    const std::size_t                           n_face_orientations,
+    Number *                                    dst_ptr,
+    const std::vector<ArrayView<const Number>> *sm_ptr,
     const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &data,
     const MatrixFreeFunctions::DoFInfo &                       dof_info,
     VectorizedArrayType *                                      values_array,
@@ -283,6 +286,7 @@ namespace internal
       n_components,
       n_face_orientations,
       dst_ptr,
+      sm_ptr,
       data,
       dof_info,
       values_array,

diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h
index a1be617d7bc7ebbb7e7fefd6debee41457d67d54..9d65dc211fda4ddfb0d46dee0ada0bb6e4cadbf6 100644 (file)
--- a/include/deal.II/matrix_free/fe_evaluation.h
+++ b/include/deal.II/matrix_free/fe_evaluation.h
@@ -1170,10 +1170,14 @@ protected:
    */
   template <typename VectorType, typename VectorOperation>
   void
-  read_write_operation(const VectorOperation &                        operation,
-                       const std::array<VectorType *, n_components_> &vectors,
-                       const std::bitset<VectorizedArrayType::size()> &mask,
-                       const bool apply_constraints = true) const;
+  read_write_operation(
+    const VectorOperation &                        operation,
+    const std::array<VectorType *, n_components_> &vectors,
+    const std::array<
+      const std::vector<ArrayView<const typename VectorType::value_type>> *,
+      n_components_> &                              vectors_sm,
+    const std::bitset<VectorizedArrayType::size()> &mask,
+    const bool apply_constraints = true) const;
 
   /**
    * A unified function to read from and write into vectors based on the given
@@ -1185,8 +1189,11 @@ protected:
   template <typename VectorType, typename VectorOperation>
   void
   read_write_operation_contiguous(
-    const VectorOperation &                         operation,
-    const std::array<VectorType *, n_components_> & vectors,
+    const VectorOperation &                        operation,
+    const std::array<VectorType *, n_components_> &vectors,
+    const std::array<
+      const std::vector<ArrayView<const typename VectorType::value_type>> *,
+      n_components_> &                              vectors_sm,
     const std::bitset<VectorizedArrayType::size()> &mask) const;
 
   /**
@@ -4181,10 +4188,14 @@ template <int dim,
 template <typename VectorType, typename VectorOperation>
 inline void
 FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
-  read_write_operation(const VectorOperation &                        operation,
-                       const std::array<VectorType *, n_components_> &src,
-                       const std::bitset<VectorizedArrayType::size()> &mask,
-                       const bool apply_constraints) const
+  read_write_operation(
+    const VectorOperation &                        operation,
+    const std::array<VectorType *, n_components_> &src,
+    const std::array<
+      const std::vector<ArrayView<const typename VectorType::value_type>> *,
+      n_components_> &                              src_sm,
+    const std::bitset<VectorizedArrayType::size()> &mask,
+    const bool                                      apply_constraints) const
 {
   // Case 1: No MatrixFree object given, simple case because we do not need to
   // process constraints and need not care about vectorization -> go to
@@ -4226,7 +4237,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
         [this->cell] >=
       internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::contiguous)
     {
-      read_write_operation_contiguous(operation, src, mask);
+      read_write_operation_contiguous(operation, src, src_sm, mask);
       return;
     }
 
@@ -4579,10 +4590,15 @@ template <typename VectorType, typename VectorOperation>
 inline void
 FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
   read_write_operation_contiguous(
-    const VectorOperation &                         operation,
-    const std::array<VectorType *, n_components_> & src,
+    const VectorOperation &                        operation,
+    const std::array<VectorType *, n_components_> &src,
+    const std::array<
+      const std::vector<ArrayView<const typename VectorType::value_type>> *,
+      n_components_> &                              vectors_sm,
     const std::bitset<VectorizedArrayType::size()> &mask) const
 {
+  (void)vectors_sm; // TODO: use it
+
   // This functions processes the functions read_dof_values,
   // distribute_local_to_global, and set_dof_values with the same code for
   // contiguous cell indices (DG case). The distinction between these three
@@ -4837,30 +4853,67 @@ namespace internal
     return nullptr;
   }
 
-  template <int n_components, typename VectorType>
-  std::array<typename internal::BlockVectorSelector<
-               typename std::remove_const<VectorType>::type,
-               IsBlockVector<typename std::remove_const<VectorType>::type>::
-                 value>::BaseVectorType *,
-             n_components>
-  get_vector_data(VectorType &src, const unsigned int first_index)
+  template <typename VectorType,
+            typename std::enable_if<has_shared_vector_data<VectorType>::value,
+                                    VectorType>::type * = nullptr>
+  const std::vector<ArrayView<const typename VectorType::value_type>> *
+  get_shared_vector_data(VectorType &vec /*other parameters?*/)
   {
-    // select between block vectors and non-block vectors. Note that the number
-    // of components is checked in the internal data
+    return vec.shared_vector_data();
+  }
+
+  template <typename VectorType,
+            typename std::enable_if<!has_shared_vector_data<VectorType>::value,
+                                    VectorType>::type * = nullptr>
+  const std::vector<ArrayView<const typename VectorType::value_type>> *
+  get_shared_vector_data(VectorType &)
+  {
+    return nullptr;
+  }
+
+  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 *,
-               n_components>
+               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>> *,
+      n_components>>
+  get_vector_data(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
+    std::pair<
+      std::array<typename internal::BlockVectorSelector<
+                   typename std::remove_const<VectorType>::type,
+                   IsBlockVector<typename std::remove_const<VectorType>::type>::
+                     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>> *,
+        n_components>>
       src_data;
+
     for (unsigned int d = 0; d < n_components; ++d)
-      src_data[d] = internal::BlockVectorSelector<
+      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);
 
+    for (unsigned int d = 0; d < n_components; ++d)
+      src_data.second[d] = get_shared_vector_data(*src_data.first[d]);
+
     return src_data;
   }
 } // namespace internal
@@ -4882,7 +4935,8 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
 
   internal::VectorReader<Number, VectorizedArrayType> reader;
   read_write_operation(reader,
-                       src_data,
+                       src_data.first,
+                       src_data.second,
                        std::bitset<VectorizedArrayType::size()>().flip(),
                        true);
 
@@ -4908,7 +4962,8 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
 
   internal::VectorReader<Number, VectorizedArrayType> reader;
   read_write_operation(reader,
-                       src_data,
+                       src_data.first,
+                       src_data.second,
                        std::bitset<VectorizedArrayType::size()>().flip(),
                        false);
 
@@ -4942,7 +4997,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
 
   internal::VectorDistributorLocalToGlobal<Number, VectorizedArrayType>
     distributor;
-  read_write_operation(distributor, dst_data, mask);
+  read_write_operation(distributor, dst_data.first, dst_data.second, mask);
 }
 
 
@@ -4968,7 +5023,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
     internal::get_vector_data<n_components_>(dst, first_index);
 
   internal::VectorSetter<Number, VectorizedArrayType> setter;
-  read_write_operation(setter, dst_data, mask);
+  read_write_operation(setter, dst_data.first, dst_data.second, mask);
 }
 
 
@@ -4995,7 +5050,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
     internal::get_vector_data<n_components_>(dst, first_index);
 
   internal::VectorSetter<Number, VectorizedArrayType> setter;
-  read_write_operation(setter, dst_data, mask, false);
+  read_write_operation(setter, dst_data.first, dst_data.second, mask, false);
 }
 
 
@@ -8472,6 +8527,10 @@ FEFaceEvaluation<dim,
       "Only EvaluationFlags::values and EvaluationFlags::gradients are supported."));
 
   const auto fu = [&]() {
+    // TODO
+    const auto shared_vector_data =
+      internal::get_shared_vector_data(input_vector);
+
     if (this->dof_access_index ==
           internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
         this->is_interior_face == false)
@@ -8488,6 +8547,7 @@ FEFaceEvaluation<dim,
           n_components,
           VectorizedArrayType::size(),
           internal::get_beginning<Number>(input_vector),
+          shared_vector_data,
           *this->data,
           *this->dof_info,
           this->begin_values(),
@@ -8528,6 +8588,7 @@ FEFaceEvaluation<dim,
               n_components,
               1,
               internal::get_beginning<Number>(input_vector),
+              shared_vector_data,
               *this->data,
               *this->dof_info,
               this->begin_values(),
@@ -8551,6 +8612,7 @@ FEFaceEvaluation<dim,
               gather_evaluate(n_components,
                               1,
                               internal::get_beginning<Number>(input_vector),
+                              shared_vector_data,
                               *this->data,
                               *this->dof_info,
                               this->begin_values(),
@@ -8637,6 +8699,9 @@ FEFaceEvaluation<dim,
           this->is_interior_face == false) == false,
          ExcNotImplemented());
 
+  // TODO
+  const auto shared_vector_data = internal::get_shared_vector_data(destination);
+
   // TODO: this copying should not be necessary once we have introduced
   // an internal-data structure
   std::array<unsigned int, VectorizedArrayType::size()> cells_            = {};
@@ -8657,6 +8722,7 @@ FEFaceEvaluation<dim,
             n_components,
             1,
             internal::get_beginning<Number>(destination),
+            shared_vector_data,
             *this->data,
             *this->dof_info,
             this->begin_dof_values(),
@@ -8688,6 +8754,7 @@ FEFaceEvaluation<dim,
             integrate_scatter(n_components,
                               1,
                               internal::get_beginning<Number>(destination),
+                              shared_vector_data,
                               *this->data,
                               *this->dof_info,
                               this->begin_dof_values(),

diff --git a/include/deal.II/matrix_free/type_traits.h b/include/deal.II/matrix_free/type_traits.h
index cbf49bd2bfb120debec4dee41bae8191d37bb698..b57c05049af3e9d13a9b8f9a8894b8dd9090b315 100644 (file)
--- a/include/deal.II/matrix_free/type_traits.h
+++ b/include/deal.II/matrix_free/type_traits.h
@@ -166,6 +166,29 @@ namespace internal
   const bool has_begin<T>::value;
 
 
+  // same as above to check
+  // ... T::shared_vector_data() const
+  template <typename T>
+  struct has_shared_vector_data
+  {
+  private:
+    static void
+    detect(...);
+
+    template <typename U>
+    static decltype(std::declval<U const>().shared_vector_data())
+    detect(const U &);
+
+  public:
+    static const bool value =
+      !std::is_same<void, decltype(detect(std::declval<T>()))>::value;
+  };
+
+  // We need to have a separate declaration for static const members
+  template <typename T>
+  const bool has_shared_vector_data<T>::value;
+
+
   // type trait for vector T and Number to see if
   // we can do vectorized load/save.
   // for VectorReader and VectorDistributorLocalToGlobal we assume that