From: Bruno Turcksin <bruno.turcksin@gmail.com>
Date: Fri, 10 Feb 2017 14:27:10 +0000 (-0500)
Subject: Instantiate LinearAlgebra::EpetraWrappers::Vector.
X-Git-Tag: v9.0.0-rc1~1682^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F3957%2Fhead;p=dealii.git

Instantiate LinearAlgebra::EpetraWrappers::Vector.
---

diff --git a/cmake/config/template-arguments.in b/cmake/config/template-arguments.in
index 0ef726287e..312a009259 100644
--- a/cmake/config/template-arguments.in
+++ b/cmake/config/template-arguments.in
@@ -46,6 +46,7 @@ SERIAL_VECTORS := { Vector<double>;
 
                     @DEAL_II_EXPAND_TRILINOS_VECTOR@;
                     @DEAL_II_EXPAND_TRILINOS_MPI_VECTOR@;
+                    @DEAL_II_EXPAND_EPETRA_VECTOR@;
                     @DEAL_II_EXPAND_PETSC_VECTOR@;
                     @DEAL_II_EXPAND_PETSC_MPI_VECTOR@;
 
@@ -73,6 +74,7 @@ REAL_SERIAL_VECTORS := { Vector<double>;
 
                     @DEAL_II_EXPAND_TRILINOS_VECTOR@;
                     @DEAL_II_EXPAND_TRILINOS_MPI_VECTOR@;
+                    @DEAL_II_EXPAND_EPETRA_VECTOR@;
                     @DEAL_II_EXPAND_PETSC_VECTOR_REAL@;
                     @DEAL_II_EXPAND_PETSC_MPI_VECTOR_REAL@;
 
@@ -94,6 +96,7 @@ REAL_NONBLOCK_VECTORS := { Vector<double>;
 
                     @DEAL_II_EXPAND_TRILINOS_VECTOR@;
                     @DEAL_II_EXPAND_TRILINOS_MPI_VECTOR@;
+                    @DEAL_II_EXPAND_EPETRA_VECTOR@;
                     @DEAL_II_EXPAND_PETSC_VECTOR_REAL@;
                     @DEAL_II_EXPAND_PETSC_MPI_VECTOR_REAL@;
                   }
@@ -108,6 +111,7 @@ EXTERNAL_SEQUENTIAL_VECTORS := { @DEAL_II_EXPAND_TRILINOS_VECTOR@;
 // wrappers for MPI vectors (PETSc/Trilinos) 
 EXTERNAL_PARALLEL_VECTORS := { @DEAL_II_EXPAND_TRILINOS_MPI_VECTOR@;
                                @DEAL_II_EXPAND_TRILINOS_MPI_BLOCKVECTOR@;
+                               @DEAL_II_EXPAND_EPETRA_VECTOR@;
                                @DEAL_II_EXPAND_PETSC_MPI_VECTOR@;
                                @DEAL_II_EXPAND_PETSC_MPI_BLOCKVECTOR@
                              }
@@ -126,6 +130,7 @@ VECTORS_WITH_MATRIX := { Vector<double>;
 
                     @DEAL_II_EXPAND_TRILINOS_VECTOR@;
                     @DEAL_II_EXPAND_TRILINOS_MPI_VECTOR@;
+                    @DEAL_II_EXPAND_EPETRA_VECTOR@;
                   }
 
 // DoFHandler and hp::DoFHandler templates
diff --git a/cmake/configure/configure_2_trilinos.cmake b/cmake/configure/configure_2_trilinos.cmake
index 67dcedee07..e00ed23534 100644
--- a/cmake/configure/configure_2_trilinos.cmake
+++ b/cmake/configure/configure_2_trilinos.cmake
@@ -194,6 +194,7 @@ MACRO(FEATURE_TRILINOS_CONFIGURE_EXTERNAL)
   SET(DEAL_II_EXPAND_TRILINOS_BLOCK_SPARSITY_PATTERN "TrilinosWrappers::BlockSparsityPattern")
   SET(DEAL_II_EXPAND_TRILINOS_MPI_BLOCKVECTOR "TrilinosWrappers::MPI::BlockVector")
   SET(DEAL_II_EXPAND_TRILINOS_MPI_VECTOR "TrilinosWrappers::MPI::Vector")
+  SET(DEAL_II_EXPAND_EPETRA_VECTOR "LinearAlgebra::EpetraWrappers::Vector")
 ENDMACRO()
 
 
diff --git a/include/deal.II/algorithms/operator.templates.h b/include/deal.II/algorithms/operator.templates.h
index 7b82ebe63f..26d1893639 100644
--- a/include/deal.II/algorithms/operator.templates.h
+++ b/include/deal.II/algorithms/operator.templates.h
@@ -19,6 +19,7 @@
 
 #include <deal.II/algorithms/operator.h>
 #include <deal.II/base/logstream.h>
+#include <deal.II/lac/vector_element_access.h>
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -54,7 +55,8 @@ namespace Algorithms
             if (v == nullptr) continue;
             deallog << vectors.name(i);
             for (unsigned int j=0; j<v->size(); ++j)
-              deallog << ' ' << (*v)(j);
+              deallog << ' ' << ::dealii::internal::ElementAccess<VectorType>::get(
+                        *v, j);
             deallog << std::endl;
           }
         deallog << std::endl;
@@ -67,7 +69,8 @@ namespace Algorithms
             const VectorType *v = vectors.try_read_ptr<VectorType>(i);
             if (v == nullptr) continue;
             for (unsigned int j=0; j<v->size(); ++j)
-              (*os) << ' ' << (*v)(j);
+              (*os) << ' ' << ::dealii::internal::ElementAccess<VectorType>::get(
+                      *v, j);
           }
         (*os) << std::endl;
       }
