Enable matrix-free loops with ArrayView arguments

diff --git a/include/deal.II/matrix_free/type_traits.h b/include/deal.II/matrix_free/type_traits.h
index 9ed8fadf558ec22627030e7750cfa829137df619..7b7ebb893660517b4ed39a8aae2d5c0cf59a79dc 100644 (file)
--- a/include/deal.II/matrix_free/type_traits.h
+++ b/include/deal.II/matrix_free/type_traits.h
@@ -28,13 +28,14 @@
 DEAL_II_NAMESPACE_OPEN
 
 #ifndef DOXYGEN
+
 namespace internal
 {
   //
   // type traits for FEEvaluation
   //
 
-  // a helper type-trait that leverage SFINAE to figure out if type T has
+  // a helper type-trait leveraging SFINAE to figure out if type T has
   // ... T::local_element() const
   template <typename T>
   using local_element_t = decltype(std::declval<const T>().local_element(0));
@@ -44,7 +45,7 @@ namespace internal
 
 
 
-  // a helper type-trait that leverage SFINAE to figure out if type T has
+  // a helper type-trait leveraging SFINAE to figure out if type T has
   // void T::add_local_element(const uint, const typename T::value_type)
   template <typename T>
   using add_local_element_t =
@@ -56,7 +57,7 @@ namespace internal
 
 
 
-  // a helper type-trait that leverage SFINAE to figure out if type T has
+  // a helper type-trait leveraging SFINAE to figure out if type T has
   // void T::set_local_element(const uint, const typename T::value_type)
   template <typename T>
   using set_local_element_t =
@@ -67,7 +68,6 @@ namespace internal
     is_supported_operation<set_local_element_t, T>;
 
 
-
   // same as above to check
   // bool T::partitioners_are_compatible(const Utilities::MPI::Partitioner &)
   // const
@@ -103,27 +103,6 @@ namespace internal
     is_supported_operation<shared_vector_data_t, T>;
 
 
-
-  // type trait for vector T and Number to see if
-  // we can do vectorized load/save.
-  // for VectorReader and VectorDistributorLocalToGlobal we assume that
-  // if both begin() and local_element()
-  // exist, then begin() + offset == local_element(offset)
-  template <typename T, typename Number>
-  struct is_vectorizable
-  {
-    static const bool value =
-      has_begin<T> &&
-      (has_local_element<T> ||
-       is_serial_vector<std::remove_const_t<T>>::value) &&
-      std::is_same_v<typename T::value_type, Number>;
-  };
-
-  // We need to have a separate declaration for static const members
-  template <typename T, typename Number>
-  const bool is_vectorizable<T, Number>::value;
-
-
   //
   // type-traits for Matrix-Free
   //
@@ -176,17 +155,6 @@ namespace internal
   constexpr bool has_communication_block_size =
     is_supported_operation<communication_block_size_t, T>;
 
-
-
-  // type trait for vector T to see if
-  // we need to do any data exchange for this vector type at all.
-  // is_serial_vector<> would have been enough, but in some circumstances
-  // (like calculation of diagonals for matrix-free operators)
-  // a dummy InVector == unsigned int is provided.
-  // Thus we have to treat this case as well.
-  template <class T, class IsSerialVectorNotSpecialized = void>
-  using not_parallel_vector_t = std::bool_constant<is_serial_vector<T>::value>;
-
   /**
    * A predicate stating whether something is a vector type. We test this
    * by seeing whether the `is_serial_vector` type is declared for the
@@ -212,6 +180,54 @@ namespace internal
   constexpr bool is_not_parallel_vector =
     (is_supported_operation<is_vector_type, VectorType> == false) ||
     (is_supported_operation<is_serial_vector_type, VectorType> == true);
+
+
+
+  /**
+   * A type trait for vector T indicating serial vectors with access through
+   * operator[], which are the deal.II serial vectors and ArrayView<Number>
+   */
+  template <typename VectorType>
+  struct is_serial_vector_or_array
+  {
+    static const bool value =
+      is_supported_operation<is_serial_vector_type, VectorType>;
+  };
+
+  /**
+   * A type trait for vector T indicating serial vectors with access through
+   * operator[], which are the deal.II serial vectors and ArrayView<Number>
+   */
+  template <typename Number>
+  struct is_serial_vector_or_array<dealii::ArrayView<Number>>
+  {
+    static const bool value = true;
+  };
+
+  // We need to have a separate declaration for static const members
+  template <typename VectorType>
+  const bool is_serial_vector_or_array<VectorType>::value;
+
+
+
+  // type trait for vector T and Number to see if
+  // we can do vectorized load/save.
+  // for VectorReader and VectorDistributorLocalToGlobal we assume that
+  // if both begin() and local_element()
+  // exist, then begin() + offset == local_element(offset)
+  template <typename T, typename Number>
+  struct is_vectorizable
+  {
+    static const bool value =
+      has_begin<T> &&
+      (has_local_element<T> || is_serial_vector_or_array<T>::value) &&
+      std::is_same_v<typename T::value_type, Number>;
+  };
+
+  // We need to have a separate declaration for static const members
+  template <typename T, typename Number>
+  const bool is_vectorizable<T, Number>::value;
+
 } // namespace internal
 #endif
 

