// $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
* function is only callable for
* active cells.
*
- * The input vector may be either a
- * <tt>Vector<float></tt>,
+ * The input vector may be either
+ * a <tt>Vector<float></tt>,
* Vector<double>, or a
* BlockVector<double>, 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 <class InputVector, typename number>
void get_dof_values (const InputVector &values,
* Vector<double>, or a
* BlockVector<double>, 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
* Vector<double>, or a
* BlockVector<double>, 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
* 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.
*
*
* <h5>Condensing vectors</h5>
* 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.
*
*
* <h3>Avoiding explicit condensation</h3>
* 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".
*
- * <li>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".
+ * <li>
+ * 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".
* </ul>
*
* In this case, one possibility is to distribute local entries to the final
* @p condensed be zero.
*
* The @p VectorType may be a
- * Vector<float>,
- * Vector<double>,
- * BlockVector<tt><...></tt>, a PETSc
- * vector wrapper class, or any other
- * type having the same interface.
+ * Vector<float>, Vector<double>,
+ * BlockVector<tt><...></tt>, a
+ * PETSc or Trilinos vector
+ * wrapper class, or any other
+ * type having the same
+ * interface.
*/
template <class VectorType>
void condense (const VectorType &uncondensed,
/**
* Condense the given vector
- * in-place. The @p VectorType may be a
- * Vector<float>,
+ * in-place. The @p VectorType
+ * may be a Vector<float>,
* Vector<double>,
- * BlockVector<tt><...></tt>, a PETSc
- * vector wrapper class, or any other
- * type having the same interface.
+ * BlockVector<tt><...></tt>, a
+ * PETSc or Trilinos vector
+ * wrapper class, or any other
+ * type having the same
+ * interface.
*/
template <class VectorType>
void condense (VectorType &vec) const;
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<float>,
- * Vector<double>,
- * BlockVector<tt><...></tt>, 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<float>, Vector<double>,
+ * BlockVector<tt><...></tt>, a
+ * PETSc or Trilinos vector
+ * wrapper class, or any other
+ * type having the same
+ * interface.
*/
template <class VectorType>
void set_zero (VectorType &vec) const;
* @p condense.
*
* The @p VectorType may be a
- * Vector<float>,
- * Vector<double>,
- * BlockVector<tt><...></tt>, a PETSc
- * vector wrapper class, or any other
- * type having the same interface.
+ * Vector<float>, Vector<double>,
+ * BlockVector<tt><...></tt>, a
+ * PETSc or Trilinos vector
+ * wrapper class, or any other
+ * type having the same
+ * interface.
*/
template <class VectorType>
void distribute (const VectorType &condensed,
VectorType &uncondensed) const;
/**
- * Re-distribute the elements of the
- * vector in-place. The @p VectorType
- * may be a Vector<float>,
- * Vector<double>,
- * BlockVector<tt><...></tt>, 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<float>, Vector<double>,
+ * BlockVector<tt><...></tt>, a
+ * PETSc or Trilinos vector
+ * wrapper class, or any other
+ * type having the same
+ * interface.
*/
template <class VectorType>
void distribute (VectorType &vec) const;
#include <lac/block_vector.h>
#include <lac/petsc_vector.h>
#include <lac/petsc_block_vector.h>
+#include <lac/trilinos_vector.h>
#include <grid/tria.h>
#include <grid/tria_iterator.h>
#include <dofs/dof_handler.h>
* <tt>values</tt> 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 <class InputVector, typename number>
void get_function_values (const InputVector& fe_function,
* 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 <class InputVector, typename number>
void get_function_values (const InputVector &fe_function,
* @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
* 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
* 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
* 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
void
get_interpolated_dof_values (const PETScWrappers::BlockVector &in,
Vector<PetscScalar> &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<TrilinosScalar> &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<TrilinosScalar> &out) const = 0;
#endif
};
get_interpolated_dof_values (const PETScWrappers::BlockVector &in,
Vector<PetscScalar> &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<TrilinosScalar> &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<TrilinosScalar> &out) const;
+#endif
+
private:
/**
* Copy of the iterator which
get_interpolated_dof_values (const PETScWrappers::BlockVector &in,
Vector<PetscScalar> &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<TrilinosScalar> &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<TrilinosScalar> &out) const;
+#endif
+
private:
/**
* Copy of the iterator which
}
#endif
+#ifdef DEAL_II_USE_TRILINOS
+namespace TrilinosWrappers
+{
+ class SparseMatrix;
+ class Vector;
+}
+#endif
+
/**
* This class provides functions that assemble certain standard matrices for a
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
+ * <code>true</code> 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<unsigned int,double> &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
* 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
#include <lac/block_vector.h>
#include <lac/petsc_vector.h>
#include <lac/petsc_block_vector.h>
+#include <lac/trilinos_vector.h>
#include <lac/sparse_matrix.h>
#include <dofs/dof_accessor.h>
+#ifdef DEAL_II_USE_TRILINOS
+template
+void
+DoFCellAccessor<DoFHandler<deal_II_dimension> >::get_dof_values<TrilinosWrappers::Vector,double>
+(const TrilinosWrappers::Vector &, Vector<double>&) const;
+template
+void
+DoFCellAccessor<DoFHandler<deal_II_dimension> >::get_dof_values<TrilinosWrappers::Vector,float>
+(const TrilinosWrappers::Vector &, Vector<float>&) const;
+template
+void
+DoFCellAccessor<DoFHandler<deal_II_dimension> >::set_dof_values<TrilinosWrappers::Vector,double>
+(const Vector<double> &, TrilinosWrappers::Vector&) const;
+template
+void
+DoFCellAccessor<DoFHandler<deal_II_dimension> >::set_dof_values<TrilinosWrappers::Vector,float>
+(const Vector<float>&, TrilinosWrappers::Vector&) const;
+
+// template
+// void
+// DoFCellAccessor<DoFHandler<deal_II_dimension> >::get_dof_values<TrilinosWrappers::BlockVector,double>
+// (const TrilinosWrappers::BlockVector &, Vector<double>&) const;
+// template
+// void
+// DoFCellAccessor<DoFHandler<deal_II_dimension> >::get_dof_values<TrilinosWrappers::BlockVector,float>
+// (const TrilinosWrappers::BlockVector &, Vector<float>&) const;
+// template
+// void
+// DoFCellAccessor<DoFHandler<deal_II_dimension> >::set_dof_values<TrilinosWrappers::BlockVector,double>
+// (const Vector<double> &, TrilinosWrappers::BlockVector&) const;
+// template
+// void
+// DoFCellAccessor<DoFHandler<deal_II_dimension> >::set_dof_values<TrilinosWrappers::BlockVector,float>
+// (const Vector<float>&, TrilinosWrappers::BlockVector&) const;
+
+#endif
+
+
+
template
void
DoFCellAccessor<DoFHandler<deal_II_dimension> >::
#endif
+// for Trilinos vectors
+#if DEAL_II_USE_TRILINOS
+template
+void
+DoFCellAccessor<DoFHandler<deal_II_dimension> >::
+get_interpolated_dof_values<TrilinosWrappers::Vector,double>
+(const TrilinosWrappers::Vector&, Vector<double>&) const;
+template
+void
+DoFCellAccessor<DoFHandler<deal_II_dimension> >::
+set_dof_values_by_interpolation<TrilinosWrappers::Vector,double>
+(const Vector<double>&, TrilinosWrappers::Vector&) const;
+
+template
+void
+DoFCellAccessor<DoFHandler<deal_II_dimension> >::
+get_interpolated_dof_values<TrilinosWrappers::Vector,float>
+(const TrilinosWrappers::Vector&, Vector<float>&) const;
+template
+void
+DoFCellAccessor<DoFHandler<deal_II_dimension> >::
+set_dof_values_by_interpolation<TrilinosWrappers::Vector,float>
+(const Vector<float>&, TrilinosWrappers::Vector&) const;
+
+
+// template
+// void
+// DoFCellAccessor<DoFHandler<deal_II_dimension> >::
+// get_interpolated_dof_values<TrilinosWrappers::BlockVector,double>
+// (const TrilinosWrappers::BlockVector&, Vector<double>&) const;
+// template
+// void
+// DoFCellAccessor<DoFHandler<deal_II_dimension> >::
+// set_dof_values_by_interpolation<TrilinosWrappers::BlockVector,double>
+// (const Vector<double>&, TrilinosWrappers::BlockVector&) const;
+// template
+
+// template
+// void
+// DoFCellAccessor<DoFHandler<deal_II_dimension> >::
+// get_interpolated_dof_values<TrilinosWrappers::BlockVector,float>
+// (const TrilinosWrappers::BlockVector&, Vector<float>&) const;
+// template
+// void
+// DoFCellAccessor<DoFHandler<deal_II_dimension> >::
+// set_dof_values_by_interpolation<TrilinosWrappers::BlockVector,float>
+// (const Vector<float>&, TrilinosWrappers::BlockVector&) const;
+// template
+
+#endif
+
+
#if deal_II_dimension == 1
template class DoFAccessor<1, DoFHandler<1> >;
template class DoFObjectAccessor<1, DoFHandler<1> >;
(const Vector<float>&, PETScWrappers::BlockVector&) const;
#endif
+// for Trilinos vectors
+#ifdef DEAL_II_USE_TRILINOS
+template
+void
+DoFCellAccessor<hp::DoFHandler<deal_II_dimension> >::get_dof_values<TrilinosWrappers::Vector,double>
+(const TrilinosWrappers::Vector &, Vector<double>&) const;
+template
+void
+DoFCellAccessor<hp::DoFHandler<deal_II_dimension> >::get_dof_values<TrilinosWrappers::Vector,float>
+(const TrilinosWrappers::Vector &, Vector<float>&) const;
+template
+void
+DoFCellAccessor<hp::DoFHandler<deal_II_dimension> >::set_dof_values<TrilinosWrappers::Vector,double>
+(const Vector<double> &, TrilinosWrappers::Vector&) const;
+template
+void
+DoFCellAccessor<hp::DoFHandler<deal_II_dimension> >::set_dof_values<TrilinosWrappers::Vector,float>
+(const Vector<float>&, TrilinosWrappers::Vector&) const;
+
+// template
+// void
+// DoFCellAccessor<hp::DoFHandler<deal_II_dimension> >::get_dof_values<TrilinosWrappers::BlockVector,double>
+// (const TrilinosWrappers::BlockVector &, Vector<double>&) const;
+// template
+// void
+// DoFCellAccessor<hp::DoFHandler<deal_II_dimension> >::get_dof_values<TrilinosWrappers::BlockVector,float>
+// (const TrilinosWrappers::BlockVector &, Vector<float>&) const;
+// template
+// void
+// DoFCellAccessor<hp::DoFHandler<deal_II_dimension> >::set_dof_values<TrilinosWrappers::BlockVector,double>
+// (const Vector<double> &, TrilinosWrappers::BlockVector&) const;
+// template
+// void
+// DoFCellAccessor<hp::DoFHandler<deal_II_dimension> >::set_dof_values<TrilinosWrappers::BlockVector,float>
+// (const Vector<float>&, TrilinosWrappers::BlockVector&) const;
+#endif
+
+
template
void
#endif
+// for Trilinos vectors
+#ifdef DEAL_II_USE_TRILINOS
+template
+void
+DoFCellAccessor<hp::DoFHandler<deal_II_dimension> >::
+get_interpolated_dof_values<TrilinosWrappers::Vector,double>
+(const TrilinosWrappers::Vector&, Vector<double>&) const;
+template
+void
+DoFCellAccessor<hp::DoFHandler<deal_II_dimension> >::
+set_dof_values_by_interpolation<TrilinosWrappers::Vector,double>
+(const Vector<double>&, TrilinosWrappers::Vector&) const;
+
+template
+void
+DoFCellAccessor<hp::DoFHandler<deal_II_dimension> >::
+get_interpolated_dof_values<TrilinosWrappers::Vector,float>
+(const TrilinosWrappers::Vector&, Vector<float>&) const;
+template
+void
+DoFCellAccessor<hp::DoFHandler<deal_II_dimension> >::
+set_dof_values_by_interpolation<TrilinosWrappers::Vector,float>
+(const Vector<float>&, TrilinosWrappers::Vector&) const;
+
+
+// template
+// void
+// DoFCellAccessor<hp::DoFHandler<deal_II_dimension> >::
+// get_interpolated_dof_values<TrilinosWrappers::BlockVector,double>
+// (const TrilinosWrappers::BlockVector&, Vector<double>&) const;
+// template
+// void
+// DoFCellAccessor<hp::DoFHandler<deal_II_dimension> >::
+// set_dof_values_by_interpolation<TrilinosWrappers::BlockVector,double>
+// (const Vector<double>&, TrilinosWrappers::BlockVector&) const;
+
+// template
+// void
+// DoFCellAccessor<hp::DoFHandler<deal_II_dimension> >::
+// get_interpolated_dof_values<TrilinosWrappers::BlockVector,float>
+// (const TrilinosWrappers::BlockVector&, Vector<float>&) const;
+// template
+// void
+// DoFCellAccessor<hp::DoFHandler<deal_II_dimension> >::
+// set_dof_values_by_interpolation<TrilinosWrappers::BlockVector,float>
+// (const Vector<float>&, TrilinosWrappers::BlockVector&) const;
+#endif
+
+
#if deal_II_dimension == 1
template class DoFAccessor<1, hp::DoFHandler<1> >;
template class DoFObjectAccessor<1, hp::DoFHandler<1> >;
#include <lac/petsc_parallel_block_vector.h>
#include <lac/petsc_parallel_sparse_matrix.h>
#include <lac/petsc_parallel_block_sparse_matrix.h>
+#include <lac/trilinos_vector.h>
+#include <lac/trilinos_sparse_matrix.h>
#include <algorithm>
#include <numeric>
VECTOR_FUNCTIONS(PETScWrappers::MPI::BlockVector);
#endif
+#ifdef DEAL_II_USE_TRILINOS
+VECTOR_FUNCTIONS(TrilinosWrappers::Vector);
+//VECTOR_FUNCTIONS(TrilinosWrappers::BlockVector);
+#endif
+
template
void
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<SparsityPattern> (const std::vector<unsigned int> &,
SparsityPattern &,
#include <lac/block_vector.h>
#include <lac/petsc_vector.h>
#include <lac/petsc_block_vector.h>
+#include <lac/trilinos_vector.h>
#include <grid/tria_iterator.h>
#include <grid/tria_accessor.h>
#include <grid/tria_boundary.h>
#include <dofs/dof_accessor.h>
-#include <lac/vector.h>
-#include <lac/block_vector.h>
-#include <lac/petsc_vector.h>
#include <fe/mapping_q1.h>
#include <fe/fe_values.h>
#include <fe/fe.h>
#endif
+#ifdef DEAL_II_USE_TRILINOS
+
+template <int dim>
+template <typename CI>
+void
+FEValuesBase<dim>::CellIterator<CI>::
+get_interpolated_dof_values (const TrilinosWrappers::Vector &in,
+ Vector<TrilinosScalar> &out) const
+{
+ cell->get_interpolated_dof_values (in, out);
+}
+
+
+
+// template <int dim>
+// template <typename CI>
+// void
+// FEValuesBase<dim>::CellIterator<CI>::
+// get_interpolated_dof_values (const TrilinosWrappers::BlockVector &in,
+// Vector<TrilinosScalar> &out) const
+// {
+// cell->get_interpolated_dof_values (in, out);
+// }
+
+#endif
+
/* ---------------- FEValuesBase<dim>::TriaCellIterator --------- */
#endif
+#ifdef DEAL_II_USE_TRILINOS
+
+template <int dim>
+void
+FEValuesBase<dim>::TriaCellIterator::
+get_interpolated_dof_values (const TrilinosWrappers::Vector &,
+ Vector<TrilinosScalar> &) const
+{
+ Assert (false, ExcMessage (message_string));
+}
+
+
+
+// template <int dim>
+// void
+// FEValuesBase<dim>::TriaCellIterator::
+// get_interpolated_dof_values (const TrilinosWrappers::BlockVector &,
+// Vector<TrilinosScalar> &) const
+// {
+// Assert (false, ExcMessage (message_string));
+// }
+
+#endif
+
#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
+
# include <lac/petsc_vector.h>
#endif
+#ifdef DEAL_II_USE_TRILINOS
+# include <lac/trilinos_sparse_matrix.h>
+# include <lac/trilinos_vector.h>
+#endif
+
#include <algorithm>
DEAL_II_NAMESPACE_OPEN
#endif
+
+#ifdef DEAL_II_USE_TRILINOS
+
+namespace TrilinosWrappers
+{
+ template <typename TrilinosMatrix, typename TrilinosVector>
+ void
+ apply_boundary_values (const std::map<unsigned int,double> &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<unsigned int, unsigned int> 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<local_range.second; ++i)
+ if (matrix.diag_element(i) != 0)
+ {
+ average_nonzero_diagonal_entry = std::fabs(matrix.diag_element(i));
+ break;
+ }
+
+ // figure out which rows of the matrix we
+ // have to eliminate on this processor
+ std::vector<unsigned int> constrained_rows;
+ for (std::map<unsigned int,double>::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<unsigned int> indices;
+ std::vector<TrilinosScalar> solution_values;
+ for (std::map<unsigned int,double>::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<solution_values.size(); ++i)
+ solution_values[i] *= average_nonzero_diagonal_entry;
+
+ right_hand_side.set (indices, solution_values);
+
+ // clean up
+ matrix.compress ();
+ solution.compress ();
+ right_hand_side.compress ();
+ }
+}
+
+
+
+void
+MatrixTools::
+apply_boundary_values (const std::map<unsigned int,double> &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<unsigned int,double> &boundary_values,