From 6b3068af0a43d61b74710b0d87fbf9fc4dae8323 Mon Sep 17 00:00:00 2001 From: Bruno Turcksin Date: Fri, 10 Feb 2017 09:27:10 -0500 Subject: [PATCH] Instantiate LinearAlgebra::EpetraWrappers::Vector. --- cmake/config/template-arguments.in | 5 + cmake/configure/configure_2_trilinos.cmake | 1 + .../deal.II/algorithms/operator.templates.h | 7 +- include/deal.II/dofs/dof_accessor.templates.h | 70 +++++++++- .../fe/fe_tools_extrapolate.templates.h | 17 ++- .../fe/fe_tools_interpolate.templates.h | 20 ++- include/deal.II/lac/constraint_matrix.h | 9 +- .../deal.II/lac/constraint_matrix.templates.h | 13 +- include/deal.II/lac/trilinos_epetra_vector.h | 4 +- include/deal.II/lac/vector_element_access.h | 120 ++++++++++++++++++ include/deal.II/multigrid/mg_transfer.h | 17 +++ .../deal.II/multigrid/mg_transfer.templates.h | 1 + .../deal.II/numerics/vector_tools.templates.h | 42 ++++-- source/algorithms/operator.cc | 1 - source/base/time_stepping.cc | 1 + source/fe/fe_values.cc | 3 +- source/fe/mapping_fe_field.cc | 3 +- source/lac/trilinos_sparse_matrix.cc | 95 ++++++++++++-- source/multigrid/mg_level_global_transfer.cc | 1 + source/multigrid/mg_transfer_prebuilt.cc | 1 + source/numerics/data_out_dof_data.cc | 3 +- source/numerics/point_value_history.cc | 4 +- source/numerics/solution_transfer.cc | 14 +- source/numerics/vector_tools_mean_value.cc | 13 ++ 24 files changed, 403 insertions(+), 62 deletions(-) create mode 100644 include/deal.II/lac/vector_element_access.h 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; @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; @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; @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; @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 #include +#include DEAL_II_NAMESPACE_OPEN @@ -54,7 +55,8 @@ namespace Algorithms if (v == nullptr) continue; deallog << vectors.name(i); for (unsigned int j=0; jsize(); ++j) - deallog << ' ' << (*v)(j); + deallog << ' ' << ::dealii::internal::ElementAccess::get( + *v, j); deallog << std::endl; } deallog << std::endl; @@ -67,7 +69,8 @@ namespace Algorithms const VectorType *v = vectors.try_read_ptr(i); if (v == nullptr) continue; for (unsigned int j=0; jsize(); ++j) - (*os) << ' ' << (*v)(j); + (*os) << ' ' << ::dealii::internal::ElementAccess::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 #include #include +#include #include #include #include @@ -1479,6 +1480,66 @@ namespace internal Assert (next == dof_indices.end (), ExcInternalError ()); } + + + template + 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 sort_indices(const types::global_dof_index *v_begin, + const types::global_dof_index *v_end) + { + // initialize original index locations + std::vector 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 + 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 sorted_indices_pos = sort_indices(cache_begin, + cache_end); + const unsigned int cache_size = cache_end-cache_begin; + std::vector 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 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::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::set_dof_values ->get_cell_cache_start (this->present_index, this->get_fe().dofs_per_cell); for (unsigned int i=0; iget_fe().dofs_per_cell; ++i, ++cache) - values(*cache) = local_values(i); + internal::ElementAccess::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::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 + void reinit_distributed(const DoFHandler &dh, + LinearAlgebra::EpetraWrappers::Vector &vector) + { + const parallel::distributed::Triangulation *parallel_tria = + dynamic_cast*>(&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 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::add(u2_local(i), + dofs[i], u2); + ::dealii::internal::ElementAccess::add(1, + dofs[i], touch_count); } } } @@ -227,9 +229,16 @@ namespace FETools for (types::global_dof_index i=0; i(touch_count(i)) != typename OutVector::value_type(0), + Assert(static_cast( + ::dealii::internal::ElementAccess::get( + touch_count, i)) != typename OutVector::value_type(0), ExcInternalError()); - u2(i) /= touch_count(i); + + + const double val = ::dealii::internal::ElementAccess::get( + u2, i); + ::dealii::internal::ElementAccess::set( + val/::dealii::internal::ElementAccess::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::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 #include +#include #include #include @@ -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::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::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::set(0., idx, vec); } } @@ -895,10 +895,11 @@ ConstraintMatrix::distribute (VectorType &vec) const new_value = it->inhomogeneity; for (unsigned int i=0; ientries.size(); ++i) new_value += (static_cast - (ghosted_vector(it->entries[i].first)) * + (internal::ElementAccess::get( + ghosted_vector, it->entries[i].first)) * it->entries[i].second); AssertIsFinite(new_value); - vec(it->line) = new_value; + internal::ElementAccess::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; ientries.size(); ++i) new_value += (static_cast - (vec(next_constraint->entries[i].first)) * + (internal::ElementAccess::get( + vec, next_constraint->entries[i].first))* next_constraint->entries[i].second); AssertIsFinite(new_value); - vec(next_constraint->line) = new_value; + internal::ElementAccess::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 &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 : std_cxx11::false_type +struct is_serial_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_NAMESPACE_OPEN + +namespace internal +{ + template + 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 + inline + void ElementAccess::add(typename VectorType::value_type value, + types::global_dof_index i, VectorType &V) + { + V(i) += value; + } + + + + template + inline + void ElementAccess::set(typename VectorType::value_type value, + types::global_dof_index i, VectorType &V) + { + V(i) = value; + } + + + + template + inline + typename VectorType::value_type + ElementAccess::get(const VectorType &V, types::global_dof_index i) + { + return V(i); + } + + + +#ifdef DEAL_II_WITH_TRILINOS + template <> + inline + void ElementAccess:: + 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(i)); + + vector[0][trilinos_i] += value; + } + + + + template <> + inline + void ElementAccess:: + 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(i)); + + vector[0][trilinos_i] = value; + } + + + template <> + inline + double + ElementAccess:: + 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(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 { @@ -88,6 +89,22 @@ namespace internal }; + template <> + struct MatrixSelector + { + typedef ::dealii::TrilinosWrappers::SparsityPattern Sparsity; + typedef ::dealii::TrilinosWrappers::SparseMatrix Matrix; + + template + 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 { 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 +#include #include #include #include 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 #include #include +#include #include #include #include @@ -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::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::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::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::get( + data_2, i); + + ::dealii::internal::ElementAccess::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::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::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 and that wouldn't // necessarily go together with the matrix and other functions for (unsigned int i=0; i::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::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::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::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::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 #include #include +#include #include #include #include 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 #include #include +#include #include #include #include @@ -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::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::update_internal_dofs dof_cell->get_dof_indices(data.local_dof_indices); for (unsigned int i=0; i::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 + double *begin(VectorType &V) + { + return V.begin(); + } + + template <> + double *begin(LinearAlgebra::EpetraWrappers::Vector &V) + { + return V.trilinos_vector()[0]; + } + + template + const double *begin(const VectorType &V) + { + return V.begin(); + } + + template <> + const double *begin(const LinearAlgebra::EpetraWrappers::Vector &V) + { + return V.trilinos_vector()[0]; + } + + template + double *end(VectorType &V) + { + return V.end(); + } + + template <> + double *end(LinearAlgebra::EpetraWrappers::Vector &V) + { + return V.trilinos_vector()[0]+V.trilinos_vector().MyLength(); + } + + template + 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(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(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(src.begin()), + const_cast(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(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(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(src.begin()), + const_cast(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 temp_vector; - temp_vector.reinit(dst.end()-dst.begin(), true); - dealii::VectorView src_view(src.end()-src.begin(), src.begin()); - dealii::VectorView dst_view(dst.end()-dst.begin(), dst.begin()); + temp_vector.reinit(internal::end(dst)-internal::begin(dst), true); + dealii::VectorView src_view(internal::end(src)-internal::begin(src), + internal::begin(src)); + dealii::VectorView dst_view(internal::end(dst)-internal::begin(dst), + internal::begin(dst)); vmult (temp_vector, static_cast&>(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 temp_vector; - temp_vector.reinit(dst.end()-dst.begin(), true); - dealii::VectorView src_view(src.end()-src.begin(), src.begin()); - dealii::VectorView dst_view(dst.end()-dst.begin(), dst.begin()); + temp_vector.reinit(internal::end(dst)-internal::begin(dst), true); + dealii::VectorView src_view(internal::end(src)-internal::begin(src), + internal::begin(src)); + dealii::VectorView dst_view(internal::end(dst)-internal::begin(dst), + internal::begin(dst)); Tvmult (temp_vector, static_cast&>(src_view)); if (dst_view.size() > 0) dst_view += temp_vector; @@ -3406,6 +3461,9 @@ namespace TrilinosWrappers SparseMatrix::vmult (dealii::LinearAlgebra::distributed::Vector &, const dealii::LinearAlgebra::distributed::Vector &) 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 &, const dealii::LinearAlgebra::distributed::Vector &) 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 &, const dealii::LinearAlgebra::distributed::Vector &) 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 &, const dealii::LinearAlgebra::distributed::Vector &) 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 #include #include +#include #include #include #include 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 #include #include +#include #include #include #include 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 #include #include +#include #include #include #include @@ -522,7 +523,7 @@ namespace internal get_vector_element (const VectorType &vector, const unsigned int cell_number) { - return vector[cell_number]; + return internal::ElementAccess::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 #include #include +#include #include @@ -589,7 +590,8 @@ void PointValueHistory 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::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 #include #include +#include #include #include #include @@ -34,8 +35,6 @@ DEAL_II_NAMESPACE_OPEN - - template SolutionTransfer:: SolutionTransfer(const DoFHandlerType &dof) @@ -158,7 +157,8 @@ SolutionTransfer::refine_interpolate Assert(dofs_per_cell==(*pointerstruct->second.indices_ptr).size(), ExcInternalError()); for (unsigned int i=0; isecond.indices_ptr)[i]); + local_values(i)=internal::ElementAccess::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 &all_in, { tmp.reinit (in_size, true); for (unsigned int i=0; i::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 &all_in, for (unsigned int i=0; i::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 &) + { + Assert(false,ExcNotImplemented()); + } +} +#endif + // ---------------------------- explicit instantiations -------------------- #include "vector_tools_mean_value.inst" -- 2.39.5