diff --git a/include/deal.II/dofs/dof_accessor.templates.h b/include/deal.II/dofs/dof_accessor.templates.h
index 7c45c37a3c..fa19f0c392 100644
--- a/include/deal.II/dofs/dof_accessor.templates.h
+++ b/include/deal.II/dofs/dof_accessor.templates.h
@@ -20,6 +20,7 @@
 #include <deal.II/base/config.h>
 #include <deal.II/base/types.h>
 #include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/read_write_vector.h>
 #include <deal.II/dofs/dof_accessor.h>
 #include <deal.II/dofs/dof_levels.h>
 #include <deal.II/dofs/dof_faces.h>
@@ -1479,6 +1480,66 @@ namespace internal
         Assert (next == dof_indices.end (), ExcInternalError ());
       }
 
+
+
+      template <typename InputVector, typename ForwardIterator>
+      static
+      void extract_subvector_to(const InputVector &values,
+                                const types::global_dof_index *cache,
+                                const types::global_dof_index *cache_end,
+                                ForwardIterator local_values_begin)
+      {
+        values.extract_subvector_to (cache, cache_end, local_values_begin);
+      }
+
+
+
+#ifdef DEAL_II_WITH_TRILINOS
+      static
+      std::vector<unsigned int> sort_indices(const types::global_dof_index *v_begin,
+                                             const types::global_dof_index *v_end)
+      {
+        // initialize original index locations
+        std::vector<unsigned int> idx(v_end-v_begin);
+        std::iota(idx.begin(), idx.end(), 0);
+
+        // sort indices based on comparing values in v
+        std::sort(idx.begin(), idx.end(),
+                  [&v_begin](unsigned int i1, unsigned int i2)
+        {
+          return *(v_begin+i1) < *(v_begin+i2);
+        });
+
+        return idx;
+      }
+
+
+
+      template <typename ForwardIterator>
+      static
+      void extract_subvector_to(const LinearAlgebra::EpetraWrappers::Vector &values,
+                                const types::global_dof_index *cache_begin,
+                                const types::global_dof_index *cache_end,
+                                ForwardIterator local_values_begin)
+      {
+        std::vector<unsigned int> sorted_indices_pos = sort_indices(cache_begin,
+                                                                    cache_end);
+        const unsigned int cache_size = cache_end-cache_begin;
+        std::vector<types::global_dof_index> cache_indices(cache_size);
+        for (unsigned int i=0; i < cache_size; ++i)
+          cache_indices[i] = *(cache_begin+sorted_indices_pos[i]);
+
+        IndexSet index_set(cache_indices.back() + 1);
+        index_set.add_indices(cache_indices.begin(), cache_indices.end());
+        index_set.compress();
+        LinearAlgebra::ReadWriteVector<double> read_write_vector(index_set);
+        read_write_vector.import(values, VectorOperation::insert);
+
+        // Copy the elements from read_write_vector and reorder them.
+        for (unsigned int i=0; i<cache_size; ++i, ++local_values_begin)
+          *local_values_begin = read_write_vector[sorted_indices_pos[i]];
+      }
+#endif
     };
   }
 }
@@ -3530,10 +3591,8 @@ DoFCellAccessor<DoFHandlerType,level_dof_access>::get_dof_values
   const types::global_dof_index *cache
     = this->dof_handler->levels[this->present_level]
       ->get_cell_cache_start (this->present_index, this->get_fe().dofs_per_cell);
-
-  values.extract_subvector_to (cache,
-                               cache + this->get_fe().dofs_per_cell,
-                               local_values_begin);
+  dealii::internal::DoFAccessor::Implementation::extract_subvector_to(
+    values, cache, cache + this->get_fe().dofs_per_cell, local_values_begin);
 }
 
 
@@ -3595,7 +3654,8 @@ DoFCellAccessor<DoFHandlerType,level_dof_access>::set_dof_values
       ->get_cell_cache_start (this->present_index, this->get_fe().dofs_per_cell);
 
   for (unsigned int i=0; i<this->get_fe().dofs_per_cell; ++i, ++cache)
-    values(*cache) = local_values(i);
+    internal::ElementAccess<OutputVector>::set(local_values(i),
+                                               *cache, values);
 }
 
 
diff --git a/include/deal.II/fe/fe_tools_extrapolate.templates.h b/include/deal.II/fe/fe_tools_extrapolate.templates.h
index a3671bfe5b..eee5051ef4 100644
--- a/include/deal.II/fe/fe_tools_extrapolate.templates.h
+++ b/include/deal.II/fe/fe_tools_extrapolate.templates.h
@@ -684,7 +684,8 @@ namespace FETools
                   const bool on_refined_neighbor
                     = (dofs_on_refined_neighbors.find(indices[j])!=dofs_on_refined_neighbors.end());
                   if (!(on_refined_neighbor && dofs_on_refined_neighbors[indices[j]]>dealii_cell->level()))
-                    u(indices[j]) = local_values(j);
+                    ::dealii::internal::ElementAccess<OutVector>::set(
+                      local_values(j), indices[j], u);
                 }
             }
         }
@@ -1450,6 +1451,20 @@ namespace FETools
         const IndexSet locally_owned_dofs = dh.locally_owned_dofs();
         vector.reinit(locally_owned_dofs, parallel_tria->get_communicator());
       }
+
+
+
+      template <int dim, int spacedim>
+      void reinit_distributed(const DoFHandler<dim, spacedim> &dh,
+                              LinearAlgebra::EpetraWrappers::Vector &vector)
+      {
+        const parallel::distributed::Triangulation<dim,spacedim> *parallel_tria =
+          dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>(&dh.get_tria());
+        Assert (parallel_tria !=0, ExcNotImplemented());
+
+        const IndexSet locally_owned_dofs = dh.locally_owned_dofs();
+        vector.reinit(locally_owned_dofs, parallel_tria->get_communicator());
+      }
 #endif //DEAL_II_WITH_TRILINOS
 
       template <int dim, int spacedim, typename Number>
