]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Make available a whole bunch of functions in deal.II that work on Trilinos vectors...
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 21 Aug 2008 22:33:41 +0000 (22:33 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 21 Aug 2008 22:33:41 +0000 (22:33 +0000)
git-svn-id: https://svn.dealii.org/trunk@16650 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/dofs/dof_accessor.h
deal.II/deal.II/include/dofs/dof_constraints.h
deal.II/deal.II/include/fe/fe_values.h
deal.II/deal.II/include/numerics/matrices.h
deal.II/deal.II/source/dofs/dof_accessor.cc
deal.II/deal.II/source/dofs/dof_constraints.cc
deal.II/deal.II/source/fe/fe_values.cc
deal.II/deal.II/source/numerics/matrices.all_dimensions.cc

index 0f98d882851c4a5d77eded51e16d436f64d55e99..8c9cc193661b5ccedefe737cc201756eb79cb0c2 100644 (file)
@@ -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<DH::dimension,DH>
                                      * function is only callable for
                                      * active cells.
                                      *
-                                     * The input vector may be either a
-                                     * <tt>Vector<float></tt>,
+                                     * The input vector may be either
+                                     * <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,
@@ -1115,7 +1118,7 @@ class DoFCellAccessor :  public DoFObjectAccessor<DH::dimension,DH>
                                      * 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
@@ -1244,7 +1247,7 @@ class DoFCellAccessor :  public DoFObjectAccessor<DH::dimension,DH>
                                      * 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
index 153f8a1068b6ca1ab1451275024934d3cf9c34de..d04118ec4103bbe156fedb07e5f937c28efd5ccb 100644 (file)
@@ -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.
  * 
  * 
  * <h5>Condensing vectors</h5>
@@ -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.
  *
  * 
  * <h3>Avoiding explicit condensation</h3>
@@ -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".
  *
- * <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
@@ -790,11 +793,12 @@ class ConstraintMatrix : public Subscriptor
                                      * @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,
@@ -802,12 +806,14 @@ class ConstraintMatrix : public Subscriptor
 
                                     /**
                                      * 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;
@@ -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<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;
@@ -1053,24 +1061,27 @@ class ConstraintMatrix : public Subscriptor
                                      * @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;
index cc843b40ca1096d69a83253b823f22f95b350e1a..c5e0dd5d096cc6efbae510bba65f1af0ccadf665 100644 (file)
@@ -27,6 +27,7 @@
 #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>
@@ -1125,15 +1126,18 @@ class FEValuesBase : protected FEValuesData<dim>,
                                      * <tt>values</tt> object already has the
                                      * correct size. 
                                      *
-                                     * The actual data type of the input
-                                     * vector may be either a Vector&lt;T&gt;,
-                                     * BlockVector&lt;T&gt;, 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&lt;T&gt;,
+                                     * BlockVector&lt;T&gt;, 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,
@@ -1149,15 +1153,18 @@ class FEValuesBase : protected FEValuesData<dim>,
                                      * applied to multi-component
                                      * elements.
                                      *
-                                     * The actual data type of the input
-                                     * vector may be either a Vector&lt;T&gt;,
-                                     * BlockVector&lt;T&gt;, 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&lt;T&gt;,
+                                     * BlockVector&lt;T&gt;, 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,
@@ -1328,15 +1335,18 @@ class FEValuesBase : protected FEValuesData<dim>,
                                      * @p gradients object already has the
                                      * right size.
                                      *
-                                     * The actual data type of the input
-                                     * vector may be either a Vector&lt;T&gt;,
-                                     * BlockVector&lt;T&gt;, 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&lt;T&gt;,
+                                     * BlockVector&lt;T&gt;, 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<dim>,
                                      * but applied to multi-component
                                      * elements.
                                      *
-                                     * The actual data type of the input
-                                     * vector may be either a Vector&lt;T&gt;,
-                                     * BlockVector&lt;T&gt;, 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&lt;T&gt;,
+                                     * BlockVector&lt;T&gt;, 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<dim>,
                                      * get_function_hessians()
                                      * function.
                                      * 
-                                     * The actual data type of the input
-                                     * vector may be either a Vector&lt;T&gt;,
-                                     * BlockVector&lt;T&gt;, 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&lt;T&gt;,
+                                     * BlockVector&lt;T&gt;, 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<dim>,
                                      * name, but applies to
                                      * vector-valued finite elements.
                                      *
-                                     * The actual data type of the input
-                                     * vector may be either a Vector&lt;T&gt;,
-                                     * BlockVector&lt;T&gt;, 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&lt;T&gt;,
+                                     * BlockVector&lt;T&gt;, 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<dim>,
         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
     };
 
@@ -2013,6 +2056,31 @@ class FEValuesBase : protected FEValuesData<dim>,
         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
@@ -2175,6 +2243,31 @@ class FEValuesBase : protected FEValuesData<dim>,
         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
index b4c05d17c734620d827a90c0d9d70017ab3ba4ec..476afa305d2b2f388fd5ecd99d8db76357d17ebe 100644 (file)
@@ -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
+                                     * <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
@@ -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
index ab2ad3ed282fce04263ed94243d877f5b8161dfb..29b11f86cd78f7269a9a46b52f8ce049b7b80d88 100644 (file)
@@ -16,6 +16,7 @@
 #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>
@@ -646,6 +647,45 @@ DoFCellAccessor<DoFHandler<deal_II_dimension> >::set_dof_values<PETScWrappers::M
 
 
 
+#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> >::
@@ -807,6 +847,58 @@ set_dof_values_by_interpolation<PETScWrappers::MPI::BlockVector,float>
 #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> >;
@@ -957,6 +1049,44 @@ DoFCellAccessor<hp::DoFHandler<deal_II_dimension> >::set_dof_values<PETScWrapper
 (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
@@ -1098,6 +1228,55 @@ set_dof_values_by_interpolation<PETScWrappers::BlockVector,float>
 #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> >;
index c0b7e1a09afc13c5102b5419f3db3afae413c0c0..fe9d36ba8abcabf9dcbe2a6510c558c66307d741 100644 (file)
@@ -31,6 +31,8 @@
 #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>
@@ -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<SparsityPattern> (const std::vector<unsigned int> &,
                                              SparsityPattern                 &,
index 1c5de4a66561e82c0156abd9e6bd841268b5e4df..34122c3c4959ce509e321907a5f00c24c238fcf0 100644 (file)
 #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>
@@ -144,6 +142,32 @@ get_interpolated_dof_values (const PETScWrappers::BlockVector &in,
 
 #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 --------- */
 
@@ -252,6 +276,30 @@ get_interpolated_dof_values (const PETScWrappers::BlockVector &,
 
 #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
+
 
 
 
@@ -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
+
index 654b871f5841dd7231011ef69d8a273264962791..78da5349246314fb3f83c1ad502c24f494454dd6 100644 (file)
 #  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
@@ -677,6 +682,148 @@ apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
 #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,

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.