From: David Wells <drwells@email.unc.edu>
Date: Thu, 1 Jun 2023 20:29:03 +0000 (-0400)
Subject: Make all vectors inherit from ReadVector.
X-Git-Tag: relicensing~820^2~8
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e690a7779e86cd545cc8e50f2a7a20a104426f06;p=dealii.git

Make all vectors inherit from ReadVector.

This gives a common interface to all vectors for use by FEValues.
---

diff --git a/include/deal.II/lac/block_vector_base.h b/include/deal.II/lac/block_vector_base.h
index edfbbf1b9e..33a9c266c1 100644
--- a/include/deal.II/lac/block_vector_base.h
+++ b/include/deal.II/lac/block_vector_base.h
@@ -26,6 +26,7 @@
 
 #include <deal.II/lac/block_indices.h>
 #include <deal.II/lac/exceptions.h>
+#include <deal.II/lac/read_vector.h>
 #include <deal.II/lac/vector.h>
 #include <deal.II/lac/vector_operation.h>
 
@@ -438,7 +439,8 @@ namespace internal
  * @ref GlossBlockLA "Block (linear algebra)"
  */
 template <class VectorType>
-class BlockVectorBase : public Subscriptor
+class BlockVectorBase : public Subscriptor,
+                        public ReadVector<typename VectorType::value_type>
 {
 public:
   /**
@@ -547,8 +549,8 @@ public:
    * Return dimension of the vector. This is the sum of the dimensions of all
    * components.
    */
-  std::size_t
-  size() const;
+  size_type
+  size() const override;
 
   /**
    * Return local dimension of the vector. This is the sum of the local
@@ -651,6 +653,10 @@ public:
   extract_subvector_to(const std::vector<size_type> &indices,
                        std::vector<OtherNumber> &    values) const;
 
+  virtual void
+  extract_subvector_to(const ArrayView<const types::global_dof_index> &indices,
+                       ArrayView<value_type> &entries) const override;
+
   /**
    * Instead of getting individual elements of a vector via operator(),
    * this function allows getting a whole set of elements at once. In
@@ -1439,7 +1445,7 @@ namespace internal
 
 
 template <class VectorType>
-inline std::size_t
+inline typename BlockVectorBase<VectorType>::size_type
 BlockVectorBase<VectorType>::size() const
 {
   return block_indices.total_size();
@@ -2146,6 +2152,22 @@ BlockVectorBase<VectorType>::extract_subvector_to(
 
 
 
+template <typename VectorType>
+inline void
+BlockVectorBase<VectorType>::extract_subvector_to(
+  const ArrayView<const types::global_dof_index> &indices,
+  ArrayView<value_type> &                         entries) const
+{
+  AssertDimension(indices.size(), entries.size());
+  for (unsigned int i = 0; i < indices.size(); ++i)
+    {
+      AssertIndexRange(indices[i], size());
+      entries[i] = (*this)[indices[i]];
+    }
+}
+
+
+
 template <typename VectorType>
 template <typename ForwardIterator, typename OutputIterator>
 inline void
diff --git a/include/deal.II/lac/la_parallel_vector.h b/include/deal.II/lac/la_parallel_vector.h
index 8777d0ce28..1ff552f623 100644
--- a/include/deal.II/lac/la_parallel_vector.h
+++ b/include/deal.II/lac/la_parallel_vector.h
@@ -27,6 +27,7 @@
 #include <deal.II/base/partitioner.h>
 #include <deal.II/base/subscriptor.h>
 
+#include <deal.II/lac/read_vector.h>
 #include <deal.II/lac/vector_operation.h>
 #include <deal.II/lac/vector_space_vector.h>
 #include <deal.II/lac/vector_type_traits.h>
@@ -248,6 +249,7 @@ namespace LinearAlgebra
      */
     template <typename Number, typename MemorySpace = MemorySpace::Host>
     class Vector : public ::dealii::LinearAlgebra::VectorSpaceVector<Number>,
+                   public ::dealii::ReadVector<Number>,
                    public Subscriptor
     {
     public:
@@ -1144,6 +1146,14 @@ namespace LinearAlgebra
       extract_subvector_to(const std::vector<size_type> &indices,
                            std::vector<OtherNumber> &    values) const;
 
+      /**
+       * Extract a range of elements all at once.
+       */
+      virtual void
+      extract_subvector_to(
+        const ArrayView<const types::global_dof_index> &indices,
+        ArrayView<Number> &elements) const override;
+
       /**
        * Instead of getting individual elements of a vector via operator(),
        * this function allows getting a whole set of elements at once. In
diff --git a/include/deal.II/lac/la_parallel_vector.templates.h b/include/deal.II/lac/la_parallel_vector.templates.h
index 848aa0446a..791326918a 100644
--- a/include/deal.II/lac/la_parallel_vector.templates.h
+++ b/include/deal.II/lac/la_parallel_vector.templates.h
@@ -1754,6 +1754,22 @@ namespace LinearAlgebra
 
 
 
+    template <typename Number, typename MemorySpaceType>
+    void
+    Vector<Number, MemorySpaceType>::extract_subvector_to(
+      const ArrayView<const types::global_dof_index> &indices,
+      ArrayView<Number> &                             elements) const
+    {
+      AssertDimension(indices.size(), elements.size());
+      for (unsigned int i = 0; i < indices.size(); ++i)
+        {
+          AssertIndexRange(indices[i], size());
+          elements[i] = (*this)[indices[i]];
+        }
+    }
+
+
+
     template <typename Number, typename MemorySpaceType>
     bool
     Vector<Number, MemorySpaceType>::all_zero() const
diff --git a/include/deal.II/lac/petsc_vector_base.h b/include/deal.II/lac/petsc_vector_base.h
index 8dfe3ebf33..3c768589d2 100644
--- a/include/deal.II/lac/petsc_vector_base.h
+++ b/include/deal.II/lac/petsc_vector_base.h
@@ -247,7 +247,7 @@ namespace PETScWrappers
    *
    * @ingroup PETScWrappers
    */
-  class VectorBase : public Subscriptor
+  class VectorBase : public ReadVector<PetscScalar>, public Subscriptor
   {
   public:
     /**
@@ -355,7 +355,7 @@ namespace PETScWrappers
      * Return the global dimension of the vector.
      */
     size_type
-    size() const;
+    size() const override;
 
     /**
      * Return the local dimension of the vector, i.e. the number of elements
@@ -494,6 +494,14 @@ namespace PETScWrappers
     extract_subvector_to(const std::vector<size_type> &indices,
                          std::vector<PetscScalar> &    values) const;
 
+    /**
+     * Extract a range of elements all at once.
+     */
+    virtual void
+    extract_subvector_to(
+      const ArrayView<const types::global_dof_index> &indices,
+      ArrayView<PetscScalar> &                        elements) const override;
+
     /**
      * Instead of getting individual elements of a vector via operator(),
      * this function allows getting a whole set of elements at once. In
@@ -1192,6 +1200,16 @@ namespace PETScWrappers
     extract_subvector_to(indices.begin(), indices.end(), values.begin());
   }
 
+  inline void
+  VectorBase::extract_subvector_to(
+    const ArrayView<const types::global_dof_index> &indices,
+    ArrayView<PetscScalar> &                        elements) const
+  {
+    AssertDimension(indices.size(), elements.size());
+    extract_subvector_to(indices.begin(), indices.end(), elements.begin());
+  }
+
+
   template <typename ForwardIterator, typename OutputIterator>
   inline void
   VectorBase::extract_subvector_to(const ForwardIterator indices_begin,
diff --git a/include/deal.II/lac/read_vector.h b/include/deal.II/lac/read_vector.h
new file mode 100644
index 0000000000..1a61c83097
--- /dev/null
+++ b/include/deal.II/lac/read_vector.h
@@ -0,0 +1,64 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2023 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii_lac_read_vector_h
+#define dealii_lac_read_vector_h
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/array_view.h>
+#include <deal.II/base/types.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+/**
+ * @addtogroup Vectors
+ * @{
+ */
+
+/**
+ * Base class for providing read-only access to vector elements.
+ *
+ * deal.II supports a large number of vector classes, including both its own
+ * serial and parallel vectors as well as vector classes from external
+ * libraries like PETSc and Trilinos. ReadVector is a common base class for
+ * all vector classes and defines a minimal interface for efficiently
+ * accessing vector elements.
+ */
+template <typename Number>
+class ReadVector
+{
+public:
+  using size_type = types::global_dof_index;
+
+  /**
+   * Return the size of the vector.
+   */
+  virtual size_type
+  size() const = 0;
+
+  /**
+   * Extract a subset of the vector specified by @p indices into the output
+   * array @p elements.
+   */
+  virtual void
+  extract_subvector_to(const ArrayView<const types::global_dof_index> &indices,
+                       ArrayView<Number> &elements) const = 0;
+};
+
+/** @} */
+
+DEAL_II_NAMESPACE_CLOSE
+#endif
diff --git a/include/deal.II/lac/read_write_vector.h b/include/deal.II/lac/read_write_vector.h
index c3464594cf..f25ddb3761 100644
--- a/include/deal.II/lac/read_write_vector.h
+++ b/include/deal.II/lac/read_write_vector.h
@@ -27,6 +27,7 @@
 #include <deal.II/base/template_constraints.h>
 #include <deal.II/base/types.h>
 
+#include <deal.II/lac/read_vector.h>
 #include <deal.II/lac/vector_operation.h>
 
 #include <cstdlib>
@@ -133,7 +134,7 @@ namespace LinearAlgebra
    * get the first index of the largest range.
    */
   template <typename Number>
-  class ReadWriteVector : public Subscriptor
+  class ReadWriteVector : public Subscriptor, public ReadVector<Number>
   {
   public:
     /**
@@ -537,7 +538,7 @@ namespace LinearAlgebra
      * number returned by the current function.
      */
     size_type
-    size() const;
+    size() const override;
 
     /**
      * This function returns the number of elements stored. It is smaller or
@@ -655,6 +656,14 @@ namespace LinearAlgebra
     extract_subvector_to(const std::vector<size_type> &indices,
                          std::vector<Number2> &        values) const;
 
+    /**
+     * Extract a range of elements all at once.
+     */
+    virtual void
+    extract_subvector_to(
+      const ArrayView<const types::global_dof_index> &indices,
+      ArrayView<Number> &                             entries) const override;
+
     /**
      * Instead of getting individual elements of a vector via operator(),
      * this function allows getting a whole set of elements at once. In
@@ -1080,6 +1089,19 @@ namespace LinearAlgebra
 
 
 
+  template <typename Number>
+  void
+  ReadWriteVector<Number>::extract_subvector_to(
+    const ArrayView<const types::global_dof_index> &indices,
+    ArrayView<Number> &                             entries) const
+  {
+    AssertDimension(indices.size(), entries.size());
+    for (unsigned int i = 0; i < indices.size(); ++i)
+      entries[i] = (*this)[indices[i]];
+  }
+
+
+
   template <typename Number>
   template <typename ForwardIterator, typename OutputIterator>
   inline void
diff --git a/include/deal.II/lac/trilinos_epetra_vector.h b/include/deal.II/lac/trilinos_epetra_vector.h
index dba315e093..b7a630548b 100644
--- a/include/deal.II/lac/trilinos_epetra_vector.h
+++ b/include/deal.II/lac/trilinos_epetra_vector.h
@@ -25,6 +25,7 @@
 #  include <deal.II/base/mpi_stub.h>
 #  include <deal.II/base/subscriptor.h>
 
+#  include <deal.II/lac/read_vector.h>
 #  include <deal.II/lac/trilinos_epetra_communication_pattern.h>
 #  include <deal.II/lac/vector_operation.h>
 #  include <deal.II/lac/vector_space_vector.h>
@@ -61,9 +62,14 @@ namespace LinearAlgebra
      * @ingroup TrilinosWrappers
      * @ingroup Vectors
      */
-    class Vector : public VectorSpaceVector<double>, public Subscriptor
+    class Vector : public VectorSpaceVector<double>,
+                   public ReadVector<double>,
+                   public Subscriptor
     {
     public:
+      using value_type = double;
+      using size_type  = types::global_dof_index;
+
       /**
        * Constructor. Create a vector of dimension zero.
        */
@@ -103,6 +109,14 @@ namespace LinearAlgebra
       reinit(const VectorSpaceVector<double> &V,
              const bool omit_zeroing_entries = false) override;
 
+      /**
+       * Extract a range of elements all at once.
+       */
+      virtual void
+      extract_subvector_to(
+        const ArrayView<const types::global_dof_index> &indices,
+        ArrayView<double> &elements) const override;
+
       /**
        * Copy function. This function takes a Vector and copies all the
        * elements. The Vector will have the same parallel distribution as @p
diff --git a/include/deal.II/lac/trilinos_tpetra_vector.h b/include/deal.II/lac/trilinos_tpetra_vector.h
index 00db7ae697..3ad3cd01f2 100644
--- a/include/deal.II/lac/trilinos_tpetra_vector.h
+++ b/include/deal.II/lac/trilinos_tpetra_vector.h
@@ -25,6 +25,7 @@
 #  include <deal.II/base/mpi_stub.h>
 #  include <deal.II/base/subscriptor.h>
 
+#  include <deal.II/lac/read_vector.h>
 #  include <deal.II/lac/trilinos_tpetra_communication_pattern.h>
 #  include <deal.II/lac/vector_operation.h>
 #  include <deal.II/lac/vector_space_vector.h>
@@ -113,7 +114,9 @@ namespace LinearAlgebra
      * @ingroup Vectors
      */
     template <typename Number>
-    class Vector : public VectorSpaceVector<Number>, public Subscriptor
+    class Vector : public VectorSpaceVector<Number>,
+                   public ReadVector<Number>,
+                   public Subscriptor
     {
     public:
       using value_type = Number;
@@ -159,6 +162,14 @@ namespace LinearAlgebra
       reinit(const VectorSpaceVector<Number> &V,
              const bool omit_zeroing_entries = false) override;
 
+      /**
+       * Extract a range of elements all at once.
+       */
+      virtual void
+      extract_subvector_to(
+        const ArrayView<const types::global_dof_index> &indices,
+        ArrayView<Number> &elements) const override;
+
       /**
        * Copy function. This function takes a Vector and copies all the
        * elements. The Vector will have the same parallel distribution as @p
diff --git a/include/deal.II/lac/trilinos_tpetra_vector.templates.h b/include/deal.II/lac/trilinos_tpetra_vector.templates.h
index 3636ba7b3d..92b9908818 100644
--- a/include/deal.II/lac/trilinos_tpetra_vector.templates.h
+++ b/include/deal.II/lac/trilinos_tpetra_vector.templates.h
@@ -118,6 +118,36 @@ namespace LinearAlgebra
 
 
 
+    template <typename Number>
+    void
+    Vector::extract_subvector_to(
+      const ArrayView<const types::global_dof_index> &indices,
+      ArrayView<Number> &                             elements) const
+    {
+      AssertDimension(indices.size(), elements.size());
+      const auto &vector = trilinos_vector();
+      const auto &map    = vector.Map();
+
+#  if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
+      auto vector_2d = vector.template getLocalView<Kokkos::HostSpace>(
+        Tpetra::Access::ReadOnly);
+#  else
+      vector.template sync<Kokkos::HostSpace>();
+      auto vector_2d = vector.template getLocalView<Kokkos::HostSpace>();
+#  endif
+      auto vector_1d = Kokkos::subview(vector_2d, Kokkos::ALL(), 0);
+
+      for (unsigned int i = 0; i < indices.size(); ++i)
+        {
+          AssertIndexRange(indices[i], size());
+          const trilinos_i = map->getLocalElement(
+            static_cast<TrilinosWrappers::types::int_type>(indices[i]));
+          elements[i] = vector_1d(trilinos_i);
+        }
+    }
+
+
+
     template <typename Number>
     Vector<Number> &
     Vector<Number>::operator=(const Vector<Number> &V)
diff --git a/include/deal.II/lac/trilinos_vector.h b/include/deal.II/lac/trilinos_vector.h
index 1d77616a1e..a27f0857ed 100644
--- a/include/deal.II/lac/trilinos_vector.h
+++ b/include/deal.II/lac/trilinos_vector.h
@@ -26,6 +26,7 @@
 #  include <deal.II/base/subscriptor.h>
 
 #  include <deal.II/lac/exceptions.h>
+#  include <deal.II/lac/read_vector.h>
 #  include <deal.II/lac/vector.h>
 #  include <deal.II/lac/vector_operation.h>
 #  include <deal.II/lac/vector_type_traits.h>
@@ -394,7 +395,7 @@ namespace TrilinosWrappers
      * @ingroup TrilinosWrappers
      * @ingroup Vectors
      */
-    class Vector : public Subscriptor
+    class Vector : public Subscriptor, public ReadVector<TrilinosScalar>
     {
     public:
       /**
@@ -760,7 +761,7 @@ namespace TrilinosWrappers
        * Return the global dimension of the vector.
        */
       size_type
-      size() const;
+      size() const override;
 
       /**
        * Return the local dimension of the vector, i.e. the number of elements
@@ -1012,6 +1013,13 @@ namespace TrilinosWrappers
       extract_subvector_to(const std::vector<size_type> &indices,
                            std::vector<TrilinosScalar> & values) const;
 
+      /**
+       * Extract a range of elements all at once.
+       */
+      virtual void
+      extract_subvector_to(const ArrayView<const size_type> &indices,
+                           ArrayView<TrilinosScalar> &elements) const override;
+
       /**
        * Instead of getting individual elements of a vector via operator(),
        * this function allows getting a whole set of elements at once. In
@@ -1560,6 +1568,20 @@ namespace TrilinosWrappers
 
 
 
+    inline void
+    Vector::extract_subvector_to(const ArrayView<const size_type> &indices,
+                                 ArrayView<TrilinosScalar> &elements) const
+    {
+      AssertDimension(indices.size(), elements.size());
+      for (unsigned int i = 0; i < indices.size(); ++i)
+        {
+          AssertIndexRange(indices[i], size());
+          elements[i] = (*this)[indices[i]];
+        }
+    }
+
+
+
     template <typename ForwardIterator, typename OutputIterator>
     inline void
     Vector::extract_subvector_to(ForwardIterator       indices_begin,
diff --git a/include/deal.II/lac/vector.h b/include/deal.II/lac/vector.h
index dcd98babf4..adb6a1deb4 100644
--- a/include/deal.II/lac/vector.h
+++ b/include/deal.II/lac/vector.h
@@ -25,6 +25,7 @@
 #include <deal.II/base/numbers.h>
 #include <deal.II/base/subscriptor.h>
 
+#include <deal.II/lac/read_vector.h>
 #include <deal.II/lac/vector_operation.h>
 #include <deal.II/lac/vector_type_traits.h>
 
@@ -105,7 +106,7 @@ namespace parallel
  * in the manual).
  */
 template <typename Number>
-class Vector : public Subscriptor
+class Vector : public Subscriptor, public ReadVector<Number>
 {
 public:
   /**
@@ -662,6 +663,13 @@ public:
   extract_subvector_to(const std::vector<size_type> &indices,
                        std::vector<OtherNumber> &    values) const;
 
+  /**
+   * Extract a range of elements all at once.
+   */
+  virtual void
+  extract_subvector_to(const ArrayView<const types::global_dof_index> &indices,
+                       ArrayView<Number> &elements) const override;
+
   /**
    * Instead of getting individual elements of a vector via operator(),
    * this function allows getting a whole set of elements at once. In
@@ -958,8 +966,8 @@ public:
   /**
    * Return dimension of the vector.
    */
-  size_type
-  size() const;
+  virtual size_type
+  size() const override;
 
   /**
    * Return local dimension of the vector. Since this vector does not support
diff --git a/include/deal.II/lac/vector.templates.h b/include/deal.II/lac/vector.templates.h
index 79920ed76e..52cba5a014 100644
--- a/include/deal.II/lac/vector.templates.h
+++ b/include/deal.II/lac/vector.templates.h
@@ -557,6 +557,22 @@ Vector<Number>::add_and_dot(const Number          a,
 
 
 
+template <typename Number>
+void
+Vector<Number>::extract_subvector_to(
+  const ArrayView<const types::global_dof_index> &indices,
+  ArrayView<Number> &                             elements) const
+{
+  AssertDimension(indices.size(), elements.size());
+  for (unsigned int i = 0; i < indices.size(); ++i)
+    {
+      AssertIndexRange(indices[i], size());
+      elements[i] = (*this)[indices[i]];
+    }
+}
+
+
+
 template <typename Number>
 Vector<Number> &
 Vector<Number>::operator+=(const Vector<Number> &v)
diff --git a/source/lac/trilinos_epetra_vector.cc b/source/lac/trilinos_epetra_vector.cc
index 599a0566d6..9a69a287c8 100644
--- a/source/lac/trilinos_epetra_vector.cc
+++ b/source/lac/trilinos_epetra_vector.cc
@@ -99,6 +99,26 @@ namespace LinearAlgebra
 
 
 
+    void
+    Vector::extract_subvector_to(
+      const ArrayView<const types::global_dof_index> &indices,
+      ArrayView<double> &                             elements) const
+    {
+      AssertDimension(indices.size(), elements.size());
+      const auto &vector = trilinos_vector();
+      const auto &map    = vector.Map();
+
+      for (unsigned int i = 0; i < indices.size(); ++i)
+        {
+          AssertIndexRange(indices[i], size());
+          const auto trilinos_i =
+            map.LID(static_cast<TrilinosWrappers::types::int_type>(indices[i]));
+          elements[i] = vector[0][trilinos_i];
+        }
+    }
+
+
+
     Vector &
     Vector::operator=(const Vector &V)
     {