diff --git a/include/deal.II/fe/fe_tools_interpolate.templates.h b/include/deal.II/fe/fe_tools_interpolate.templates.h
index fc987af8d0..7fb69f1961 100644
--- a/include/deal.II/fe/fe_tools_interpolate.templates.h
+++ b/include/deal.II/fe/fe_tools_interpolate.templates.h
@@ -201,8 +201,10 @@ namespace FETools
               const types::global_dof_index gdi = dofs[i];
               if (u2_elements.is_element(gdi))
                 {
-                  u2(dofs[i])+=u2_local(i);
-                  touch_count(dofs[i]) += 1;
+                  ::dealii::internal::ElementAccess<OutVector>::add(u2_local(i),
+                                                                    dofs[i], u2);
+                  ::dealii::internal::ElementAccess<OutVector>::add(1,
+                                                                    dofs[i], touch_count);
                 }
             }
         }
@@ -227,9 +229,16 @@ namespace FETools
     for (types::global_dof_index i=0; i<dof2.n_dofs(); ++i)
       if (locally_owned_dofs.is_element(i))
         {
-          Assert(static_cast<typename OutVector::value_type>(touch_count(i)) != typename OutVector::value_type(0),
+          Assert(static_cast<typename OutVector::value_type>(
+                   ::dealii::internal::ElementAccess<OutVector>::get(
+                     touch_count, i)) != typename OutVector::value_type(0),
                  ExcInternalError());
-          u2(i) /= touch_count(i);
+
+
+          const double val = ::dealii::internal::ElementAccess<OutVector>::get(
+                               u2, i);
+          ::dealii::internal::ElementAccess<OutVector>::set(
+            val/::dealii::internal::ElementAccess<OutVector>::get(touch_count,i), i, u2);
         }
 
     // finish the work on parallel vectors
@@ -639,7 +648,8 @@ namespace FETools
         cell2->get_dof_indices(dofs);
         for (unsigned int i=0; i<n2; ++i)
           {
-            u2(dofs[i])+=u2_local(i);
+            ::dealii::internal::ElementAccess<OutVector>::add(u2_local(i),
+                                                              dofs[i], u2);
           }
 
         ++cell1;
diff --git a/include/deal.II/lac/constraint_matrix.h b/include/deal.II/lac/constraint_matrix.h
index 38eb9432fd..148640b06a 100644
--- a/include/deal.II/lac/constraint_matrix.h
+++ b/include/deal.II/lac/constraint_matrix.h
@@ -24,6 +24,7 @@
 #include <deal.II/base/template_constraints.h>
 
 #include <deal.II/lac/vector.h>
+#include <deal.II/lac/vector_element_access.h>
 
 #include <vector>
 #include <set>
@@ -1658,14 +1659,16 @@ void ConstraintMatrix::distribute_local_to_global (
         ++local_vector_begin, ++local_indices_begin)
     {
       if (is_constrained(*local_indices_begin) == false)
-        global_vector(*local_indices_begin) += *local_vector_begin;
+        internal::ElementAccess<VectorType>::add(*local_vector_begin,
+                                                 *local_indices_begin, global_vector);
       else
         {
           const ConstraintLine &position =
             lines[lines_cache[calculate_line_index(*local_indices_begin)]];
           for (size_type j=0; j<position.entries.size(); ++j)
-            global_vector(position.entries[j].first)
-            += *local_vector_begin * position.entries[j].second;
+            internal::ElementAccess<VectorType>::add(*local_vector_begin *
+                                                     position.entries[j].second, position.entries[j].first,
+                                                     global_vector);
         }
     }
 }
diff --git a/include/deal.II/lac/constraint_matrix.templates.h b/include/deal.II/lac/constraint_matrix.templates.h
index fdc1b1029d..5f5f979fb7 100644
--- a/include/deal.II/lac/constraint_matrix.templates.h
+++ b/include/deal.II/lac/constraint_matrix.templates.h
@@ -531,7 +531,7 @@ namespace internal
               continue;
             size_type idx = *it - shift;
             if (idx<vec.size() && locally_owned.is_element(idx))
-              vec(idx) = 0.;
+              internal::ElementAccess<VEC>::set(0., idx, vec);
           }
       }
 
@@ -895,10 +895,11 @@ ConstraintMatrix::distribute (VectorType &vec) const
             new_value = it->inhomogeneity;
             for (unsigned int i=0; i<it->entries.size(); ++i)
               new_value += (static_cast<typename VectorType::value_type>
-                            (ghosted_vector(it->entries[i].first)) *
+                            (internal::ElementAccess<VectorType>::get(
+                               ghosted_vector, it->entries[i].first)) *
                             it->entries[i].second);
             AssertIsFinite(new_value);
-            vec(it->line) = new_value;
+            internal::ElementAccess<VectorType>::set(new_value, it->line, vec);
           }
 
       // now compress to communicate the entries that we added to
@@ -923,10 +924,12 @@ ConstraintMatrix::distribute (VectorType &vec) const
           new_value = next_constraint->inhomogeneity;
           for (unsigned int i=0; i<next_constraint->entries.size(); ++i)
             new_value += (static_cast<typename VectorType::value_type>
-                          (vec(next_constraint->entries[i].first)) *
+                          (internal::ElementAccess<VectorType>::get(
+                             vec, next_constraint->entries[i].first))*
                           next_constraint->entries[i].second);
           AssertIsFinite(new_value);
-          vec(next_constraint->line) = new_value;
+          internal::ElementAccess<VectorType>::set(new_value, next_constraint->line,
+                                                   vec);
         }
     }
 }
diff --git a/include/deal.II/lac/trilinos_epetra_vector.h b/include/deal.II/lac/trilinos_epetra_vector.h
index 1745f94768..803cc2c230 100644
--- a/include/deal.II/lac/trilinos_epetra_vector.h
+++ b/include/deal.II/lac/trilinos_epetra_vector.h
@@ -227,7 +227,7 @@ namespace LinearAlgebra
                                  const VectorSpaceVector<double> &W);
       /**
        * This function always returns false and is present only for backward
-       * compatibily.
+       * compatibility.
        */
       bool has_ghost_elements() const;
 
@@ -342,7 +342,7 @@ namespace LinearAlgebra
  * Declare dealii::LinearAlgebra::EpetraWrappers::Vector as distributed vector.
  */
 template <>
