]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implement function get_constraint_entries that allows queries to the constraints...
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Wed, 11 Aug 2010 13:22:01 +0000 (13:22 +0000)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Wed, 11 Aug 2010 13:22:01 +0000 (13:22 +0000)
git-svn-id: https://svn.dealii.org/trunk@21641 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/block_indices.h
deal.II/lac/include/lac/block_matrix_base.h
deal.II/lac/include/lac/constraint_matrix.h
deal.II/lac/include/lac/constraint_matrix.templates.h
deal.II/lac/source/constraint_matrix.cc

index 70da23a366c2d031e0375bcadcc3e238312045d0..85a2c1c0b33b401b332651b4e03c083d784ca700 100644 (file)
@@ -164,6 +164,89 @@ class BlockIndices : public Subscriptor
 
 
 
+template <typename MatrixType> class BlockMatrixBase;
+template <typename SparsityType> class BlockSparsityPatternBase;
+template <typename number>     class BlockSparseMatrixEZ;
+
+/**
+ * A class that can be used to determine whether a given type is a block
+ * matrix type or not. For example,
+ * @code
+ *   IsBlockMatrix<SparseMatrix<double> >::value
+ * @endcode
+ * has the value false, whereas
+ * @code
+ *   IsBlockMatrix<BlockSparseMatrix<double> >::value
+ * @endcode
+ * is true. This is sometimes useful in template contexts where we may
+ * want to do things differently depending on whether a template type
+ * denotes a regular or a block matrix type.
+ *
+ * @see @ref GlossBlockLA "Block (linear algebra)"
+ * @author Wolfgang Bangerth, 2009
+ */
+template <typename MatrixType>
+struct IsBlockMatrix
+{
+  private:
+    struct yes_type { char c[1]; };
+    struct no_type  { char c[2]; };
+
+                                    /**
+                                     * Overload returning true if the class
+                                     * is derived from BlockMatrixBase,
+                                     * which is what block matrices do
+                                     * (with the exception of
+                                     * BlockSparseMatrixEZ).
+                                     */
+    template <typename T>
+    static yes_type check_for_block_matrix (const BlockMatrixBase<T> *);
+
+                                    /**
+                                     * Overload returning true if the class
+                                     * is derived from
+                                     * BlockSparsityPatternBase, which is
+                                     * what block sparsity patterns do.
+                                     */
+    template <typename T>
+    static yes_type check_for_block_matrix (const BlockSparsityPatternBase<T> *);
+
+                                    /**
+                                     * Overload for BlockSparseMatrixEZ,
+                                     * which is the only block matrix not
+                                     * derived from BlockMatrixBase at the
+                                     * time of writing this class.
+                                     */
+    template <typename T>
+    static yes_type check_for_block_matrix (const BlockSparseMatrixEZ<T> *);
+
+                                    /**
+                                     * Catch all for all other potential
+                                     * matrix types that are not block
+                                     * matrices.
+                                     */
+    static no_type check_for_block_matrix (...);
+
+  public:
+                                    /**
+                                     * A statically computable value that
+                                     * indicates whether the template
+                                     * argument to this class is a block
+                                     * matrix (in fact whether the type is
+                                     * derived from BlockMatrixBase<T>).
+                                     */
+    static const bool value = (sizeof(check_for_block_matrix
+                                     ((MatrixType*)0))
+                              ==
+                              sizeof(yes_type));
+};
+
+
+// instantiation of the static member
+template <typename MatrixType>
+const bool IsBlockMatrix<MatrixType>::value;
+
+
 /* ---------------------- template and inline functions ------------------- */
 
 inline
index 8c63c06eea6da6b3cde66c9a854e84c50f181ee4..963f0eaa15a2816393c5655eac680dc8c86b9fc1 100644 (file)
@@ -28,87 +28,6 @@ DEAL_II_NAMESPACE_OPEN
 
 
 template <typename> class MatrixIterator;
-template <typename MatrixType> class BlockMatrixBase;
-template <typename SparsityType> class BlockSparsityPatternBase;
-template <typename number>     class BlockSparseMatrixEZ;
-
-/**
- * A class that can be used to determine whether a given type is a block
- * matrix type or not. For example,
- * @code
- *   IsBlockMatrix<SparseMatrix<double> >::value
- * @endcode
- * has the value false, whereas
- * @code
- *   IsBlockMatrix<BlockSparseMatrix<double> >::value
- * @endcode
- * is true. This is sometimes useful in template contexts where we may
- * want to do things differently depending on whether a template type
- * denotes a regular or a block matrix type.
- *
- * @see @ref GlossBlockLA "Block (linear algebra)"
- * @author Wolfgang Bangerth, 2009
- */
-template <typename MatrixType>
-struct IsBlockMatrix
-{
-  private:
-    struct yes_type { char c[1]; };
-    struct no_type  { char c[2]; };
-
-                                    /**
-                                     * Overload returning true if the class
-                                     * is derived from BlockMatrixBase,
-                                     * which is what block matrices do
-                                     * (with the exception of
-                                     * BlockSparseMatrixEZ).
-                                     */
-    template <typename T>
-    static yes_type check_for_block_matrix (const BlockMatrixBase<T> *);
-
-                                    /**
-                                     * Overload returning true if the class
-                                     * is derived from
-                                     * BlockSparsityPatternBase, which is
-                                     * what block sparsity patterns do.
-                                     */
-    template <typename T>
-    static yes_type check_for_block_matrix (const BlockSparsityPatternBase<T> *);
-
-                                    /**
-                                     * Overload for BlockSparseMatrixEZ,
-                                     * which is the only block matrix not
-                                     * derived from BlockMatrixBase at the
-                                     * time of writing this class.
-                                     */
-    template <typename T>
-    static yes_type check_for_block_matrix (const BlockSparseMatrixEZ<T> *);
-
-                                    /**
-                                     * Catch all for all other potential
-                                     * matrix types that are not block
-                                     * matrices.
-                                     */
-    static no_type check_for_block_matrix (...);
-
-  public:
-                                    /**
-                                     * A statically computable value that
-                                     * indicates whether the template
-                                     * argument to this class is a block
-                                     * matrix (in fact whether the type is
-                                     * derived from BlockMatrixBase<T>).
-                                     */
-    static const bool value = (sizeof(check_for_block_matrix
-                                     ((MatrixType*)0))
-                              ==
-                              sizeof(yes_type));
-};
-
-
-// instantiation of the static member
-template <typename MatrixType>
-const bool IsBlockMatrix<MatrixType>::value;
 
 
 