diff --git a/include/deal.II/matrix_free/vector_access_internal.h b/include/deal.II/matrix_free/vector_access_internal.h
index b371c650e32e57c31ba2c082fa8ca291306b337c..b156b25a2bc2983cf06b11932835b0d1394ccedb 100644 (file)
--- a/include/deal.II/matrix_free/vector_access_internal.h
+++ b/include/deal.II/matrix_free/vector_access_internal.h
@@ -35,30 +35,28 @@ namespace internal
 {
   // below we use type-traits from matrix-free/type_traits.h
 
-  // access to generic const vectors that have operator ().
-  // FIXME: this is wrong for Trilinos/PETSc MPI vectors
-  // where we should first do Partitioner::local_to_global()
-  template <
-    typename VectorType,
-    std::enable_if_t<!has_local_element<VectorType>, VectorType> * = nullptr>
+
+
+  // access to serial const vectors that have operator[].
+  template <typename VectorType,
+            std::enable_if_t<is_serial_vector_or_array<VectorType>::value,
+                             VectorType> * = nullptr>
   inline typename VectorType::value_type
   vector_access(const VectorType &vec, const unsigned int entry)
   {
-    return vec(entry);
+    return vec[entry];
   }
 
 
 
-  // access to generic non-const vectors that have operator ().
-  // FIXME: this is wrong for Trilinos/PETSc MPI vectors
-  // where we should first do Partitioner::local_to_global()
-  template <
-    typename VectorType,
-    std::enable_if_t<!has_local_element<VectorType>, VectorType> * = nullptr>
+  // access to serial non-const vectors that have operator[].
+  template <typename VectorType,
+            std::enable_if_t<is_serial_vector_or_array<VectorType>::value,
+                             VectorType> * = nullptr>
   inline typename VectorType::value_type &
   vector_access(VectorType &vec, const unsigned int entry)
   {
-    return vec(entry);
+    return vec[entry];
   }
 
 
@@ -136,7 +134,7 @@ namespace internal
                            const types::global_dof_index          entry,
                            const typename VectorType::value_type &val)
   {
-    vec(entry) += val;
+    vec[entry] += val;
   }
 
 
@@ -483,7 +481,7 @@ namespace internal
                        const VectorType             &vec,
                        Number                       &res) const
     {
-      res = vec(index);
+      res = vec[index];
     }
 
 
@@ -953,7 +951,7 @@ namespace internal
                        VectorType                   &vec,
                        Number                       &res) const
     {
-      vec(index) = res;
+      vec[index] = res;
     }
 
 

diff --git a/tests/matrix_free/fe_evaluation_type_traits.cc b/tests/matrix_free/fe_evaluation_type_traits.cc
index 927804c12f031f71463ce8f3867d8b43b7f8c665..a217ab5a4042194dc9655c94838a3a4e59c2fda5 100644 (file)
--- a/tests/matrix_free/fe_evaluation_type_traits.cc
+++ b/tests/matrix_free/fe_evaluation_type_traits.cc
@@ -71,14 +71,14 @@ public:
   using value_type = Number;
 
   Number
-  operator()(const unsigned int) const
+  operator[](const unsigned int) const
   {
     deallog << "Dummy2::operator() const" << std::endl;
     return Number();
   }
 
   Number &
-  operator()(const unsigned int)
+  operator[](const unsigned int)
   {
     deallog << "Dummy2::operator()" << std::endl;
     return dummy;
@@ -112,7 +112,7 @@ main()
   // we expect local_element() to be called
   deallog << "internal::vector_access:" << std::endl;
   internal::vector_access(dummy, 0);
-  internal::vector_access(dummy2, 0);
+  // internal::vector_access(dummy2, 0);
 
 
   // now check has_partitioners_are_compatible:

diff --git a/tests/matrix_free/fe_evaluation_type_traits.with_trilinos=true.output b/tests/matrix_free/fe_evaluation_type_traits.with_trilinos=true.output
index b7e821fc996fcc8f6507c1abff7ea3df6fcfb9db..ed4d98ded68c8c1c4885e18a13e927dd970573f0 100644 (file)
--- a/tests/matrix_free/fe_evaluation_type_traits.with_trilinos=true.output
+++ b/tests/matrix_free/fe_evaluation_type_traits.with_trilinos=true.output
@@ -7,7 +7,6 @@ DEAL::Dummy2 = 0
 DEAL::Vector = 0
 DEAL::internal::vector_access:
 DEAL::Dummy::local_element()
-DEAL::Dummy2::operator()
 DEAL::has_partitioners_are_compatible:
 DEAL::LinearAlgebra::distributed::Vector = 1
 DEAL::TrilinosWrappers::MPI::Vector = 0