-struct is_serial_vector<LinearAlgebra::EpetraWrappers::Vector> : std_cxx11::false_type
+struct is_serial_vector<LinearAlgebra::EpetraWrappers::Vector> : std::false_type
 {
 };
 
diff --git a/include/deal.II/lac/vector_element_access.h b/include/deal.II/lac/vector_element_access.h
new file mode 100644
index 0000000000..1ad1054ad2
--- /dev/null
+++ b/include/deal.II/lac/vector_element_access.h
@@ -0,0 +1,120 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii__vector_element_access_h
+#define dealii__vector_element_access_h
+
+
+#include <deal.II/lac/trilinos_epetra_vector.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace internal
+{
+  template <typename VectorType>
+  struct ElementAccess
+  {
+  public:
+    static void add(typename VectorType::value_type value,
+                    types::global_dof_index i, VectorType &V);
+
+    static void set(typename VectorType::value_type value,
+                    types::global_dof_index i, VectorType &V);
+
+    static typename VectorType::value_type get(const VectorType &V,
+                                               types::global_dof_index i);
+  };
+
+
+
+  template <typename VectorType>
+  inline
+  void ElementAccess<VectorType>::add(typename VectorType::value_type value,
+                                      types::global_dof_index i, VectorType &V)
+  {
+    V(i) += value;
+  }
+
+
+
+  template <typename VectorType>
+  inline
+  void ElementAccess<VectorType>::set(typename VectorType::value_type value,
+                                      types::global_dof_index i, VectorType &V)
+  {
+    V(i) = value;
+  }
+
+
+
+  template <typename VectorType>
+  inline
+  typename VectorType::value_type
+  ElementAccess<VectorType>::get(const VectorType &V, types::global_dof_index i)
+  {
+    return V(i);
+  }
+
+
+
+#ifdef DEAL_II_WITH_TRILINOS
+  template <>
+  inline
+  void ElementAccess<LinearAlgebra::EpetraWrappers::Vector>::
+  add(double value, types::global_dof_index i, LinearAlgebra::EpetraWrappers::Vector &V)
+  {
+    // Extract local indices in the vector.
+    Epetra_FEVector vector = V.trilinos_vector();
+    TrilinosWrappers::types::int_type trilinos_i =
+      vector.Map().LID(static_cast<TrilinosWrappers::types::int_type>(i));
+
+    vector[0][trilinos_i] += value;
+  }
+
+
+
+  template <>
+  inline
+  void ElementAccess<LinearAlgebra::EpetraWrappers::Vector>::
+  set(double value, types::global_dof_index i, LinearAlgebra::EpetraWrappers::Vector &V)
+  {
+    // Extract local indices in the vector.
+    Epetra_FEVector vector = V.trilinos_vector();
+    TrilinosWrappers::types::int_type trilinos_i =
+      vector.Map().LID(static_cast<TrilinosWrappers::types::int_type>(i));
+
+    vector[0][trilinos_i] = value;
+  }
+
+
+  template <>
+  inline
+  double
+  ElementAccess<LinearAlgebra::EpetraWrappers::Vector>::
+  get(const LinearAlgebra::EpetraWrappers::Vector &V, types::global_dof_index i)
+  {
+    // Extract local indices in the vector.
+    Epetra_FEVector vector = V.trilinos_vector();
+    TrilinosWrappers::types::int_type trilinos_i =
+      vector.Map().LID(static_cast<TrilinosWrappers::types::int_type>(i));
+
+    return vector[0][trilinos_i];
+  }
+#endif
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
diff --git a/include/deal.II/multigrid/mg_transfer.h b/include/deal.II/multigrid/mg_transfer.h
index b102d3fa46..b585b17560 100644
--- a/include/deal.II/multigrid/mg_transfer.h
+++ b/include/deal.II/multigrid/mg_transfer.h
@@ -72,6 +72,7 @@ namespace internal
     }
 
   };
+
   template <>
   struct MatrixSelector<dealii::TrilinosWrappers::MPI::Vector>
   {
@@ -88,6 +89,22 @@ namespace internal
 
   };
 
+  template <>
+  struct MatrixSelector<dealii::LinearAlgebra::EpetraWrappers::Vector>
+  {
+    typedef ::dealii::TrilinosWrappers::SparsityPattern Sparsity;
+    typedef ::dealii::TrilinosWrappers::SparseMatrix Matrix;
+
+    template <typename SparsityPatternType, typename DoFHandlerType>
+    static void reinit(Matrix &matrix, Sparsity &, int level, const SparsityPatternType &sp, DoFHandlerType &dh)
+    {
+      matrix.reinit(dh.locally_owned_mg_dofs(level+1),
+                    dh.locally_owned_mg_dofs(level),
+                    sp, MPI_COMM_WORLD, true);
+    }
+
+  };
+
   template <>
   struct MatrixSelector<dealii::TrilinosWrappers::Vector>
   {
diff --git a/include/deal.II/multigrid/mg_transfer.templates.h b/include/deal.II/multigrid/mg_transfer.templates.h
index 68d7b87ced..1fdc9232b9 100644
--- a/include/deal.II/multigrid/mg_transfer.templates.h
+++ b/include/deal.II/multigrid/mg_transfer.templates.h
@@ -18,6 +18,7 @@
 #define dealii__mg_transfer_templates_h
 
 #include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/lac/trilinos_epetra_vector.h>
 #include <deal.II/lac/sparse_matrix.h>
 #include <deal.II/lac/vector_view.h>
 #include <deal.II/grid/tria_iterator.h>
diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h
index dea18e79b1..86be34fa0d 100644
--- a/include/deal.II/numerics/vector_tools.templates.h
+++ b/include/deal.II/numerics/vector_tools.templates.h
@@ -30,6 +30,7 @@
 #include <deal.II/lac/la_parallel_block_vector.h>
 #include <deal.II/lac/petsc_vector.h>
 #include <deal.II/lac/petsc_block_vector.h>
+#include <deal.II/lac/trilinos_epetra_vector.h>
 #include <deal.II/lac/trilinos_vector.h>
 #include <deal.II/lac/trilinos_block_vector.h>
 #include <deal.II/lac/sparse_matrix.h>
@@ -264,8 +265,9 @@ namespace VectorTools
                       if (component_mask[component] == true)
                         {
                           const unsigned int rep_dof=dof_to_rep_index_table[fe_index][i];
-                          vec(dofs_on_cell[i])
-                            = function_values_system[fe_index][rep_dof](component);
+                          ::dealii::internal::ElementAccess<VectorType>::set(
+                            function_values_system[fe_index][rep_dof](component),
+                            dofs_on_cell[i], vec);
                         }
                     }
                 }
