*/
namespace EpetraWrappers
{
+ /**
+ * @cond internal
+ */
+
+ /**
+ * A namespace for internal implementation details of the TrilinosWrapper
+ * members.
+ *
+ * @ingroup TrilinosWrappers
+ */
+ namespace internal
+ {
+ /**
+ * Declare type for container size.
+ */
+ using size_type = dealii::types::global_dof_index;
+
+ /**
+ * Declare type for container value type.
+ */
+ using value_type = double;
+
+ /**
+ * This class implements a wrapper for accessing the Trilinos Epetra
+ * vector in the same way as we access deal.II objects: it is initialized
+ * with a vector and an element within it, and has a conversion operator
+ * to extract the scalar value of this element. It also has a variety of
+ * assignment operator for writing to this one element.
+ *
+ * @ingroup TrilinosWrappers
+ */
+ class VectorReference
+ {
+ private:
+ /**
+ * Constructor. It is made private so as to only allow the actual vector
+ * class to create it.
+ */
+ VectorReference(Vector &vector, const size_type index);
+
+ public:
+ /**
+ * Copy constructor.
+ */
+ VectorReference(const VectorReference &) = default;
+
+ /**
+ * This looks like a copy operator, but does something different than
+ * usual. In particular, it does not copy the member variables of this
+ * reference. Rather, it handles the situation where we have two vectors
+ * @p v and @p w, and assign elements like in <tt>v(i)=w(i)</tt>. Here,
+ * both left and right hand side of the assignment have data type
+ * VectorReference, but what we really mean is to assign the vector
+ * elements represented by the two references. This operator implements
+ * this operation. Note also that this allows us to make the assignment
+ * operator const.
+ */
+ const VectorReference &
+ operator=(const VectorReference &r) const;
+
+ /**
+ * Same as above but for non-const reference objects.
+ */
+ VectorReference &
+ operator=(const VectorReference &r);
+
+ /**
+ * Set the referenced element of the vector to <tt>s</tt>.
+ */
+ const VectorReference &
+ operator=(const value_type &s) const;
+
+ /**
+ * Add <tt>s</tt> to the referenced element of the vector->
+ */
+ const VectorReference &
+ operator+=(const value_type &s) const;
+
+ /**
+ * Subtract <tt>s</tt> from the referenced element of the vector->
+ */
+ const VectorReference &
+ operator-=(const value_type &s) const;
+
+ /**
+ * Multiply the referenced element of the vector by <tt>s</tt>.
+ */
+ const VectorReference &
+ operator*=(const value_type &s) const;
+
+ /**
+ * Divide the referenced element of the vector by <tt>s</tt>.
+ */
+ const VectorReference &
+ operator/=(const value_type &s) const;
+
+ /**
+ * Convert the reference to an actual value, i.e. return the value of
+ * the referenced element of the vector.
+ */
+ operator value_type() const;
+
+ /**
+ * Exception
+ */
+ DeclException1(ExcTrilinosError,
+ int,
+ << "An error with error number " << arg1
+ << " occurred while calling a Trilinos function");
+
+ /*
+ * Access to a an element that is not (locally-)owned.
+ *
+ * @ingroup Exceptions
+ */
+ DeclException4(
+ ExcAccessToNonLocalElement,
+ size_type,
+ size_type,
+ size_type,
+ size_type,
+ << "You are trying to access element " << arg1
+ << " of a distributed vector, but this element is not stored "
+ << "on the current processor. Note: There are " << arg2
+ << " elements stored "
+ << "on the current processor from within the range [" << arg3 << ','
+ << arg4 << "] but Trilinos vectors need not store contiguous "
+ << "ranges on each processor, and not every element in "
+ << "this range may in fact be stored locally."
+ << "\n\n"
+ << "A common source for this kind of problem is that you "
+ << "are passing a 'fully distributed' vector into a function "
+ << "that needs read access to vector elements that correspond "
+ << "to degrees of freedom on ghost cells (or at least to "
+ << "'locally active' degrees of freedom that are not also "
+ << "'locally owned'). You need to pass a vector that has these "
+ << "elements as ghost entries.");
+
+ private:
+ /**
+ * Point to the vector we are referencing.
+ */
+ Vector &vector;
+
+ /**
+ * Index of the referenced element of the vector.
+ */
+ const size_type index;
+
+ // Make the vector class a friend, so that it can create objects of the
+ // present type.
+ friend class ::dealii::LinearAlgebra::EpetraWrappers::Vector;
+ }; // class VectorReference
+
+ } // namespace internal
+ /**
+ * @endcond
+ */
+
+
/**
* This class implements a wrapper to the Trilinos distributed vector
* class Epetra_FEVector. This class requires Trilinos to be compiled
* @ingroup TrilinosWrappers
* @ingroup Vectors
*/
- class Vector : public ReadVector<double>, public Subscriptor
+ class Vector : public ReadVector<internal::value_type>, public Subscriptor
{
public:
- using value_type = double;
- using size_type = types::global_dof_index;
- using real_type = double;
+ using value_type = internal::value_type;
+ using size_type = types::global_dof_index;
+ using reference = internal::VectorReference;
+ using const_reference = const internal::VectorReference;
+
+ /**
+ * @name 1: Basic Object-handling
+ */
+ /** @{ */
/**
* Constructor. Create a vector of dimension zero.
import_elements(V, operation, communication_pattern);
}
+ /** @} */
+
+
+ /**
+ * @name 2: Data-Access
+ */
+ /** @{ */
+
+ /**
+ * Provide access to a given element, both read and write.
+ *
+ * When using a vector distributed with MPI, this operation only makes
+ * sense for elements that are actually present on the calling processor.
+ * Otherwise, an exception is thrown.
+ */
+ reference
+ operator()(const size_type index);
+
+ /**
+ * Provide read-only access to an element.
+ *
+ * When using a vector distributed with MPI, this operation only makes
+ * sense for elements that are actually present on the calling processor.
+ * Otherwise, an exception is thrown.
+ */
+ value_type
+ operator()(const size_type index) const;
+
+ /**
+ * Provide access to a given element, both read and write.
+ *
+ * Exactly the same as operator().
+ */
+ reference
+ operator[](const size_type index);
+
+ /**
+ * Provide read-only access to an element.
+ *
+ * Exactly the same as operator().
+ */
+ value_type
+ operator[](const size_type index) const;
+
+ /** @} */
+
+ /**
+ * @name 3: Modification of vectors
+ */
+ /** @{ */
+
/**
* Multiply the entire vector by a fixed factor.
*/
bool
all_zero() const;
+ /** @} */
+
+
+ /**
+ * @name 4: Scalar products, norms and related operations
+ */
+ /** @{ */
+
/**
* Return the mean value of the element of this vector.
*/
*/
double
add_and_dot(const double a, const Vector &V, const Vector &W);
+
+ /** @} */
+
+ /**
+ * @name 5: Scalar products, norms and related operations
+ */
+ /** @{ */
+
/**
* This function always returns false and is present only for backward
* compatibility.
::dealii::IndexSet
locally_owned_elements() const;
+ /** @} */
+
+
+ /**
+ * @name 6: Mixed stuff
+ */
+ /** @{ */
+
/**
* Compress the underlying representation of the Trilinos object, i.e.
* flush the buffers of the vector object if it has any. This function is
std::size_t
memory_consumption() const;
+ /** @} */
+
/**
* The vectors have different partitioning, i.e. their IndexSet objects
* don't represent the same indices.
* source_stored_elements IndexSet and the current vector.
*/
std::shared_ptr<const CommunicationPattern> epetra_comm_pattern;
+
+ // Make the reference class a friend.
+ friend class internal::VectorReference;
};
+# ifndef DOXYGEN
+
+ // VectorReference
+ namespace internal
+ {
+ inline VectorReference::VectorReference(Vector &vector,
+ const size_type index)
+ : vector(vector)
+ , index(index)
+ {}
+
+
+
+ inline const VectorReference &
+ VectorReference::operator=(const VectorReference &r) const
+ {
+ // as explained in the class
+ // documentation, this is not the copy
+ // operator. so simply pass on to the
+ // "correct" assignment operator
+ *this = static_cast<value_type>(r);
+
+ return *this;
+ }
+
+
+
+ inline VectorReference &
+ VectorReference::operator=(const VectorReference &r)
+ {
+ // as above
+ *this = static_cast<value_type>(r);
+
+ return *this;
+ }
+
+
+
+ inline const VectorReference &
+ VectorReference::operator=(const value_type &value) const
+ {
+ (*vector.vector)[0][index] = value;
+
+ return *this;
+ }
+
+
+
+ inline const VectorReference &
+ VectorReference::operator+=(const value_type &value) const
+ {
+ value_type new_value = static_cast<value_type>(*this) + value;
+ (*vector.vector)[0][index] = new_value;
+
+ return *this;
+ }
+
+
+
+ inline const VectorReference &
+ VectorReference::operator-=(const value_type &value) const
+ {
+ value_type new_value = static_cast<value_type>(*this) - value;
+ (*vector.vector)[0][index] = new_value;
+
+ return *this;
+ }
+
+
+
+ inline const VectorReference &
+ VectorReference::operator*=(const value_type &value) const
+ {
+ value_type new_value = static_cast<value_type>(*this) * value;
+ (*vector.vector)[0][index] = new_value;
+
+ return *this;
+ }
+
+
+
+ inline const VectorReference &
+ VectorReference::operator/=(const value_type &value) const
+ {
+ value_type new_value = static_cast<value_type>(*this) / value;
+ (*vector.vector)[0][index] = new_value;
+
+ return *this;
+ }
+ } // namespace internal
+
+# endif /* DOXYGEN */
+
+
+ inline internal::VectorReference
+ Vector::operator()(const size_type index)
+ {
+ return internal::VectorReference(*this, index);
+ }
+
+ inline internal::VectorReference
+ Vector::operator[](const size_type index)
+ {
+ return operator()(index);
+ }
+
+ inline Vector::value_type
+ Vector::operator[](const size_type index) const
+ {
+ return operator()(index);
+ }
+
inline bool
Vector::has_ghost_elements() const
--- /dev/null
+// ------------------------------------------------------------------------
+//
+// SPDX-License-Identifier: LGPL-2.1-or-later
+// Copyright (C) 2017 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// Part of the source code is dual licensed under Apache-2.0 WITH
+// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
+// governing the source code and code contributions can be found in
+// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
+//
+// ------------------------------------------------------------------------
+
+
+
+// LinearAlgebra::EpetraWrappers::internal::VectorReference had its non-const
+// copy operator return a const reference. That was non-intuitive but, because
+// of the particular semantics of copying vector reference objects, made no
+// difference. Either way, we now check these semantics.
+
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/lac/trilinos_epetra_vector.h>
+
+#include <iostream>
+
+#include "../tests.h"
+
+
+void
+test()
+{
+ LinearAlgebra::EpetraWrappers::Vector v;
+ v.reinit(complete_index_set(3), MPI_COMM_WORLD);
+ v(0) = 0;
+ v(1) = 1;
+ v(2) = 2;
+
+ LinearAlgebra::EpetraWrappers::internal::VectorReference a(v(0));
+ LinearAlgebra::EpetraWrappers::internal::VectorReference b(v(1));
+ LinearAlgebra::EpetraWrappers::internal::VectorReference c(v(2));
+
+ // Copy the VectorReference objects. Note that operator= for these
+ // objects does not copy the *content* of the reference object, but
+ // rather assigns the value of the vector element referenced on the
+ // right to the vector element referenced on the left.
+ // So the applied operations result in the following actions:
+ // 1. We first assign c to point to where a points (c points to entry 0)
+ // 2. We next assign c to point to where b points (c points to entry 1)
+ // 3. Lastly, we assign b to point to where a points (b points to entry 0)
+ (c = a) = b;
+ b = a;
+
+ deallog << static_cast<TrilinosScalar>(a)
+ << std::endl; // should point to v(0)
+ deallog << static_cast<TrilinosScalar>(b)
+ << std::endl; // should point to v(0)
+ deallog << static_cast<TrilinosScalar>(c)
+ << std::endl; // should point to v(1)
+}
+
+
+
+int
+main(int argc, char **argv)
+{
+ Utilities::MPI::MPI_InitFinalize mpi_initialization(
+ argc, argv, testing_max_num_threads());
+ initlog();
+
+ test();
+}
projects / dealii.git / commitdiff