index 4796d5219b1b51e90d9ad608f7e08ae89c4d7029..b3296ce25c02a018dcfa3a87831b8fd5bddf1776 100644 (file)
@@ -21,7 +21,6 @@
 #include <base/template_constraints.h>
 
 #include <lac/vector.h>
-#include <lac/trilinos_vector.h>
 
 #include <vector>
 #include <map>
@@ -47,6 +46,7 @@ class BlockCompressedSimpleSparsityPattern;
 template <typename number> class SparseMatrix;
 template <typename number> class BlockSparseMatrix;
 class BlockIndices;
+template <typename MatrixType> struct IsBlockMatrix;
 
 namespace internals
 {
@@ -371,10 +371,27 @@ class ConstraintMatrix : public Subscriptor
                                      * This function copies the content of @p
                                      * constraints_in with DoFs that are
                                      * element of the IndexSet @p
-                                     * filter. Constrained dofs are
-                                     * transformed to local index space of
-                                     * the filter, and elements not present
-                                     * in the IndexSet are ignored.
+                                     * filter. Elements that are not present
+                                     * in the IndexSet are ignored. All DoFs
+                                     * will be transformed to local index
+                                     * space of the filter, both the
+                                     * constrained DoFs and the other DoFs
+                                     * these entries are constrained to. The
+                                     * local index space of the filter is a
+                                     * contiguous numbering of all (global)
+                                     * DoFs that are elements in the
+                                     * filter.
+                                     *
+                                     * If, for example, the filter represents
+                                     * the range <tt>[10,20)</tt>, and the
+                                     * constraint matrix @p constraints_in
+                                     * includes the global indices
+                                     * <tt>{7,13,14}</tt>, the indices
+                                     * <tt>{3,4}</tt> are added to the
+                                     * calling constraint matrix (since 13
+                                     * and 14 are elements in the filter and
+                                     * element 13 is the fourth element in
+                                     * the index, and 14 is the fifth).
                                      *
                                      * This function provides an easy way to
                                      * create a ConstraintMatrix for certain
@@ -744,6 +761,15 @@ class ConstraintMatrix : public Subscriptor
                                      */
     bool has_inhomogeneities () const;
 
+                                     /**
+                                      * Returns a pointer to the the vector of
+                                      * entries if a line is constrained, and a
+                                      * zero pointer in case the dof is not
+                                      * constrained.
+                                      */
+    const std::vector<std::pair<unsigned int,double> >*
+    get_constraint_entries (unsigned int line) const;
+
                                     /**
                                      * Print the constraint lines. Mainly
                                      * for debugging purposes.
@@ -1918,15 +1944,6 @@ class ConstraintMatrix : public Subscriptor
                          const Vector<double>                 &local_vector,
                          const std::vector<unsigned int>      &local_dof_indices,
                          const FullMatrix<double>             &local_matrix) const;
-
-#ifdef DEAL_II_USE_TRILINOS
-//TODO: Make use of the following member thread safe
-                                     /**
-                                      * This vector is used to import data
-                                      * within the distribute function.
-                                      */
-    mutable boost::scoped_ptr<TrilinosWrappers::MPI::Vector> vec_distribute;
-#endif
 };
 
 
@@ -1951,9 +1968,6 @@ ConstraintMatrix::ConstraintMatrix (const ConstraintMatrix &constraint_matrix)
                 lines_cache (constraint_matrix.lines_cache),
                 local_lines (constraint_matrix.local_lines),
                 sorted (constraint_matrix.sorted)
-#ifdef DEAL_II_USE_TRILINOS
-               ,vec_distribute ()
-#endif
 {}
 
 
@@ -2069,6 +2083,41 @@ ConstraintMatrix::is_inhomogeneously_constrained (const unsigned int index) cons
 
 
 
+inline
+const std::vector<std::pair<unsigned int,double> >*
+ConstraintMatrix::get_constraint_entries (unsigned int line) const
+{
+  if (is_constrained(line))
+    return &lines[lines_cache[calculate_line_index(line)]].entries;
+  else
+    return 0;
+}
+
+
+
+inline unsigned int
+ConstraintMatrix::calculate_line_index (const unsigned int line) const
+{
+                                  //IndexSet is unused (serial case)
+  if (!local_lines.size())
+    return line;
+
+  Assert(local_lines.is_element(line),
+        ExcRowNotStoredHere(line));
+
+  return local_lines.index_within_set(line);
+}
+
+
+
+inline bool
+ConstraintMatrix::can_store_line(unsigned int line_index) const
+{
+  return !local_lines.size() || local_lines.is_element(line_index);
+}
+
+
+
 template <class InVector, class OutVector>
 inline
 void
@@ -2105,14 +2154,8 @@ void ConstraintMatrix::
          const ConstraintLine& position =
            lines[lines_cache[calculate_line_index(*local_indices_begin)]];
          for (unsigned int j=0; j<position.entries.size(); ++j)