@@ -282,8 +284,9 @@ namespace VectorTools
                   // values to the global
                   // vector
                   for (unsigned int i=0; i<fe[fe_index].dofs_per_cell; ++i)
-                    vec(dofs_on_cell[i])
-                      = function_values_scalar[fe_index][dof_to_rep_index_table[fe_index][i]];
+                    ::dealii::internal::ElementAccess<VectorType>::set(
+                      function_values_scalar[fe_index][dof_to_rep_index_table[fe_index][i]],
+                      dofs_on_cell[i], vec);
                 }
             }
         }
@@ -337,7 +340,8 @@ namespace VectorTools
         // distribute cell vector
         for (unsigned int j=0; j<dof_2.get_fe().dofs_per_cell; ++j)
           {
-            data_2(local_dof_indices[j]) += cell_data_2(j);
+            ::dealii::internal::ElementAccess<OutVector>::add(cell_data_2(j),
+                                                              local_dof_indices[j], data_2);
 
             // count, how often we have
             // added to this dof
@@ -355,7 +359,11 @@ namespace VectorTools
         Assert (touch_count[i] != 0,
                 ExcInternalError());
 
-        data_2(i) /= touch_count[i];
+        const double val = ::dealii::internal::ElementAccess<OutVector>::get(
+                             data_2, i);
+
+        ::dealii::internal::ElementAccess<OutVector>::set(
+          val/touch_count[i], i, data_2);
       }
   }
 
@@ -495,7 +503,9 @@ namespace VectorTools
                     if ( component_mask[component] )
                       {
                         const unsigned int rep_dof = dof_to_rep_index_table[fe_index][i];
-                        dst(dofs_on_cell[i])       = function_values_system[fe_index][rep_dof](component);
+                        ::dealii::internal::ElementAccess<VectorType>::set(
+                          function_values_system[fe_index][rep_dof](component),
+                          dofs_on_cell[i], dst);
                       }
                   }
               }
@@ -509,8 +519,9 @@ namespace VectorTools
                  0);
 
                 for (unsigned int i = 0; i < fe[fe_index].dofs_per_cell; ++i)
-                  dst(dofs_on_cell[i]) = function_values_scalar[fe_index]
-                                         [dof_to_rep_index_table[fe_index][i]];
+                  ::dealii::internal::ElementAccess<VectorType>::set(
+                    function_values_scalar[fe_index][dof_to_rep_index_table[fe_index][i]],
+                    dofs_on_cell[i], dst);
               }
           }
 
@@ -890,7 +901,8 @@ namespace VectorTools
       // it may be of another type than Vector<double> and that wouldn't
       // necessarily go together with the matrix and other functions
       for (unsigned int i=0; i<vec.size(); ++i)
-        vec_result(i) = vec(i);
+        ::dealii::internal::ElementAccess<VectorType>::set(vec(i), i,
+                                                           vec_result);
     }
 
 
@@ -1090,7 +1102,8 @@ namespace VectorTools
       const IndexSet locally_owned_dofs = dof.locally_owned_dofs();
       IndexSet::ElementIterator it = locally_owned_dofs.begin();
       for (; it!=locally_owned_dofs.end(); ++it)
-        vec_result(*it) = work_result(*it);
+        ::dealii::internal::ElementAccess<VectorType>::set(work_result(*it),
+                                                           *it, vec_result);
       vec_result.compress(VectorOperation::insert);
     }
 
@@ -1235,7 +1248,7 @@ namespace VectorTools
       const IndexSet locally_owned_dofs = dof.locally_owned_dofs();
       IndexSet::ElementIterator it = locally_owned_dofs.begin();
       for (; it!=locally_owned_dofs.end(); ++it)
-        vec_result(*it) = vec(*it);
+        ::dealii::internal::ElementAccess<VectorType>::set(vec(*it), *it, vec_result);
       vec_result.compress(VectorOperation::insert);
     }
 
@@ -1310,7 +1323,7 @@ namespace VectorTools
       const IndexSet locally_owned_dofs = dof.locally_owned_dofs();
       IndexSet::ElementIterator it = locally_owned_dofs.begin();
       for (; it!=locally_owned_dofs.end(); ++it)
-        vec_result(*it) = vec(*it);
+        ::dealii::internal::ElementAccess<VectorType>::set(vec(*it), *it, vec_result);
       vec_result.compress(VectorOperation::insert);
     }
   }
@@ -7933,7 +7946,8 @@ namespace VectorTools
               {
                 unsigned int comp = fe.system_to_component_index(q).first;
                 if (fe_mask[comp])
-                  vector(dofs[q]) = points[q][fe_to_real[comp]];
+                  ::dealii::internal::ElementAccess<VectorType>::set(
+                    points[q][fe_to_real[comp]], dofs[q], vector);
               }
           }
       }
diff --git a/source/algorithms/operator.cc b/source/algorithms/operator.cc
index d71ae60074..884778c336 100644
--- a/source/algorithms/operator.cc
+++ b/source/algorithms/operator.cc
@@ -51,7 +51,6 @@ namespace Algorithms
     notifications.clear();
   }
 
-
 #include "operator.inst"
 }
 
