]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Make all vectors inherit from ReadVector.
authorDavid Wells <drwells@email.unc.edu>
Thu, 1 Jun 2023 20:29:03 +0000 (16:29 -0400)
committerDavid Wells <drwells@email.unc.edu>
Sat, 1 Jul 2023 14:19:08 +0000 (10:19 -0400)
This gives a common interface to all vectors for use by FEValues.

13 files changed:
include/deal.II/lac/block_vector_base.h
include/deal.II/lac/la_parallel_vector.h
include/deal.II/lac/la_parallel_vector.templates.h
include/deal.II/lac/petsc_vector_base.h
include/deal.II/lac/read_vector.h [new file with mode: 0644]
include/deal.II/lac/read_write_vector.h
include/deal.II/lac/trilinos_epetra_vector.h
include/deal.II/lac/trilinos_tpetra_vector.h
include/deal.II/lac/trilinos_tpetra_vector.templates.h
include/deal.II/lac/trilinos_vector.h
include/deal.II/lac/vector.h
include/deal.II/lac/vector.templates.h
source/lac/trilinos_epetra_vector.cc

index edfbbf1b9e344ee1e841763743b1981943416866..33a9c266c13f190552380e967301459df3920b61 100644 (file)
@@ -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
index 8777d0ce2822850723a4c6d80bc0a04e11fb92d4..1ff552f623c6ed4b4f8e2275cc506458adf269bb 100644 (file)
@@ -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
index 848aa0446a8f67e7fbee49c9691c2b4ea951e935..791326918a6b91edb37ef2c08988677cca86eb22 100644 (file)
@@ -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
index 8dfe3ebf33715bc45b5b56e7133d2dcf6c6a04c8..3c768589d2399f9e2b608d42da65c83d33272758 100644 (file)
@@ -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 (file)
index 0000000..1a61c83
--- /dev/null
@@ -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
index c3464594cfc719b25150d9eb95f37bf0dbdeaedd..f25ddb3761dec57c5ba42172aa6cb931b884f627 100644 (file)
@@ -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
index dba315e093fbf0849bec6270cd51f2cd7b62e4df..b7a630548bfd27056179b681928be37df7e4cd75 100644 (file)
@@ -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
index 00db7ae697dea330ef048a64f19fbca1cbbb836f..3ad3cd01f2e85cdec57915f39ca63f580ca2db04 100644 (file)
@@ -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
index 3636ba7b3d7aa2f41334d5e8cee41be772fa7aa7..92b990881861684c023b2f11d5077a006fee49cc 100644 (file)
@@ -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)
index 1d77616a1eaf2db00f752b670a39d630dbf026d7..a27f0857edc5a3712fad8b7566b3d43586e487bf 100644 (file)
@@ -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,
index dcd98babf4e2bb2c088feb45f940cafacfe23ea4..adb6a1deb4448abc68a46c8886e3bd8aea79a4ec 100644 (file)
@@ -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
index 79920ed76e060a2b76a801b38265c2654e4e685f..52cba5a014cb5fb28e4a6b9c456770ca1a707f31 100644 (file)
@@ -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)
index 599a0566d60a65cceafd9f858341a3f9351f0dc1..9a69a287c8dcf09ed31c1830d03030ba86bdaa8f 100644 (file)
@@ -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)
     {

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.