-           {
-             Assert (!(!local_lines.size()
-                       || local_lines.is_element(position.entries[j].first))
-                     || is_constrained(position.entries[j].first) == false,
-                     ExcMessage ("Tried to distribute to a fixed dof."));
-             global_vector(position.entries[j].first)
-               += *local_vector_begin * position.entries[j].second;
-           }
+           global_vector(position.entries[j].first)
+             += *local_vector_begin * position.entries[j].second;
        }
     }
 }
@@ -2139,12 +2182,8 @@ void ConstraintMatrix::get_dof_values (const VectorType  &global_vector,
            lines[lines_cache[calculate_line_index(*local_indices_begin)]];
          typename VectorType::value_type value = position.inhomogeneity;
          for (unsigned int j=0; j<position.entries.size(); ++j)
-           {
-             Assert (is_constrained(position.entries[j].first) == false,
-                     ExcMessage ("Tried to distribute to a fixed dof."));
-             value += (global_vector(position.entries[j].first) *
-                       position.entries[j].second);
-           }
+           value += (global_vector(position.entries[j].first) *
+                     position.entries[j].second);
          *local_vector_begin = value;
        }
     }
@@ -2152,27 +2191,65 @@ void ConstraintMatrix::get_dof_values (const VectorType  &global_vector,
 
 
 
-inline unsigned int
-ConstraintMatrix::calculate_line_index (const unsigned int line) const
+template <typename MatrixType>
+inline
+void
+ConstraintMatrix::
+distribute_local_to_global (const FullMatrix<double>        &local_matrix,
+                            const std::vector<unsigned int> &local_dof_indices,
+                            MatrixType                      &global_matrix) const
 {
-                                  //IndexSet is unused (serial case)
-  if (!local_lines.size())
-    return line;
+                                   // create a dummy and hand on to the
+                                   // function actually implementing this
+                                   // feature in the cm.templates.h file.
+  Vector<double> dummy(0);
+  distribute_local_to_global (local_matrix, dummy, local_dof_indices,
+                              global_matrix, dummy,
+                              internal::bool2type<IsBlockMatrix<MatrixType>::value>());
+}
 
-  Assert(local_lines.is_element(line),
-        ExcRowNotStoredHere(line));
 
-  return local_lines.index_within_set(line);
-}
 
-inline bool
-ConstraintMatrix::can_store_line(unsigned int line_index) const
+template <typename MatrixType, typename VectorType>
+inline
+void
+ConstraintMatrix::
+distribute_local_to_global (const FullMatrix<double>        &local_matrix,
+                            const Vector<double>            &local_vector,
+                            const std::vector<unsigned int> &local_dof_indices,
+                            MatrixType                      &global_matrix,
+                            VectorType                      &global_vector) const
 {
-  return !local_lines.size() || local_lines.is_element(line_index);
+                                   // enter the internal function with the
+                                   // respective block information set, the
+                                   // actual implementation follows in the
+                                   // cm.templates.h file.
+  distribute_local_to_global (local_matrix, local_vector, local_dof_indices,
+                              global_matrix, global_vector,
+                              internal::bool2type<IsBlockMatrix<MatrixType>::value>());
 }
 
 
 
+template <typename SparsityType>
+inline
+void
+ConstraintMatrix::
+add_entries_local_to_global (const std::vector<unsigned int> &local_dof_indices,
+                             SparsityType                    &sparsity_pattern,
+                             const bool                       keep_constrained_entries,
+                             const Table<2,bool>             &dof_mask) const
+{
+                                   // enter the internal function with the
+                                   // respective block information set, the
+                                   // actual implementation follows in the
+                                   // cm.templates.h file.
+  add_entries_local_to_global (local_dof_indices, sparsity_pattern,
+                               keep_constrained_entries, dof_mask,
+                               internal::bool2type<IsBlockMatrix<SparsityType>::value>());
+}
+
+
 DEAL_II_NAMESPACE_CLOSE
 
 #endif