diff --git a/source/base/time_stepping.cc b/source/base/time_stepping.cc
index 74c61291c5..372c4dce94 100644
--- a/source/base/time_stepping.cc
+++ b/source/base/time_stepping.cc
@@ -22,6 +22,7 @@
 #include <deal.II/lac/petsc_block_vector.h>
 #include <deal.II/lac/petsc_parallel_vector.h>
 #include <deal.II/lac/petsc_parallel_block_vector.h>
+#include <deal.II/lac/trilinos_epetra_vector.h>
 #include <deal.II/lac/trilinos_vector.h>
 #include <deal.II/lac/trilinos_block_vector.h>
 #include <deal.II/lac/trilinos_parallel_block_vector.h>
diff --git a/source/fe/fe_values.cc b/source/fe/fe_values.cc
index af64a397cb..a0ec6532a0 100644
--- a/source/fe/fe_values.cc
+++ b/source/fe/fe_values.cc
@@ -22,6 +22,7 @@
 #include <deal.II/lac/la_vector.h>
 #include <deal.II/lac/la_parallel_vector.h>
 #include <deal.II/lac/la_parallel_block_vector.h>
+#include <deal.II/lac/vector_element_access.h>
 #include <deal.II/lac/petsc_vector.h>
 #include <deal.II/lac/petsc_block_vector.h>
 #include <deal.II/lac/trilinos_vector.h>
@@ -47,7 +48,7 @@ namespace
   get_vector_element (const VectorType &vector,
                       const types::global_dof_index cell_number)
   {
-    return vector[cell_number];
+    return internal::ElementAccess<VectorType>::get(vector,cell_number);
   }
 
 
diff --git a/source/fe/mapping_fe_field.cc b/source/fe/mapping_fe_field.cc
index 88432c17b7..cf56a5d51c 100644
--- a/source/fe/mapping_fe_field.cc
+++ b/source/fe/mapping_fe_field.cc
@@ -1986,7 +1986,8 @@ MappingFEField<dim,spacedim,VectorType,DoFHandlerType>::update_internal_dofs
   dof_cell->get_dof_indices(data.local_dof_indices);
 
   for (unsigned int i=0; i<data.local_dof_values.size(); ++i)
-    data.local_dof_values[i] = (*euler_vector)(data.local_dof_indices[i]);
+    data.local_dof_values[i] = internal::ElementAccess<VectorType>::get(
+                                 *euler_vector, data.local_dof_indices[i]);
 }
 
 // explicit instantiations
diff --git a/source/lac/trilinos_sparse_matrix.cc b/source/lac/trilinos_sparse_matrix.cc
index e0622ba46a..f97722f38f 100644
--- a/source/lac/trilinos_sparse_matrix.cc
+++ b/source/lac/trilinos_sparse_matrix.cc
@@ -113,6 +113,57 @@ namespace TrilinosWrappers
 #endif
   }
 
