From: bangerth Date: Thu, 21 Aug 2008 22:33:41 +0000 (+0000) Subject: Make available a whole bunch of functions in deal.II that work on Trilinos vectors... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=13dfedfd632cdfc7a9e7b8e3d6a3143ecfd579a6;p=dealii-svn.git Make available a whole bunch of functions in deal.II that work on Trilinos vectors and matrices. git-svn-id: https://svn.dealii.org/trunk@16650 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/dofs/dof_accessor.h b/deal.II/deal.II/include/dofs/dof_accessor.h index 0f98d88285..8c9cc19366 100644 --- a/deal.II/deal.II/include/dofs/dof_accessor.h +++ b/deal.II/deal.II/include/dofs/dof_accessor.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -1072,16 +1072,19 @@ class DoFCellAccessor : public DoFObjectAccessor * function is only callable for * active cells. * - * The input vector may be either a - * Vector, + * The input vector may be either + * a Vector, * Vector, or a * BlockVector, or a - * PETSc vector if deal.II is compiled to - * support PETSc. It is in the - * responsibility of the caller to assure - * that the types of the numbers stored - * in input and output vectors are - * compatible and with similar accuracy. + * PETSc or Trilinos vector if + * deal.II is compiled to support + * these libraries. It is in the + * responsibility of the caller + * to assure that the types of + * the numbers stored in input + * and output vectors are + * compatible and with similar + * accuracy. */ template void get_dof_values (const InputVector &values, @@ -1115,7 +1118,7 @@ class DoFCellAccessor : public DoFObjectAccessor * Vector, or a * BlockVector, or a * PETSc vector if deal.II is compiled to - * support PETSc. It is in the + * support these libraries. It is in the * responsibility of the caller to assure * that the types of the numbers stored * in input and output vectors are @@ -1244,7 +1247,7 @@ class DoFCellAccessor : public DoFObjectAccessor * Vector, or a * BlockVector, or a * PETSc vector if deal.II is compiled to - * support PETSc. It is in the + * support these libraries. It is in the * responsibility of the caller to assure * that the types of the numbers stored * in input and output vectors are diff --git a/deal.II/deal.II/include/dofs/dof_constraints.h b/deal.II/deal.II/include/dofs/dof_constraints.h index 153f8a1068..d04118ec41 100644 --- a/deal.II/deal.II/include/dofs/dof_constraints.h +++ b/deal.II/deal.II/include/dofs/dof_constraints.h @@ -178,17 +178,18 @@ class BlockIndices; * arguments refer to the first possibility above, those taking only one do * their job in-place and refer to the second possibility. * - * The condensation functions exist for different argument types. The in-place - * functions (i.e. those following the second way) exist for arguments of type - * SparsityPattern, SparseMatrix and BlockSparseMatrix. Note that there are no - * versions for arguments of type PETScWrappers::SparseMatrix() or any of the - * other PETSc matrix wrapper classes. This is due to the fact that it is - * relatively hard to get a representation of the sparsity structure of PETSc - * matrices, and to modify them; this holds in particular, if the matrix is - * actually distributed across a cluster of computers. If you want to use - * PETSc matrices, you can either copy an already condensed deal.II matrix, or - * build the PETSc matrix in the already condensed form, see the discussion - * below. + * The condensation functions exist for different argument types. The + * in-place functions (i.e. those following the second way) exist for + * arguments of type SparsityPattern, SparseMatrix and + * BlockSparseMatrix. Note that there are no versions for arguments of + * type PETScWrappers::SparseMatrix() or any of the other PETSc or + * Trilinos matrix wrapper classes. This is due to the fact that it is + * relatively hard to get a representation of the sparsity structure + * of PETSc matrices, and to modify them; this holds in particular, if + * the matrix is actually distributed across a cluster of + * computers. If you want to use PETSc matrices, you can either copy + * an already condensed deal.II matrix, or build the PETSc matrix in + * the already condensed form, see the discussion below. * * *
Condensing vectors
@@ -199,10 +200,11 @@ class BlockIndices; * object has been condensed, further condensation operations don't change it * any more. * - * In contrast to the matrix condensation functions, the vector condensation - * functions exist in variants for PETSc vectors. However, using them is - * typically expensive, and should be avoided. You should use the same - * techniques as mentioned above to avoid their use. + * In contrast to the matrix condensation functions, the vector + * condensation functions exist in variants for PETSc and Trilinos + * vectors. However, using them is typically expensive, and should be + * avoided. You should use the same techniques as mentioned above to + * avoid their use. * * *