index c1cb2362a643446ba65c359c8af4ca45e7eec15d..88ffafa09bb9efac16063be879a9ba215babdfbd 100644 (file)
@@ -803,63 +803,6 @@ distribute_local_to_global (const Vector<double>            &local_vector,
 
 
 
-template <typename MatrixType>
-void
-ConstraintMatrix::
-distribute_local_to_global (const FullMatrix<double>        &local_matrix,
-                            const std::vector<unsigned int> &local_dof_indices,
-                            MatrixType                      &global_matrix) const
-{
-                                  // create a dummy and hand on to the
-                                  // function actually implementing this
-                                  // feature further down.
-  Vector<double> dummy(0);
-  distribute_local_to_global (local_matrix, dummy, local_dof_indices,
-                             global_matrix, dummy,
-                             internal::bool2type<IsBlockMatrix<MatrixType>::value>());
-}
-
-
-
-template <typename MatrixType, typename VectorType>
-void
-ConstraintMatrix::
-distribute_local_to_global (const FullMatrix<double>        &local_matrix,
-                           const Vector<double>            &local_vector,
-                            const std::vector<unsigned int> &local_dof_indices,
-                            MatrixType                      &global_matrix,
-                           VectorType                      &global_vector) const
-{
-                                  // enter the internal function with the
-                                  // respective block information set, the
-                                  // actual implementation follows further
-                                  // down.
-  distribute_local_to_global (local_matrix, local_vector, local_dof_indices,
-                             global_matrix, global_vector,
-                             internal::bool2type<IsBlockMatrix<MatrixType>::value>());
-}
-
-
-
-template <typename SparsityType>
-void
-ConstraintMatrix::
-add_entries_local_to_global (const std::vector<unsigned int> &local_dof_indices,
-                            SparsityType                    &sparsity_pattern,
-                            const bool                       keep_constrained_entries,
-                            const Table<2,bool>             &dof_mask) const
-{
-                                  // enter the internal function with the
-                                  // respective block information set, the
-                                  // actual implementation follows further
-                                  // down.
-  add_entries_local_to_global (local_dof_indices, sparsity_pattern,
-                              keep_constrained_entries, dof_mask,
-                              internal::bool2type<IsBlockMatrix<SparsityType>::value>());
-}
-
-
-
 template<class VectorType>
 void
 ConstraintMatrix::distribute (const VectorType &condensed,
@@ -1005,7 +948,7 @@ namespace internals
     local_row = in.local_row;
                                   // the constraints pointer should not
                                   // contain any data here.
-    Assert (constraint_position == numbers::invalid_unsigned_int, 
+    Assert (constraint_position == numbers::invalid_unsigned_int,
            ExcInternalError());
 
     if (in.constraint_position != numbers::invalid_unsigned_int)
@@ -1123,13 +1066,13 @@ namespace internals
                                // "Distributing" structs using access via
                                // the DataCache. Provides some
                                // specialized sort and insert functions.
-                               // 
+                               //
                                // in case there are no constraints, this is
                                // basically a list of pairs <uint,unit> with
                                // the first index being the global index and
                                // the second index the local index. The list
                                // is sorted with respect to the global index.
-                               // 
+                               //
                                // in case there are constraints, a global dof
                                // might get a contribution also because it
                                // gets data from a constrained dof. This
@@ -1246,7 +1189,7 @@ namespace internals
                   total_row_indices.back().local_row);
        ++n_inhomogeneous_rows; };
 
-                               // 
+                               //
     unsigned int inhomogeneity (unsigned int i) const
       { return total_row_indices[n_active_rows+i].local_row; };
 
@@ -1457,7 +1400,7 @@ namespace internals
                                   // the innermost loop. goes through the
                                   // origin of each global entry and finds
                                   // out which data we need to collect.
-  inline
+  static inline
   double resolve_matrix_entry (const GlobalRowsFromLocal&global_rows,
                                const GlobalRowsFromLocal&global_cols,
                               const unsigned int        i,
@@ -1474,7 +1417,7 @@ namespace internals
                                   // set the value to zero.
     if (loc_row != numbers::invalid_unsigned_int)
       {
-       col_val = ((loc_col != numbers::invalid_unsigned_int) ? 
+       col_val = ((loc_col != numbers::invalid_unsigned_int) ?
                   local_matrix(loc_row, loc_col) : 0);
 
                                   // account for indirect contributions by
@@ -1515,7 +1458,7 @@ namespace internals
                                // resulting values in val_ptr, and the
                                // corresponding column indices in col_ptr.
   template <typename number>
-  inline
+  static inline
   void
   resolve_matrix_row (const GlobalRowsFromLocal&global_rows,
                       const GlobalRowsFromLocal&global_cols,
@@ -1582,7 +1525,7 @@ namespace internals
   namespace dealiiSparseMatrix
   {
     template <typename number>
-    inline
+    static inline
     void add_value (const double value,
                    const unsigned int row,
                    const unsigned int column,
@@ -1616,7 +1559,7 @@ namespace internals
                                // operations just in place, i.e., in the
                                // respective matrix row
   template <typename number>
-  inline
+  static inline
   void
   resolve_matrix_row (const GlobalRowsFromLocal&global_rows,
                      const unsigned int        i,
@@ -1668,9 +1611,9 @@ namespace internals
              {
                const unsigned int loc_col = global_rows.local_row(j);
                const double col_val = matrix_ptr[loc_col];
-               dealiiSparseMatrix::add_value (col_val, row, 
-                                              global_rows.global_row(j), 
-                                              col_ptr, false, counter, 
+               dealiiSparseMatrix::add_value (col_val, row,
+                                              global_rows.global_row(j),
+                                              col_ptr, false, counter,
                                               val_ptr);
              }
          }
@@ -1680,7 +1623,7 @@ namespace internals
              {
                double col_val = resolve_matrix_entry (global_rows, global_rows, i, j,
                                                       loc_row, local_matrix);
-               dealiiSparseMatrix::add_value (col_val, row, 
+               dealiiSparseMatrix::add_value (col_val, row,
                                               global_rows.global_row(j), col_ptr,
                                               false, counter, val_ptr);
              }
@@ -1698,7 +1641,7 @@ namespace internals
              {
                const unsigned int loc_col = global_rows.local_row(j);
                const double col_val = matrix_ptr[loc_col];
-               dealiiSparseMatrix::add_value(col_val, row, 
+               dealiiSparseMatrix::add_value(col_val, row,
                                              global_rows.global_row(j), col_ptr,
                                              false, counter, val_ptr);
              }
@@ -1707,7 +1650,7 @@ namespace internals
              {
                const unsigned int loc_col = global_rows.local_row(j);
                const double col_val = matrix_ptr[loc_col];
-               dealiiSparseMatrix::add_value(col_val, row, 
+               dealiiSparseMatrix::add_value(col_val, row,
                                              global_rows.global_row(j), col_ptr,
                                              false, counter, val_ptr);
              }
@@ -1718,7 +1661,7 @@ namespace internals
              {
                double col_val = resolve_matrix_entry (global_rows, global_rows, i, j,
                                                       loc_row, local_matrix);
-               dealiiSparseMatrix::add_value (col_val, row, 
+               dealiiSparseMatrix::add_value (col_val, row,
                                               global_rows.global_row(j), col_ptr,
                                               false, counter, val_ptr);
              }
@@ -1728,7 +1671,7 @@ namespace internals
              {
                double col_val = resolve_matrix_entry (global_rows, global_rows, i, j,
                                                       loc_row, local_matrix);
-               dealiiSparseMatrix::add_value (col_val, row, 
+               dealiiSparseMatrix::add_value (col_val, row,
                                               global_rows.global_row(j), col_ptr,
                                               false, counter, val_ptr);
              }
@@ -1746,9 +1689,9 @@ namespace internals
          {
            const unsigned int loc_col = global_rows.local_row(j);
            const double col_val = matrix_ptr[loc_col];
-           dealiiSparseMatrix::add_value(col_val, row, 
+           dealiiSparseMatrix::add_value(col_val, row,
                                          global_rows.global_row(j), col_ptr,
-                                         row==global_rows.global_row(j), 
+                                         row==global_rows.global_row(j),
                                          counter, val_ptr);
          }
       }
@@ -1758,9 +1701,9 @@ namespace internals
          {
            double col_val = resolve_matrix_entry (global_rows, global_rows, i, j,
                                                   loc_row, local_matrix);
-           dealiiSparseMatrix::add_value (col_val, row, 
+           dealiiSparseMatrix::add_value (col_val, row,
                                           global_rows.global_row(j), col_ptr,
-                                          row==global_rows.global_row(j), 
+                                          row==global_rows.global_row(j),
                                           counter, val_ptr);
          }
       }
@@ -1772,7 +1715,7 @@ namespace internals
                                // will be added to the given global row
                                // global_rows[i] as before, now for sparsity
                                // pattern
-  inline
+  static inline
   void
   resolve_matrix_row (const GlobalRowsFromLocal     &global_rows,
                      const unsigned int             i,
@@ -2306,8 +2249,9 @@ ConstraintMatrix::distribute_local_to_global (
        {
          unsigned int * col_ptr = &cols[0];
          number * val_ptr = &vals[0];
-         resolve_matrix_row (global_rows, global_rows, i, 0, n_actual_dofs,
-                             local_matrix, col_ptr, val_ptr);
+         internals::resolve_matrix_row (global_rows, global_rows, i, 0,
+                                        n_actual_dofs,
+                                        local_matrix, col_ptr, val_ptr);
          const unsigned int n_values = col_ptr - &cols[0];
          Assert (n_values == (unsigned int)(val_ptr - &vals[0]),
                  ExcInternalError());
@@ -2315,8 +2259,8 @@ ConstraintMatrix::distribute_local_to_global (
            global_matrix.add(row, n_values, &cols[0], &vals[0], false, true);
        }
       else
-       resolve_matrix_row (global_rows, i, 0, n_actual_dofs,
-                           local_matrix, sparse_matrix);
+       internals::resolve_matrix_row (global_rows, i, 0, n_actual_dofs,
+                                      local_matrix, sparse_matrix);
 
                                   // now to the vectors. besides doing the
                                   // same job as we did above (i.e.,
@@ -2381,8 +2325,9 @@ ConstraintMatrix::distribute_local_to_global (
                                   // written into the matrix row.
          unsigned int * col_ptr = &cols[0];
          number * val_ptr = &vals[0];
-         resolve_matrix_row (global_rows, global_cols, i, 0, n_actual_col_dofs,
-                             local_matrix, col_ptr, val_ptr);
+         internals::resolve_matrix_row (global_rows, global_cols, i, 0,
+                                        n_actual_col_dofs,
+                                        local_matrix, col_ptr, val_ptr);
          const unsigned int n_values = col_ptr - &cols[0];
          Assert (n_values == (unsigned int)(val_ptr - &vals[0]),
                  ExcInternalError());
@@ -2476,8 +2421,9 @@ distribute_local_to_global (const FullMatrix<double>        &local_matrix,
                {
                  unsigned int * col_ptr = &cols[0];
                  number * val_ptr = &vals[0];
-                 resolve_matrix_row (global_rows, global_rows, i, start_block,
-                                     end_block, local_matrix, col_ptr, val_ptr);
+                 internals::resolve_matrix_row (global_rows, global_rows, i,
+                                                start_block, end_block,
+                                                local_matrix, col_ptr, val_ptr);
                  const unsigned int n_values = col_ptr - &cols[0];
                  Assert (n_values == (unsigned int)(val_ptr - &vals[0]),
                          ExcInternalError());
@@ -2492,8 +2438,8 @@ distribute_local_to_global (const FullMatrix<double>        &local_matrix,
                    = dynamic_cast<SparseMatrix<number> *>(&global_matrix.block(block,
                                                                                block_col));
                  Assert (sparse_matrix != 0, ExcInternalError());
-                 resolve_matrix_row (global_rows, i, start_block,
-                                     end_block, local_matrix, sparse_matrix);
+                 internals::resolve_matrix_row (global_rows, i, start_block,
+                                                end_block, local_matrix, sparse_matrix);
                }
            }
 
@@ -2580,8 +2526,8 @@ add_entries_local_to_global (const std::vector<unsigned int> &local_dof_indices,
     {
       std::vector<unsigned int>::iterator col_ptr = cols.begin();
       const unsigned int row = global_rows.global_row(i);
-      resolve_matrix_row (global_rows, i, 0, n_actual_dofs,
-                         dof_mask, col_ptr);
+      internals::resolve_matrix_row (global_rows, i, 0, n_actual_dofs,
+                                    dof_mask, col_ptr);
 
                                   // finally, write all the information
                                   // that accumulated under the given
@@ -2693,8 +2639,8 @@ add_entries_local_to_global (const std::vector<unsigned int> &local_dof_indices,
              const unsigned int begin_block = block_starts[block_col],
                end_block = block_starts[block_col+1];
              std::vector<unsigned int>::iterator col_ptr = cols.begin();
-             resolve_matrix_row (global_rows, i, begin_block, end_block,
-                                 dof_mask, col_ptr);
+             internals::resolve_matrix_row (global_rows, i, begin_block,
+                                            end_block, dof_mask, col_ptr);
 
              sparsity_pattern.block(block, block_col).add_entries(row,
                                                                   cols.begin(),
index 14bcc926615cb7901f9571e65ecc23d14b1a2ee7..040a854e4dde31c65f41be1db34581eba9e06c63 100644 (file)
@@ -811,13 +811,6 @@ void ConstraintMatrix::clear ()
     lines_cache.swap (tmp);
   }
 
-#ifdef DEAL_II_USE_TRILINOS
-  {
-                                  // reset distribute vector
-    vec_distribute.reset();
-  }
-#endif
-
   sorted = false;
 }
 
@@ -1967,47 +1960,42 @@ ConstraintMatrix::distribute (TrilinosWrappers::MPI::Vector &vec) const
                                   // Here we search all the indices that we
                                   // need to have read-access to - the
                                   // local nodes and all the nodes that the
-                                  // constraints indicate. Do this only at
-                                  // the first call and provide the class
-                                  // with a vector for further use.
-  if (!vec_distribute || vec_distribute->size()!=vec.size())
-    {
-      IndexSet my_indices (vec.size());
+                                  // constraints indicate.
+  IndexSet my_indices (vec.size());
+  {
+    const std::pair<unsigned int, unsigned int>
+      local_range = vec.local_range();
 
-      const std::pair<unsigned int, unsigned int>
-       local_range = vec.local_range();
+    my_indices.add_range (local_range.first, local_range.second);
 
-      my_indices.add_range (local_range.first, local_range.second);
-
-      std::set<unsigned int> individual_indices;
-      for (constraint_iterator it = begin_my_constraints;
-          it != end_my_constraints; ++it)
-       for (unsigned int i=0; i<it->entries.size(); ++i)
-         if ((it->entries[i].first < local_range.first)
-             ||
-             (it->entries[i].first >= local_range.second))
-           individual_indices.insert (it->entries[i].first);
+    std::set<unsigned int> individual_indices;
+    for (constraint_iterator it = begin_my_constraints;
+        it != end_my_constraints; ++it)
+      for (unsigned int i=0; i<it->entries.size(); ++i)
+       if ((it->entries[i].first < local_range.first)
+           ||
+           (it->entries[i].first >= local_range.second))
+         individual_indices.insert (it->entries[i].first);
 
-      my_indices.add_indices (individual_indices.begin(),
-                             individual_indices.end());
+    my_indices.add_indices (individual_indices.begin(),
+                           individual_indices.end());
+  }
 
 #ifdef DEAL_II_COMPILER_SUPPORTS_MPI
-      const Epetra_MpiComm *mpi_comm
-       = dynamic_cast<const Epetra_MpiComm*>(&vec.trilinos_vector().Comm());
+  const Epetra_MpiComm *mpi_comm
+    = dynamic_cast<const Epetra_MpiComm*>(&vec.trilinos_vector().Comm());
 
-      Assert (mpi_comm != 0, ExcInternalError());
+  Assert (mpi_comm != 0, ExcInternalError());
 
-      vec_distribute.reset (new TrilinosWrappers::MPI::Vector
-                           (my_indices.make_trilinos_map (mpi_comm->Comm(),
-                                                          true)));
+  TrilinosWrappers::MPI::Vector vec_distribute
+    (my_indices.make_trilinos_map (mpi_comm->Comm(), true));
 #else
-      vec_distribute.reset (new TrilinosWrappers::MPI::Vector
-                           (my_indices.make_trilinos_map (MPI_COMM_WORLD,
-                                                          true)));
+  TrilinosWrappers::MPI::Vector vec_distribute
+    (my_indices.make_trilinos_map (MPI_COMM_WORLD, true));
 #endif
-    }
+
                                   // here we import the data
-  vec_distribute->reinit(vec,false,true);
+  vec_distribute.reinit(vec,false,true);
 
   for (constraint_iterator it = begin_my_constraints;
        it != end_my_constraints; ++it)
@@ -2017,7 +2005,7 @@ ConstraintMatrix::distribute (TrilinosWrappers::MPI::Vector &vec) const
                                       // different contributions
       double new_value = it->inhomogeneity;
       for (unsigned int i=0; i<it->entries.size(); ++i)
-       new_value += ((*vec_distribute)(it->entries[i].first) *
+       new_value += (vec_distribute(it->entries[i].first) *
                       it->entries[i].second);
       vec(it->line) = new_value;
     }
@@ -2072,16 +2060,16 @@ ConstraintMatrix::distribute (TrilinosWrappers::MPI::BlockVector &vec) const
     }
 
 #ifdef DEAL_II_COMPILER_SUPPORTS_MPI
-      const Epetra_MpiComm *mpi_comm
-       = dynamic_cast<const Epetra_MpiComm*>(&vec.block(0).trilinos_vector().Comm());
+  const Epetra_MpiComm *mpi_comm
+    = dynamic_cast<const Epetra_MpiComm*>(&vec.block(0).trilinos_vector().Comm());
 
-      Assert (mpi_comm != 0, ExcInternalError());
+  Assert (mpi_comm != 0, ExcInternalError());
 
-      TrilinosWrappers::MPI::Vector vec_distribute
-       (my_indices.make_trilinos_map (mpi_comm->Comm(), true));
+  TrilinosWrappers::MPI::Vector vec_distribute
+    (my_indices.make_trilinos_map (mpi_comm->Comm(), true));
 #else
-      TrilinosWrappers::MPI::Vector vec_distribute
-       (my_indices.make_trilinos_map (MPI_COMM_WORLD, true));
+  TrilinosWrappers::MPI::Vector vec_distribute
+    (my_indices.make_trilinos_map (MPI_COMM_WORLD, true));
 #endif
 
                                   // here we import the data
@@ -2446,93 +2434,107 @@ void
 ConstraintMatrix::condense<float>(BlockSparseMatrix<float> &uncondensed) const;
 
 
-#define MATRIX_FUNCTIONS(MatrixType, VectorType)       \
+#define MATRIX_VECTOR_FUNCTIONS(MatrixType, VectorType) \
+template void ConstraintMatrix:: \
+distribute_local_to_global<MatrixType,VectorType > (const FullMatrix<double>        &, \
+                                                    const Vector<double>            &, \
+                                                    const std::vector<unsigned int> &, \
+                                                    MatrixType                      &, \
+                                                    VectorType                      &, \
+                                                    internal::bool2type<false>) const
+#define MATRIX_FUNCTIONS(MatrixType) \
 template void ConstraintMatrix:: \
-distribute_local_to_global<MatrixType > (const FullMatrix<double>        &, \
-                                         const std::vector<unsigned int> &, \
-                                         MatrixType                      &) const; \
+distribute_local_to_global<MatrixType,Vector<double> > (const FullMatrix<double>        &, \
+                                                        const Vector<double>            &, \
+                                                        const std::vector<unsigned int> &, \
+                                                        MatrixType                      &, \
+                                                        Vector<double>                  &, \
+                                                        internal::bool2type<false>) const
+#define BLOCK_MATRIX_VECTOR_FUNCTIONS(MatrixType, VectorType)   \
 template void ConstraintMatrix:: \
 distribute_local_to_global<MatrixType,VectorType > (const FullMatrix<double>        &, \
-                                                   const Vector<double>            &, \
-                                                   const std::vector<unsigned int> &, \
-                                                   MatrixType                      &, \
-                                                   VectorType                      &) const
-
-MATRIX_FUNCTIONS(SparseMatrix<double>, Vector<double>);
-MATRIX_FUNCTIONS(SparseMatrix<float>,  Vector<float>);
-template void ConstraintMatrix::distribute_local_to_global<SparseMatrix<float>,Vector<double> >
-(const FullMatrix<double>        &,
- const Vector<double>            &,
- const std::vector<unsigned int> &,
- SparseMatrix<float>             &,
- Vector<double>                  &) const;
-
-
-MATRIX_FUNCTIONS(BlockSparseMatrix<double>, BlockVector<double>);
-MATRIX_FUNCTIONS(BlockSparseMatrix<float>,  BlockVector<float>);
-template void ConstraintMatrix::distribute_local_to_global<BlockSparseMatrix<float>,BlockVector<double> >
-(const FullMatrix<double>        &,
- const Vector<double>            &,
- const std::vector<unsigned int> &,
- BlockSparseMatrix<float>        &,
- BlockVector<double>             &) const;
-
-
-MATRIX_FUNCTIONS(SparseMatrixEZ<double>, Vector<double>);
-MATRIX_FUNCTIONS(SparseMatrixEZ<float>,  Vector<float>);
-
-// MATRIX_FUNCTIONS(BlockSparseMatrixEZ<double>, Vector<double>);
-// MATRIX_FUNCTIONS(BlockSparseMatrixEZ<float>,  Vector<float>);
-
+                                                    const Vector<double>            &, \
+                                                    const std::vector<unsigned int> &, \
+                                                    MatrixType                      &, \
+                                                    VectorType                      &, \
+                                                    internal::bool2type<true>) const
+#define BLOCK_MATRIX_FUNCTIONS(MatrixType)      \
+template void ConstraintMatrix:: \
+distribute_local_to_global<MatrixType,Vector<double> > (const FullMatrix<double>        &, \
+                                                        const Vector<double>            &, \
+                                                        const std::vector<unsigned int> &, \
+                                                        MatrixType                      &, \
+                                                        Vector<double>                  &, \
+                                                        internal::bool2type<true>) const
+
+MATRIX_FUNCTIONS(SparseMatrix<double>);
+MATRIX_FUNCTIONS(SparseMatrix<float>);
+MATRIX_VECTOR_FUNCTIONS(SparseMatrix<float>, Vector<float>);
+
+BLOCK_MATRIX_FUNCTIONS(BlockSparseMatrix<double>);
+BLOCK_MATRIX_FUNCTIONS(BlockSparseMatrix<float>);
+BLOCK_MATRIX_VECTOR_FUNCTIONS(BlockSparseMatrix<double>, BlockVector<double>);
+BLOCK_MATRIX_VECTOR_FUNCTIONS(BlockSparseMatrix<float>,  BlockVector<float>);
+BLOCK_MATRIX_VECTOR_FUNCTIONS(BlockSparseMatrix<float>,  BlockVector<double>);
+
+MATRIX_FUNCTIONS(SparseMatrixEZ<double>);
+MATRIX_FUNCTIONS(SparseMatrixEZ<float>);
+MATRIX_VECTOR_FUNCTIONS(SparseMatrixEZ<float>,  Vector<float>);
+
+// BLOCK_MATRIX_FUNCTIONS(BlockSparseMatrixEZ<double>);
+// BLOCK_MATRIX_VECTOR_FUNCTIONS(BlockSparseMatrixEZ<float>,  Vector<float>);
 
 #ifdef DEAL_II_USE_PETSC
-MATRIX_FUNCTIONS(PETScWrappers::SparseMatrix, PETScWrappers::Vector);
-MATRIX_FUNCTIONS(PETScWrappers::BlockSparseMatrix, PETScWrappers::BlockVector);
-MATRIX_FUNCTIONS(PETScWrappers::MPI::SparseMatrix, PETScWrappers::MPI::Vector);
-MATRIX_FUNCTIONS(PETScWrappers::MPI::BlockSparseMatrix ,PETScWrappers::MPI::BlockVector);
+MATRIX_FUNCTIONS(PETScWrappers::SparseMatrix);
+BLOCK_MATRIX_FUNCTIONS(PETScWrappers::BlockSparseMatrix);
+MATRIX_FUNCTIONS(PETScWrappers::MPI::SparseMatrix);
+BLOCK_MATRIX_FUNCTIONS(PETScWrappers::MPI::BlockSparseMatrix);
+MATRIX_VECTOR_FUNCTIONS(PETScWrappers::SparseMatrix, PETScWrappers::Vector);
+BLOCK_MATRIX_VECTOR_FUNCTIONS(PETScWrappers::BlockSparseMatrix, PETScWrappers::BlockVector);
+MATRIX_VECTOR_FUNCTIONS(PETScWrappers::MPI::SparseMatrix, PETScWrappers::MPI::Vector);
+BLOCK_MATRIX_VECTOR_FUNCTIONS(PETScWrappers::MPI::BlockSparseMatrix ,PETScWrappers::MPI::BlockVector);
 #endif
 
 #ifdef DEAL_II_USE_TRILINOS
-MATRIX_FUNCTIONS(TrilinosWrappers::SparseMatrix, TrilinosWrappers::Vector);
-MATRIX_FUNCTIONS(TrilinosWrappers::BlockSparseMatrix, TrilinosWrappers::BlockVector);
-template void ConstraintMatrix::distribute_local_to_global
-<TrilinosWrappers::SparseMatrix,TrilinosWrappers::MPI::Vector>
-  (const FullMatrix<double>        &,
-   const Vector<double>            &,
-   const std::vector<unsigned int> &,
-   TrilinosWrappers::SparseMatrix &,
-   TrilinosWrappers::MPI::Vector  &) const;
-template void ConstraintMatrix::distribute_local_to_global
-<TrilinosWrappers::BlockSparseMatrix,TrilinosWrappers::MPI::BlockVector>
-  (const FullMatrix<double>        &,
-   const Vector<double>            &,
-   const std::vector<unsigned int> &,
-   TrilinosWrappers::BlockSparseMatrix &,
-   TrilinosWrappers::MPI::BlockVector  &) const;
+MATRIX_FUNCTIONS(TrilinosWrappers::SparseMatrix);
+BLOCK_MATRIX_FUNCTIONS(TrilinosWrappers::BlockSparseMatrix);
+MATRIX_VECTOR_FUNCTIONS(TrilinosWrappers::SparseMatrix, TrilinosWrappers::Vector);
+BLOCK_MATRIX_VECTOR_FUNCTIONS(TrilinosWrappers::BlockSparseMatrix, TrilinosWrappers::BlockVector);
+MATRIX_VECTOR_FUNCTIONS(TrilinosWrappers::SparseMatrix, TrilinosWrappers::MPI::Vector);
+BLOCK_MATRIX_VECTOR_FUNCTIONS(TrilinosWrappers::BlockSparseMatrix, TrilinosWrappers::MPI::BlockVector);
 #endif
 
 
 #define SPARSITY_FUNCTIONS(SparsityType) \
   template void ConstraintMatrix::add_entries_local_to_global<SparsityType> (\
     const std::vector<unsigned int> &, \
-    SparsityType &,                   \
-    const bool,                               \
-    const Table<2,bool> &) const
+    SparsityType &,                    \
+    const bool,                        \
+    const Table<2,bool> &, \
+    internal::bool2type<false>) const
+#define BLOCK_SPARSITY_FUNCTIONS(SparsityType) \
+  template void ConstraintMatrix::add_entries_local_to_global<SparsityType> (\
+    const std::vector<unsigned int> &, \
+    SparsityType &,                    \
+    const bool,                        \
+    const Table<2,bool> &, \
+    internal::bool2type<true>) const
 
 SPARSITY_FUNCTIONS(SparsityPattern);
 SPARSITY_FUNCTIONS(CompressedSparsityPattern);
 SPARSITY_FUNCTIONS(CompressedSetSparsityPattern);
 SPARSITY_FUNCTIONS(CompressedSimpleSparsityPattern);
-SPARSITY_FUNCTIONS(BlockSparsityPattern);
-SPARSITY_FUNCTIONS(BlockCompressedSparsityPattern);
-SPARSITY_FUNCTIONS(BlockCompressedSetSparsityPattern);
-SPARSITY_FUNCTIONS(BlockCompressedSimpleSparsityPattern);
+BLOCK_SPARSITY_FUNCTIONS(BlockSparsityPattern);
+BLOCK_SPARSITY_FUNCTIONS(BlockCompressedSparsityPattern);
+BLOCK_SPARSITY_FUNCTIONS(BlockCompressedSetSparsityPattern);
+BLOCK_SPARSITY_FUNCTIONS(BlockCompressedSimpleSparsityPattern);
 
 #ifdef DEAL_II_USE_TRILINOS
 SPARSITY_FUNCTIONS(TrilinosWrappers::SparsityPattern);
-SPARSITY_FUNCTIONS(TrilinosWrappers::BlockSparsityPattern);
+BLOCK_SPARSITY_FUNCTIONS(TrilinosWrappers::BlockSparsityPattern);
 #endif
 
+
 #define ONLY_MATRIX_FUNCTIONS(MatrixType) \
   template void ConstraintMatrix::distribute_local_to_global<MatrixType > (\
   const FullMatrix<double>        &, \

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.