+  namespace internal
+  {
+    template <typename VectorType>
+    double *begin(VectorType &V)
+    {
+      return V.begin();
+    }
+
+    template <>
+    double *begin(LinearAlgebra::EpetraWrappers::Vector &V)
+    {
+      return V.trilinos_vector()[0];
+    }
+
+    template <typename VectorType>
+    const double *begin(const VectorType &V)
+    {
+      return V.begin();
+    }
+
+    template <>
+    const double *begin(const LinearAlgebra::EpetraWrappers::Vector &V)
+    {
+      return V.trilinos_vector()[0];
+    }
+
+    template <typename VectorType>
+    double *end(VectorType &V)
+    {
+      return V.end();
+    }
+
+    template <>
+    double *end(LinearAlgebra::EpetraWrappers::Vector &V)
+    {
+      return V.trilinos_vector()[0]+V.trilinos_vector().MyLength();
+    }
+
+    template <typename VectorType>
+    const double *end(const VectorType &V)
+    {
+      return V.end();
+    }
+
+    template <>
+    const double *end(const LinearAlgebra::EpetraWrappers::Vector &V)
+    {
+      return V.trilinos_vector()[0]+V.trilinos_vector().MyLength();
+    }
+  }
+
 
   namespace SparseMatrixIterators
   {
@@ -1958,15 +2009,15 @@ namespace TrilinosWrappers
     (void)dst;
 
     internal::SparseMatrix::check_vector_map_equality(*matrix, src, dst);
-    const size_type dst_local_size = dst.end() - dst.begin();
+    const size_type dst_local_size = internal::end(dst) - internal::begin(dst);
     AssertDimension (dst_local_size, static_cast<size_type>(matrix->RangeMap().NumMyPoints()));
-    const size_type src_local_size = src.end() - src.begin();
+    const size_type src_local_size = internal::end(src) - internal::begin(src);
     AssertDimension (src_local_size, static_cast<size_type>(matrix->DomainMap().NumMyPoints()));
 
-    Epetra_MultiVector tril_dst (View, matrix->RangeMap(), dst.begin(),
+    Epetra_MultiVector tril_dst (View, matrix->RangeMap(), internal::begin(dst),
                                  dst_local_size, 1);
     Epetra_MultiVector tril_src (View, matrix->DomainMap(),
-                                 const_cast<TrilinosScalar *>(src.begin()),
+                                 const_cast<TrilinosScalar *>(internal::begin(src)),
                                  src_local_size, 1);
 
     const int ierr = matrix->Multiply (false, tril_src, tril_dst);
@@ -1985,15 +2036,15 @@ namespace TrilinosWrappers
     Assert (matrix->Filled(), ExcMatrixNotCompressed());
 
     internal::SparseMatrix::check_vector_map_equality(*matrix, dst, src);
-    const size_type dst_local_size = dst.end() - dst.begin();
+    const size_type dst_local_size = internal::end(dst) - internal::begin(dst);
     AssertDimension (dst_local_size, static_cast<size_type>(matrix->DomainMap().NumMyPoints()));
-    const size_type src_local_size = src.end() - src.begin();
+    const size_type src_local_size = internal::end(src) - internal::begin(src);
     AssertDimension (src_local_size, static_cast<size_type>(matrix->RangeMap().NumMyPoints()));
 
-    Epetra_MultiVector tril_dst (View, matrix->DomainMap(), dst.begin(),
+    Epetra_MultiVector tril_dst (View, matrix->DomainMap(), internal::begin(dst),
                                  dst_local_size, 1);
     Epetra_MultiVector tril_src (View, matrix->RangeMap(),
-                                 const_cast<double *>(src.begin()),
+                                 const_cast<double *>(internal::begin(src)),
                                  src_local_size, 1);
 
     const int ierr = matrix->Multiply (true, tril_src, tril_dst);
@@ -2016,9 +2067,11 @@ namespace TrilinosWrappers
     // deal.II local vector that has this fast setting. It will be accepted in
     // vmult because it only checks the local size.
     dealii::Vector<TrilinosScalar> temp_vector;
-    temp_vector.reinit(dst.end()-dst.begin(), true);
-    dealii::VectorView<TrilinosScalar> src_view(src.end()-src.begin(), src.begin());
-    dealii::VectorView<TrilinosScalar> dst_view(dst.end()-dst.begin(), dst.begin());
+    temp_vector.reinit(internal::end(dst)-internal::begin(dst), true);
+    dealii::VectorView<TrilinosScalar> src_view(internal::end(src)-internal::begin(src),
+                                                internal::begin(src));
+    dealii::VectorView<TrilinosScalar> dst_view(internal::end(dst)-internal::begin(dst),
+                                                internal::begin(dst));
     vmult (temp_vector, static_cast<const dealii::Vector<TrilinosScalar>&>(src_view));
     if (dst_view.size() > 0)
       dst_view += temp_vector;
@@ -2039,9 +2092,11 @@ namespace TrilinosWrappers
     // deal.II local vector that has this fast setting. It will be accepted in
     // vmult because it only checks the local size.
     dealii::Vector<TrilinosScalar> temp_vector;
-    temp_vector.reinit(dst.end()-dst.begin(), true);
-    dealii::VectorView<TrilinosScalar> src_view(src.end()-src.begin(), src.begin());
-    dealii::VectorView<TrilinosScalar> dst_view(dst.end()-dst.begin(), dst.begin());
+    temp_vector.reinit(internal::end(dst)-internal::begin(dst), true);
+    dealii::VectorView<TrilinosScalar> src_view(internal::end(src)-internal::begin(src),
+                                                internal::begin(src));
+    dealii::VectorView<TrilinosScalar> dst_view(internal::end(dst)-internal::begin(dst),
+                                                internal::begin(dst));
     Tvmult (temp_vector, static_cast<const dealii::Vector<TrilinosScalar>&>(src_view));
     if (dst_view.size() > 0)
       dst_view += temp_vector;
@@ -3406,6 +3461,9 @@ namespace TrilinosWrappers
   SparseMatrix::vmult (dealii::LinearAlgebra::distributed::Vector<double> &,
                        const dealii::LinearAlgebra::distributed::Vector<double> &) const;
   template void
+  SparseMatrix::vmult (dealii::LinearAlgebra::EpetraWrappers::Vector &,
+                       const dealii::LinearAlgebra::EpetraWrappers::Vector &) const;
+  template void
   SparseMatrix::Tvmult (VectorBase &,
                         const VectorBase &) const;
   template void
@@ -3421,6 +3479,9 @@ namespace TrilinosWrappers
   SparseMatrix::Tvmult (dealii::LinearAlgebra::distributed::Vector<double> &,
                         const dealii::LinearAlgebra::distributed::Vector<double> &) const;
   template void
+  SparseMatrix::Tvmult (dealii::LinearAlgebra::EpetraWrappers::Vector &,
+                        const dealii::LinearAlgebra::EpetraWrappers::Vector &) const;
+  template void
   SparseMatrix::vmult_add (VectorBase &,
                            const VectorBase &) const;
   template void
@@ -3436,6 +3497,9 @@ namespace TrilinosWrappers
   SparseMatrix::vmult_add (dealii::LinearAlgebra::distributed::Vector<double> &,
                            const dealii::LinearAlgebra::distributed::Vector<double> &) const;
   template void
+  SparseMatrix::vmult_add (dealii::LinearAlgebra::EpetraWrappers::Vector &,
+                           const dealii::LinearAlgebra::EpetraWrappers::Vector &) const;
+  template void
   SparseMatrix::Tvmult_add (VectorBase &,
                             const VectorBase &) const;
   template void
@@ -3450,6 +3514,9 @@ namespace TrilinosWrappers
   template void
   SparseMatrix::Tvmult_add (dealii::LinearAlgebra::distributed::Vector<double> &,
                             const dealii::LinearAlgebra::distributed::Vector<double> &) const;
+  template void
+  SparseMatrix::Tvmult_add (dealii::LinearAlgebra::EpetraWrappers::Vector &,
+                            const dealii::LinearAlgebra::EpetraWrappers::Vector &) const;
 
 }
 
diff --git a/source/multigrid/mg_level_global_transfer.cc b/source/multigrid/mg_level_global_transfer.cc
index 033a7ca73d..13b4ef7cdc 100644
--- a/source/multigrid/mg_level_global_transfer.cc
+++ b/source/multigrid/mg_level_global_transfer.cc
@@ -23,6 +23,7 @@
 #include <deal.II/lac/la_parallel_block_vector.h>
 #include <deal.II/lac/petsc_vector.h>
 #include <deal.II/lac/petsc_block_vector.h>
+#include <deal.II/lac/trilinos_epetra_vector.h>
 #include <deal.II/lac/trilinos_vector.h>
 #include <deal.II/lac/trilinos_block_vector.h>
 #include <deal.II/grid/tria.h>
diff --git a/source/multigrid/mg_transfer_prebuilt.cc b/source/multigrid/mg_transfer_prebuilt.cc
index 625890b882..d2c1883847 100644
--- a/source/multigrid/mg_transfer_prebuilt.cc
+++ b/source/multigrid/mg_transfer_prebuilt.cc
@@ -23,6 +23,7 @@
 #include <deal.II/lac/la_parallel_block_vector.h>
 #include <deal.II/lac/petsc_vector.h>
 #include <deal.II/lac/petsc_block_vector.h>
+#include <deal.II/lac/trilinos_epetra_vector.h>
 #include <deal.II/lac/trilinos_vector.h>
 #include <deal.II/lac/trilinos_block_vector.h>
 #include <deal.II/lac/sparse_matrix.h>
diff --git a/source/numerics/data_out_dof_data.cc b/source/numerics/data_out_dof_data.cc
index abac0f9b42..3b8ffea364 100644
--- a/source/numerics/data_out_dof_data.cc
+++ b/source/numerics/data_out_dof_data.cc
@@ -26,6 +26,7 @@
 #include <deal.II/lac/petsc_block_vector.h>
 #include <deal.II/lac/trilinos_vector.h>
 #include <deal.II/lac/trilinos_block_vector.h>
+#include <deal.II/lac/vector_element_access.h>
 #include <deal.II/numerics/data_out.h>
 #include <deal.II/numerics/data_out_dof_data.h>
 #include <deal.II/dofs/dof_handler.h>
@@ -522,7 +523,7 @@ namespace internal
       get_vector_element (const VectorType   &vector,
                           const unsigned int  cell_number)
       {
-        return vector[cell_number];
+        return internal::ElementAccess<VectorType>::get(vector,cell_number);
       }
 
 
diff --git a/source/numerics/point_value_history.cc b/source/numerics/point_value_history.cc
index d7b41f505f..8b339902da 100644
--- a/source/numerics/point_value_history.cc
+++ b/source/numerics/point_value_history.cc
@@ -23,6 +23,7 @@
 #include <deal.II/lac/petsc_block_vector.h>
 #include <deal.II/lac/trilinos_vector.h>
 #include <deal.II/lac/trilinos_block_vector.h>
+#include <deal.II/lac/vector_element_access.h>
 
 #include <deal.II/numerics/vector_tools.h>
 
@@ -589,7 +590,8 @@ void PointValueHistory<dim>
           if (mask->second[comp])
             {
               unsigned int solution_index = point->solution_indices[comp];
-              data_store_field->second[data_store_index * n_stored + store_index].push_back (solution (solution_index));
+              data_store_field->second[data_store_index * n_stored + store_index].push_back (
+                internal::ElementAccess<VectorType>::get(solution,solution_index));
               store_index++;
             }
         }
diff --git a/source/numerics/solution_transfer.cc b/source/numerics/solution_transfer.cc
index 7c27ab681b..a619124789 100644
--- a/source/numerics/solution_transfer.cc
+++ b/source/numerics/solution_transfer.cc
@@ -1,6 +1,6 @@
 // ---------------------------------------------------------------------
 //
-// Copyright (C) 1999 - 2016 by the deal.II authors
+// Copyright (C) 1999 - 2017 by the deal.II authors
 //
 // This file is part of the deal.II library.
 //
@@ -21,6 +21,7 @@
 #include <deal.II/dofs/dof_handler.h>
 #include <deal.II/dofs/dof_accessor.h>
 #include <deal.II/fe/fe.h>
+#include <deal.II/lac/vector_element_access.h>
 #include <deal.II/lac/vector.h>
 #include <deal.II/lac/la_vector.h>
 #include <deal.II/lac/la_parallel_vector.h>
@@ -34,8 +35,6 @@
 
 DEAL_II_NAMESPACE_OPEN
 
-
-
 template<int dim, typename VectorType, typename DoFHandlerType>
 SolutionTransfer<dim, VectorType, DoFHandlerType>::
 SolutionTransfer(const DoFHandlerType &dof)
@@ -158,7 +157,8 @@ SolutionTransfer<dim, VectorType, DoFHandlerType>::refine_interpolate
           Assert(dofs_per_cell==(*pointerstruct->second.indices_ptr).size(),
                  ExcInternalError());
           for (unsigned int i=0; i<dofs_per_cell; ++i)
-            local_values(i)=in((*pointerstruct->second.indices_ptr)[i]);
+            local_values(i)=internal::ElementAccess<VectorType>::get(
+                              in,(*pointerstruct->second.indices_ptr)[i]);
           cell->set_dof_values_by_interpolation(local_values, out,
                                                 this_fe_index);
         }