Avoiding explicit condensation

@@ -219,12 +221,13 @@ class BlockIndices; * paper". This is the case discussed in the hp tutorial program, * @ref step_27 "step-27", as well as in @ref step_31 "step-31". * - *
  • There may not be a condense() - * function for the matrix you use (this is, for example, the case - * for the PETSc wrapper classes, where we have no access to the underlying - * representation of the matrix, and therefore cannot efficiently implement - * the condense() operation). This is the case discussed in - * @ref step_17 "step-17" and @ref step_18 "step-18". + *
  • + * There may not be a condense() function for the matrix you use + * (this is, for example, the case for the PETSc and Trilinos wrapper + * classes, where we have no access to the underlying representation + * of the matrix, and therefore cannot efficiently implement the + * condense() operation). This is the case discussed in @ref step_17 + * "step-17" and @ref step_18 "step-18". * * * In this case, one possibility is to distribute local entries to the final @@ -790,11 +793,12 @@ class ConstraintMatrix : public Subscriptor * @p condensed be zero. * * The @p VectorType may be a - * Vector, - * Vector, - * BlockVector<...>, a PETSc - * vector wrapper class, or any other - * type having the same interface. + * Vector, Vector, + * BlockVector<...>, a + * PETSc or Trilinos vector + * wrapper class, or any other + * type having the same + * interface. */ template void condense (const VectorType &uncondensed, @@ -802,12 +806,14 @@ class ConstraintMatrix : public Subscriptor /** * Condense the given vector - * in-place. The @p VectorType may be a - * Vector, + * in-place. The @p VectorType + * may be a Vector, * Vector, - * BlockVector<...>, a PETSc - * vector wrapper class, or any other - * type having the same interface. + * BlockVector<...>, a + * PETSc or Trilinos vector + * wrapper class, or any other + * type having the same + * interface. */ template void condense (VectorType &vec) const; @@ -1016,14 +1022,16 @@ class ConstraintMatrix : public Subscriptor const bool keep_constrained_entries = true) const; /** - * Delete hanging nodes in a vector. - * Sets all hanging node values to - * zero. The @p VectorType may be a - * Vector, - * Vector, - * BlockVector<...>, a PETSc - * vector wrapper class, or any other - * type having the same interface. + * Delete hanging nodes in a + * vector. Sets all hanging node + * values to zero. The @p + * VectorType may be a + * Vector, Vector, + * BlockVector<...>, a + * PETSc or Trilinos vector + * wrapper class, or any other + * type having the same + * interface. */ template void set_zero (VectorType &vec) const; @@ -1053,24 +1061,27 @@ class ConstraintMatrix : public Subscriptor * @p condense. * * The @p VectorType may be a - * Vector, - * Vector, - * BlockVector<...>, a PETSc - * vector wrapper class, or any other - * type having the same interface. + * Vector, Vector, + * BlockVector<...>, a + * PETSc or Trilinos vector + * wrapper class, or any other + * type having the same + * interface. */ template void distribute (const VectorType &condensed, VectorType &uncondensed) const; /** - * Re-distribute the elements of the - * vector in-place. The @p VectorType - * may be a Vector, - * Vector, - * BlockVector<...>, a PETSc - * vector wrapper class, or any other - * type having the same interface. + * Re-distribute the elements of + * the vector in-place. The @p + * VectorType may be a + * Vector, Vector, + * BlockVector<...>, a + * PETSc or Trilinos vector + * wrapper class, or any other + * type having the same + * interface. */ template void distribute (VectorType &vec) const; diff --git a/deal.II/deal.II/include/fe/fe_values.h b/deal.II/deal.II/include/fe/fe_values.h index cc843b40ca..c5e0dd5d09 100644 --- a/deal.II/deal.II/include/fe/fe_values.h +++ b/deal.II/deal.II/include/fe/fe_values.h @@ -27,6 +27,7 @@ #include #include #include +#include #include #include #include @@ -1125,15 +1126,18 @@ class FEValuesBase : protected FEValuesData, * values object already has the * correct size. * - * The actual data type of the input - * vector may be either a Vector<T>, - * BlockVector<T>, or one of the - * sequential PETSc vector wrapper classes. It - * represents a global vector of - * DoF values associated with the - * DofHandler object with - * which this FEValues - * object was last initialized. + * The actual data type of the + * input vector may be either a + * Vector<T>, + * BlockVector<T>, or one + * of the sequential PETSc or + * Trilinos vector wrapper + * classes. It represents a + * global vector of DoF values + * associated with the DofHandler + * object with which this + * FEValues object was last + * initialized. */ template void get_function_values (const InputVector& fe_function, @@ -1149,15 +1153,18 @@ class FEValuesBase : protected FEValuesData, * applied to multi-component * elements. * - * The actual data type of the input - * vector may be either a Vector<T>, - * BlockVector<T>, or one of the - * sequential PETSc vector wrapper classes. It - * represents a global vector of - * DoF values associated with the - * DofHandler object with - * which this FEValues - * object was last initialized. + * The actual data type of the + * input vector may be either a + * Vector<T>, + * BlockVector<T>, or one + * of the sequential PETSc or + * Trilinos vector wrapper + * classes. It represents a + * global vector of DoF values + * associated with the DofHandler + * object with which this + * FEValues object was last + * initialized. */ template void get_function_values (const InputVector &fe_function, @@ -1328,15 +1335,18 @@ class FEValuesBase : protected FEValuesData, * @p gradients object already has the * right size. * - * The actual data type of the input - * vector may be either a Vector<T>, - * BlockVector<T>, or one of the - * sequential PETSc vector wrapper classes. It - * represents a global vector of - * DoF values associated with the - * DofHandler object with - * which this FEValues - * object was last initialized. + * The actual data type of the + * input vector may be either a + * Vector<T>, + * BlockVector<T>, or one + * of the sequential PETSc or + * Trilinos vector wrapper + * classes. It represents a + * global vector of DoF values + * associated with the DofHandler + * object with which this + * FEValues object was last + * initialized. * * The output are the gradients * of the function represented by @@ -1375,15 +1385,18 @@ class FEValuesBase : protected FEValuesData, * but applied to multi-component * elements. * - * The actual data type of the input - * vector may be either a Vector<T>, - * BlockVector<T>, or one of the - * sequential PETSc vector wrapper classes. It - * represents a global vector of - * DoF values associated with the - * DofHandler object with - * which this FEValues - * object was last initialized. + * The actual data type of the + * input vector may be either a + * Vector<T>, + * BlockVector<T>, or one + * of the sequential PETSc or + * Trilinos vector wrapper + * classes. It represents a + * global vector of DoF values + * associated with the DofHandler + * object with which this + * FEValues object was last + * initialized. * * The output are the gradients * of the function represented by @@ -1466,15 +1479,18 @@ class FEValuesBase : protected FEValuesData, * get_function_hessians() * function. * - * The actual data type of the input - * vector may be either a Vector<T>, - * BlockVector<T>, or one of the - * sequential PETSc vector wrapper classes..It - * represents a global vector of - * DoF values associated with the - * DofHandler object with - * which this FEValues - * object was last initialized. + * The actual data type of the + * input vector may be either a + * Vector<T>, + * BlockVector<T>, or one + * of the sequential PETSc or + * Trilinos vector wrapper + * classes..It represents a + * global vector of DoF values + * associated with the DofHandler + * object with which this + * FEValues object was last + * initialized. * * The output are the second derivatives * of the function represented by @@ -1504,15 +1520,18 @@ class FEValuesBase : protected FEValuesData, * name, but applies to * vector-valued finite elements. * - * The actual data type of the input - * vector may be either a Vector<T>, - * BlockVector<T>, or one of the - * sequential PETSc vector wrapper classes. It - * represents a global vector of - * DoF values associated with the - * DofHandler object with - * which this FEValues - * object was last initialized. + * The actual data type of the + * input vector may be either a + * Vector<T>, + * BlockVector<T>, or one + * of the sequential PETSc or + * Trilinos vector wrapper + * classes. It represents a + * global vector of DoF values + * associated with the DofHandler + * object with which this + * FEValues object was last + * initialized. * * The output are the second derivatives * of the function represented by @@ -1894,6 +1913,30 @@ class FEValuesBase : protected FEValuesData, void get_interpolated_dof_values (const PETScWrappers::BlockVector &in, Vector &out) const = 0; +#endif + +#ifdef DEAL_II_USE_TRILINOS + /** + * Call + * @p get_interpolated_dof_values + * of the iterator with the + * given arguments. + */ + virtual + void + get_interpolated_dof_values (const TrilinosWrappers::Vector &in, + Vector &out) const = 0; + +// /** +// * Call +// * @p get_interpolated_dof_values +// * of the iterator with the +// * given arguments. +// */ +// virtual +// void +// get_interpolated_dof_values (const TrilinosWrappers::BlockVector &in, +// Vector &out) const = 0; #endif }; @@ -2013,6 +2056,31 @@ class FEValuesBase : protected FEValuesData, get_interpolated_dof_values (const PETScWrappers::BlockVector &in, Vector &out) const; #endif + +#ifdef DEAL_II_USE_TRILINOS + /** + * Call + * @p get_interpolated_dof_values + * of the iterator with the + * given arguments. + */ + virtual + void + get_interpolated_dof_values (const TrilinosWrappers::Vector &in, + Vector &out) const; + +// /** +// * Call +// * @p get_interpolated_dof_values +// * of the iterator with the +// * given arguments. +// */ +// virtual +// void +// get_interpolated_dof_values (const TrilinosWrappers::BlockVector &in, +// Vector &out) const; +#endif + private: /** * Copy of the iterator which @@ -2175,6 +2243,31 @@ class FEValuesBase : protected FEValuesData, get_interpolated_dof_values (const PETScWrappers::BlockVector &in, Vector &out) const; #endif + +#ifdef DEAL_II_USE_TRILINOS + /** + * Call + * @p get_interpolated_dof_values + * of the iterator with the + * given arguments. + */ + virtual + void + get_interpolated_dof_values (const TrilinosWrappers::Vector &in, + Vector &out) const; + +// /** +// * Call +// * @p get_interpolated_dof_values +// * of the iterator with the +// * given arguments. +// */ +// virtual +// void +// get_interpolated_dof_values (const TrilinosWrappers::BlockVector &in, +// Vector &out) const; +#endif + private: /** * Copy of the iterator which diff --git a/deal.II/deal.II/include/numerics/matrices.h b/deal.II/deal.II/include/numerics/matrices.h index b4c05d17c7..476afa305d 100644 --- a/deal.II/deal.II/include/numerics/matrices.h +++ b/deal.II/deal.II/include/numerics/matrices.h @@ -60,6 +60,14 @@ namespace PETScWrappers } #endif +#ifdef DEAL_II_USE_TRILINOS +namespace TrilinosWrappers +{ + class SparseMatrix; + class Vector; +} +#endif + /** * This class provides functions that assemble certain standard matrices for a @@ -1096,6 +1104,52 @@ class MatrixTools : public MatrixCreator const bool eliminate_columns = true); #endif +#ifdef DEAL_II_USE_TRILINOS + /** + * Apply dirichlet boundary conditions to + * the system matrix and vectors as + * described in the general + * documentation. This function works on + * the classes that are used to wrap + * Trilinos objects. + * + * Note that this function is not very + * efficient: it needs to alternatingly + * read and write into the matrix, a + * situation that Trilinos does not handle + * too well. In addition, we only get rid + * of rows corresponding to boundary + * nodes, but the corresponding case of + * deleting the respective columns + * (i.e. if @p eliminate_columns is @p + * true) is not presently implemented, + * and probably will never because it is + * too expensive without direct access to + * the Trilinos data structures. (This leads + * to the situation where the action + * indicates by the default value of the + * last argument is actually not + * implemented; that argument has + * true as its default value + * to stay consistent with the other + * functions of same name in this class.) + * A third reason against this function + * is that it doesn't handle the case + * where the matrix is distributed across + * an MPI system. + * + * This function is used in + * @ref step_17 "step-17" and + * @ref step_18 "step-18". + */ + static void + apply_boundary_values (const std::map &boundary_values, + TrilinosWrappers::SparseMatrix &matrix, + TrilinosWrappers::Vector &solution, + TrilinosWrappers::Vector &right_hand_side, + const bool eliminate_columns = true); +#endif + /** * Rather than applying boundary * values to the global matrix @@ -1124,6 +1178,7 @@ class MatrixTools : public MatrixCreator * access to the sparse matrix is * expensive (e.g. for block * sparse matrices, or for PETSc + * or trilinos * matrices). However, it doesn't * work as expected if there are * also hanging nodes to be diff --git a/deal.II/deal.II/source/dofs/dof_accessor.cc b/deal.II/deal.II/source/dofs/dof_accessor.cc index ab2ad3ed28..29b11f86cd 100644 --- a/deal.II/deal.II/source/dofs/dof_accessor.cc +++ b/deal.II/deal.II/source/dofs/dof_accessor.cc @@ -16,6 +16,7 @@ #include #include #include +#include #include #include @@ -646,6 +647,45 @@ DoFCellAccessor >::set_dof_values >::get_dof_values +(const TrilinosWrappers::Vector &, Vector&) const; +template +void +DoFCellAccessor >::get_dof_values +(const TrilinosWrappers::Vector &, Vector&) const; +template +void +DoFCellAccessor >::set_dof_values +(const Vector &, TrilinosWrappers::Vector&) const; +template +void +DoFCellAccessor >::set_dof_values +(const Vector&, TrilinosWrappers::Vector&) const; + +// template +// void +// DoFCellAccessor >::get_dof_values +// (const TrilinosWrappers::BlockVector &, Vector&) const; +// template +// void +// DoFCellAccessor >::get_dof_values +// (const TrilinosWrappers::BlockVector &, Vector&) const; +// template +// void +// DoFCellAccessor >::set_dof_values +// (const Vector &, TrilinosWrappers::BlockVector&) const; +// template +// void +// DoFCellAccessor >::set_dof_values +// (const Vector&, TrilinosWrappers::BlockVector&) const; + +#endif + + + template void DoFCellAccessor >:: @@ -807,6 +847,58 @@ set_dof_values_by_interpolation #endif +// for Trilinos vectors +#if DEAL_II_USE_TRILINOS +template +void +DoFCellAccessor >:: +get_interpolated_dof_values +(const TrilinosWrappers::Vector&, Vector&) const; +template +void +DoFCellAccessor >:: +set_dof_values_by_interpolation +(const Vector&, TrilinosWrappers::Vector&) const; + +template +void +DoFCellAccessor >:: +get_interpolated_dof_values +(const TrilinosWrappers::Vector&, Vector&) const; +template +void +DoFCellAccessor >:: +set_dof_values_by_interpolation +(const Vector&, TrilinosWrappers::Vector&) const; + + +// template +// void +// DoFCellAccessor >:: +// get_interpolated_dof_values +// (const TrilinosWrappers::BlockVector&, Vector&) const; +// template +// void +// DoFCellAccessor >:: +// set_dof_values_by_interpolation +// (const Vector&, TrilinosWrappers::BlockVector&) const; +// template + +// template +// void +// DoFCellAccessor >:: +// get_interpolated_dof_values +// (const TrilinosWrappers::BlockVector&, Vector&) const; +// template +// void +// DoFCellAccessor >:: +// set_dof_values_by_interpolation +// (const Vector&, TrilinosWrappers::BlockVector&) const; +// template + +#endif + + #if deal_II_dimension == 1 template class DoFAccessor<1, DoFHandler<1> >; template class DoFObjectAccessor<1, DoFHandler<1> >; @@ -957,6 +1049,44 @@ DoFCellAccessor >::set_dof_values&, PETScWrappers::BlockVector&) const; #endif +// for Trilinos vectors +#ifdef DEAL_II_USE_TRILINOS +template +void +DoFCellAccessor >::get_dof_values +(const TrilinosWrappers::Vector &, Vector&) const; +template +void +DoFCellAccessor >::get_dof_values +(const TrilinosWrappers::Vector &, Vector&) const; +template +void +DoFCellAccessor >::set_dof_values +(const Vector &, TrilinosWrappers::Vector&) const; +template +void +DoFCellAccessor >::set_dof_values +(const Vector&, TrilinosWrappers::Vector&) const; + +// template +// void +// DoFCellAccessor >::get_dof_values +// (const TrilinosWrappers::BlockVector &, Vector&) const; +// template +// void +// DoFCellAccessor >::get_dof_values +// (const TrilinosWrappers::BlockVector &, Vector&) const; +// template +// void +// DoFCellAccessor >::set_dof_values +// (const Vector &, TrilinosWrappers::BlockVector&) const; +// template +// void +// DoFCellAccessor >::set_dof_values +// (const Vector&, TrilinosWrappers::BlockVector&) const; +#endif + + template void @@ -1098,6 +1228,55 @@ set_dof_values_by_interpolation #endif +// for Trilinos vectors +#ifdef DEAL_II_USE_TRILINOS +template +void +DoFCellAccessor >:: +get_interpolated_dof_values +(const TrilinosWrappers::Vector&, Vector&) const; +template +void +DoFCellAccessor >:: +set_dof_values_by_interpolation +(const Vector&, TrilinosWrappers::Vector&) const; + +template +void +DoFCellAccessor >:: +get_interpolated_dof_values +(const TrilinosWrappers::Vector&, Vector&) const; +template +void +DoFCellAccessor >:: +set_dof_values_by_interpolation +(const Vector&, TrilinosWrappers::Vector&) const; + + +// template +// void +// DoFCellAccessor >:: +// get_interpolated_dof_values +// (const TrilinosWrappers::BlockVector&, Vector&) const; +// template +// void +// DoFCellAccessor >:: +// set_dof_values_by_interpolation +// (const Vector&, TrilinosWrappers::BlockVector&) const; + +// template +// void +// DoFCellAccessor >:: +// get_interpolated_dof_values +// (const TrilinosWrappers::BlockVector&, Vector&) const; +// template +// void +// DoFCellAccessor >:: +// set_dof_values_by_interpolation +// (const Vector&, TrilinosWrappers::BlockVector&) const; +#endif + + #if deal_II_dimension == 1 template class DoFAccessor<1, hp::DoFHandler<1> >; template class DoFObjectAccessor<1, hp::DoFHandler<1> >; diff --git a/deal.II/deal.II/source/dofs/dof_constraints.cc b/deal.II/deal.II/source/dofs/dof_constraints.cc index c0b7e1a09a..fe9d36ba8a 100644 --- a/deal.II/deal.II/source/dofs/dof_constraints.cc +++ b/deal.II/deal.II/source/dofs/dof_constraints.cc @@ -31,6 +31,8 @@ #include #include #include +#include +#include #include #include @@ -1733,6 +1735,11 @@ VECTOR_FUNCTIONS(PETScWrappers::MPI::Vector); VECTOR_FUNCTIONS(PETScWrappers::MPI::BlockVector); #endif +#ifdef DEAL_II_USE_TRILINOS +VECTOR_FUNCTIONS(TrilinosWrappers::Vector); +//VECTOR_FUNCTIONS(TrilinosWrappers::BlockVector); +#endif + template void @@ -1782,6 +1789,11 @@ MATRIX_FUNCTIONS(PETScWrappers::MPI::SparseMatrix); MATRIX_FUNCTIONS(PETScWrappers::MPI::BlockSparseMatrix); #endif +#ifdef DEAL_II_USE_TRILINOS +MATRIX_FUNCTIONS(TrilinosWrappers::SparseMatrix); +//MATRIX_FUNCTIONS(TrilinosWrappers::BlockSparseMatrix); +#endif + template void ConstraintMatrix:: add_entries_local_to_global (const std::vector &, SparsityPattern &, diff --git a/deal.II/deal.II/source/fe/fe_values.cc b/deal.II/deal.II/source/fe/fe_values.cc index 1c5de4a665..34122c3c49 100644 --- a/deal.II/deal.II/source/fe/fe_values.cc +++ b/deal.II/deal.II/source/fe/fe_values.cc @@ -18,13 +18,11 @@ #include #include #include +#include #include #include #include #include -#include -#include -#include #include #include #include @@ -144,6 +142,32 @@ get_interpolated_dof_values (const PETScWrappers::BlockVector &in, #endif +#ifdef DEAL_II_USE_TRILINOS + +template +template +void +FEValuesBase::CellIterator:: +get_interpolated_dof_values (const TrilinosWrappers::Vector &in, + Vector &out) const +{ + cell->get_interpolated_dof_values (in, out); +} + + + +// template +// template +// void +// FEValuesBase::CellIterator:: +// get_interpolated_dof_values (const TrilinosWrappers::BlockVector &in, +// Vector &out) const +// { +// cell->get_interpolated_dof_values (in, out); +// } + +#endif + /* ---------------- FEValuesBase::TriaCellIterator --------- */ @@ -252,6 +276,30 @@ get_interpolated_dof_values (const PETScWrappers::BlockVector &, #endif +#ifdef DEAL_II_USE_TRILINOS + +template +void +FEValuesBase::TriaCellIterator:: +get_interpolated_dof_values (const TrilinosWrappers::Vector &, + Vector &) const +{ + Assert (false, ExcMessage (message_string)); +} + + + +// template +// void +// FEValuesBase::TriaCellIterator:: +// get_interpolated_dof_values (const TrilinosWrappers::BlockVector &, +// Vector &) const +// { +// Assert (false, ExcMessage (message_string)); +// } + +#endif + @@ -2057,3 +2105,13 @@ DEAL_II_NAMESPACE_CLOSE #undef IN #endif +#ifdef DEAL_II_USE_TRILINOS +#define IN TrilinosWrappers::Vector +#include "fe_values.instance.h" +#undef IN + +// #define IN TrilinosWrappers::BlockVector +// #include "fe_values.instance.h" +// #undef IN +#endif + diff --git a/deal.II/deal.II/source/numerics/matrices.all_dimensions.cc b/deal.II/deal.II/source/numerics/matrices.all_dimensions.cc index 654b871f58..78da534924 100644 --- a/deal.II/deal.II/source/numerics/matrices.all_dimensions.cc +++ b/deal.II/deal.II/source/numerics/matrices.all_dimensions.cc @@ -26,6 +26,11 @@ # include #endif +#ifdef DEAL_II_USE_TRILINOS +# include +# include +#endif + #include DEAL_II_NAMESPACE_OPEN @@ -677,6 +682,148 @@ apply_boundary_values (const std::map &boundary_values, #endif + +#ifdef DEAL_II_USE_TRILINOS + +namespace TrilinosWrappers +{ + template + void + apply_boundary_values (const std::map &boundary_values, + TrilinosMatrix &matrix, + TrilinosVector &solution, + TrilinosVector &right_hand_side, + const bool eliminate_columns) + { + Assert (eliminate_columns == false, ExcNotImplemented()); + + Assert (matrix.n() == right_hand_side.size(), + ExcDimensionMismatch(matrix.n(), right_hand_side.size())); + Assert (matrix.n() == solution.size(), + ExcDimensionMismatch(matrix.n(), solution.size())); + + // if no boundary values are to be applied + // simply return + if (boundary_values.size() == 0) + return; + + const std::pair local_range + = matrix.local_range(); + Assert (local_range == right_hand_side.local_range(), + ExcInternalError()); + Assert (local_range == solution.local_range(), + ExcInternalError()); + + + // we have to read and write from this + // matrix (in this order). this will only + // work if we compress the matrix first, + // done here + matrix.compress (); + + // determine the first nonzero diagonal + // entry from within the part of the + // matrix that we can see. if we can't + // find such an entry, take one + TrilinosScalar average_nonzero_diagonal_entry = 1; + for (unsigned int i=local_range.first; i constrained_rows; + for (std::map::const_iterator + dof = boundary_values.begin(); + dof != boundary_values.end(); + ++dof) + if ((dof->first >= local_range.first) && + (dof->first < local_range.second)) + constrained_rows.push_back (dof->first); + + // then eliminate these rows and set + // their diagonal entry to what we have + // determined above. note that for trilinos + // matrices interleaving read with write + // operations is very expensive. thus, we + // here always replace the diagonal + // element, rather than first checking + // whether it is nonzero and in that case + // preserving it. this is different from + // the case of deal.II sparse matrices + // treated in the other functions. +//TODO: clear_row is not currently implemented for Trilinos + Assert (false, ExcInternalError()); +// matrix.clear_rows (constrained_rows, average_nonzero_diagonal_entry); + + // the next thing is to set right hand + // side to the wanted value. there's one + // drawback: if we write to individual + // vector elements, then we have to do + // that on all processors. however, some + // processors may not need to set + // anything because their chunk of + // matrix/rhs do not contain any boundary + // nodes. therefore, rather than using + // individual calls, we use one call for + // all elements, thereby making sure that + // all processors call this function, + // even if some only have an empty set of + // elements to set + right_hand_side.compress (); + solution.compress (); + + std::vector indices; + std::vector solution_values; + for (std::map::const_iterator + dof = boundary_values.begin(); + dof != boundary_values.end(); + ++dof) + if ((dof->first >= local_range.first) && + (dof->first < local_range.second)) + { + indices.push_back (dof->first); + solution_values.push_back (dof->second); + } + solution.set (indices, solution_values); + + // now also set appropriate values for + // the rhs + for (unsigned int i=0; i &boundary_values, + TrilinosWrappers::SparseMatrix &matrix, + TrilinosWrappers::Vector &solution, + TrilinosWrappers::Vector &right_hand_side, + const bool eliminate_columns) +{ + // simply redirect to the generic function + // used for both trilinos matrix types + TrilinosWrappers::apply_boundary_values (boundary_values, matrix, solution, + right_hand_side, eliminate_columns); +} + +#endif + + + void MatrixTools:: local_apply_boundary_values (const std::map &boundary_values,