From: David Wells 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 #include +#include #include #include @@ -438,7 +439,8 @@ namespace internal * @ref GlossBlockLA "Block (linear algebra)" */ template -class BlockVectorBase : public Subscriptor +class BlockVectorBase : public Subscriptor, + public ReadVector { 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 &indices, std::vector & values) const; + virtual void + extract_subvector_to(const ArrayView &indices, + ArrayView &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 -inline std::size_t +inline typename BlockVectorBase::size_type BlockVectorBase::size() const { return block_indices.total_size(); @@ -2146,6 +2152,22 @@ BlockVectorBase::extract_subvector_to( +template +inline void +BlockVectorBase::extract_subvector_to( + const ArrayView &indices, + ArrayView & 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 template 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 #include +#include #include #include #include @@ -248,6 +249,7 @@ namespace LinearAlgebra */ template class Vector : public ::dealii::LinearAlgebra::VectorSpaceVector, + public ::dealii::ReadVector, public Subscriptor { public: @@ -1144,6 +1146,14 @@ namespace LinearAlgebra extract_subvector_to(const std::vector &indices, std::vector & values) const; + /** + * Extract a range of elements all at once. + */ + virtual void + extract_subvector_to( + const ArrayView &indices, + ArrayView &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 + void + Vector::extract_subvector_to( + const ArrayView &indices, + ArrayView & 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 bool Vector::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, 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 &indices, std::vector & values) const; + /** + * Extract a range of elements all at once. + */ + virtual void + extract_subvector_to( + const ArrayView &indices, + ArrayView & 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 &indices, + ArrayView & elements) const + { + AssertDimension(indices.size(), elements.size()); + extract_subvector_to(indices.begin(), indices.end(), elements.begin()); + } + + template 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 + +#include +#include + +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 +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 &indices, + ArrayView &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 #include +#include #include #include @@ -133,7 +134,7 @@ namespace LinearAlgebra * get the first index of the largest range. */ template - class ReadWriteVector : public Subscriptor + class ReadWriteVector : public Subscriptor, public ReadVector { 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 &indices, std::vector & values) const; + /** + * Extract a range of elements all at once. + */ + virtual void + extract_subvector_to( + const ArrayView &indices, + ArrayView & 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 + void + ReadWriteVector::extract_subvector_to( + const ArrayView &indices, + ArrayView & entries) const + { + AssertDimension(indices.size(), entries.size()); + for (unsigned int i = 0; i < indices.size(); ++i) + entries[i] = (*this)[indices[i]]; + } + + + template template 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 # include +# include # include # include # include @@ -61,9 +62,14 @@ namespace LinearAlgebra * @ingroup TrilinosWrappers * @ingroup Vectors */ - class Vector : public VectorSpaceVector, public Subscriptor + class Vector : public VectorSpaceVector, + public ReadVector, + 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 &V, const bool omit_zeroing_entries = false) override; + /** + * Extract a range of elements all at once. + */ + virtual void + extract_subvector_to( + const ArrayView &indices, + ArrayView &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 # include +# include # include # include # include @@ -113,7 +114,9 @@ namespace LinearAlgebra * @ingroup Vectors */ template - class Vector : public VectorSpaceVector, public Subscriptor + class Vector : public VectorSpaceVector, + public ReadVector, + public Subscriptor { public: using value_type = Number; @@ -159,6 +162,14 @@ namespace LinearAlgebra reinit(const VectorSpaceVector &V, const bool omit_zeroing_entries = false) override; + /** + * Extract a range of elements all at once. + */ + virtual void + extract_subvector_to( + const ArrayView &indices, + ArrayView &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 + void + Vector::extract_subvector_to( + const ArrayView &indices, + ArrayView & 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( + Tpetra::Access::ReadOnly); +# else + vector.template sync(); + auto vector_2d = vector.template getLocalView(); +# 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(indices[i])); + elements[i] = vector_1d(trilinos_i); + } + } + + + template Vector & Vector::operator=(const Vector &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 # include +# include # include # include # include @@ -394,7 +395,7 @@ namespace TrilinosWrappers * @ingroup TrilinosWrappers * @ingroup Vectors */ - class Vector : public Subscriptor + class Vector : public Subscriptor, public ReadVector { 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 &indices, std::vector & values) const; + /** + * Extract a range of elements all at once. + */ + virtual void + extract_subvector_to(const ArrayView &indices, + ArrayView &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 &indices, + ArrayView &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 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 #include +#include #include #include @@ -105,7 +106,7 @@ namespace parallel * in the manual). */ template -class Vector : public Subscriptor +class Vector : public Subscriptor, public ReadVector { public: /** @@ -662,6 +663,13 @@ public: extract_subvector_to(const std::vector &indices, std::vector & values) const; + /** + * Extract a range of elements all at once. + */ + virtual void + extract_subvector_to(const ArrayView &indices, + ArrayView &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::add_and_dot(const Number a, +template +void +Vector::extract_subvector_to( + const ArrayView &indices, + ArrayView & 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 Vector & Vector::operator+=(const Vector &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 &indices, + ArrayView & 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(indices[i])); + elements[i] = vector[0][trilinos_i]; + } + } + + + Vector & Vector::operator=(const Vector &V) {