@@ -467,7 +467,8 @@ interpolate (const std::vector<VectorType> &all_in,
                 {
                   tmp.reinit (in_size, true);
                   for (unsigned int i=0; i<in_size; ++i)
-                    tmp(i) = all_in[j]((*indexptr)[i]);
+                    tmp(i) = internal::ElementAccess<VectorType>::get(all_in[j],
+                                                                      (*indexptr)[i]);
 
                   cell->set_dof_values_by_interpolation (tmp, all_out[j],
                                                          old_fe_index);
@@ -520,7 +521,8 @@ interpolate (const std::vector<VectorType> &all_in,
 
 
                   for (unsigned int i=0; i<dofs_per_cell; ++i)
-                    all_out[j](dofs[i])=(*data)(i);
+                    internal::ElementAccess<VectorType>::set((*data)(i), dofs[i],
+                                                             all_out[j]);
                 }
             }
           // undefined status
diff --git a/source/numerics/vector_tools_mean_value.cc b/source/numerics/vector_tools_mean_value.cc
index 43b8e51e98..890ff4d2ce 100644
--- a/source/numerics/vector_tools_mean_value.cc
+++ b/source/numerics/vector_tools_mean_value.cc
@@ -18,6 +18,19 @@
 
 DEAL_II_NAMESPACE_OPEN
 
+#ifdef DEAL_II_WITH_TRILINOS
+namespace VectorTools
+{
+  template <>
+  void
+  subtract_mean_value(LinearAlgebra::EpetraWrappers::Vector &,
+                      const std::vector<bool> &)
+  {
+    Assert(false,ExcNotImplemented());
+  }
+}
+#endif
+
 // ---------------------------- explicit instantiations --------------------
 #include "vector_tools_mean_value.inst"