patch by Alexander Grayver for condense()

diff --git a/include/deal.II/lac/constraint_matrix.h b/include/deal.II/lac/constraint_matrix.h
index 75b91077b78f036d4d6f49aa05e1fd204e665a2e..282fa05f95154485b46ad06660074565a9faa251 100644 (file)
--- a/include/deal.II/lac/constraint_matrix.h
+++ b/include/deal.II/lac/constraint_matrix.h
@@ -641,26 +641,6 @@ public:
   template <typename number>
   void condense (BlockSparseMatrix<number> &matrix) const;
 
-  /**
-   * Condense the given vector @p uncondensed into @p condensed. It is the
-   * user's responsibility to guarantee that all entries of @p condensed be
-   * zero. Note that this function does not take any inhomogeneity into
-   * account and throws an exception in case there are any
-   * inhomogeneities. Use the function using both a matrix and vector for that
-   * case.
-   *
-   * 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.
-   *
-   * @deprecated The functions converting an uncondensed matrix into
-   * its condensed form are deprecated. Use the functions doing the
-   * in-place condensation leaving the size of the linear system unchanged.
-   */
-  template <class VectorType>
-  void condense (const VectorType &uncondensed,
-                 VectorType       &condensed) const DEAL_II_DEPRECATED;
-
   /**
    * Condense the given vector in-place. The @p VectorType may be a
    * Vector<float>, Vector<double>, BlockVector<tt><...></tt>, a PETSc or
@@ -669,10 +649,23 @@ public:
    * account and throws an exception in case there are any
    * inhomogeneities. Use the function using both a matrix and vector for that
    * case.
+   * Note: this function does not work for MPI vectors because it would require 
+   * ghosted vector in that case, however all ghosted vectors are read-only objects.
+   * Instead, use the analogue with two input arguments.
    */
   template <class VectorType>
   void condense (VectorType &vec) const;
 
+  /**
+   * Same as above, but accepts two input arguments. If called in parallel, one 
+   * vector is supposed to be a ghosted vector and the other vector without ghost
+   * elements. The function copies and condenses values from @p vec_ghosted into 
+   * @p vec.
+   */
+  template <class VectorType>
+  void condense (const VectorType &vec_ghosted,
+                 VectorType       &vec) const;
+
   /**
    * Condense a given matrix and a given vector. The associated matrix struct
    * should be condensed and compressed. It is the user's responsibility to

diff --git a/include/deal.II/lac/constraint_matrix.templates.h b/include/deal.II/lac/constraint_matrix.templates.h
index a1a97f6a953b4100e6e9f479cd1b6e88832fcf69..d54eca03e6bb337daa0eb2d1703c2f0a6a348948 100644 (file)
--- a/include/deal.II/lac/constraint_matrix.templates.h
+++ b/include/deal.II/lac/constraint_matrix.templates.h
@@ -73,82 +73,12 @@ ConstraintMatrix::condense (BlockSparseMatrix<number> &uncondensed) const
 
 template<class VectorType>
 void
-ConstraintMatrix::condense (const VectorType &uncondensed,
-                            VectorType       &condensed) const
+ConstraintMatrix::condense (const VectorType &vec_ghosted,
+                            VectorType       &vec) const
 {
   Assert (sorted == true, ExcMatrixNotClosed());
-  AssertDimension (condensed.size()+n_constraints(), uncondensed.size());
-
-  // store for each line of the
-  // vector its new line number after
-  // compression. If the shift is -1,
-  // this line will be condensed away
-  std::vector<int> new_line;
-
-  new_line.reserve (uncondensed.size());
-
-  std::vector<ConstraintLine>::const_iterator next_constraint = lines.begin();
-  size_type                                   shift           = 0;
-  size_type n_rows = uncondensed.size();
-
-  if (next_constraint == lines.end())
-    // if no constraint is to be handled
-    for (size_type row=0; row!=n_rows; ++row)
-      new_line.push_back (row);
-  else
-    for (size_type row=0; row!=n_rows; ++row)
-      if (row == next_constraint->line)
-        {
-          // this line is constrained
-          new_line.push_back (-1);
-          // note that @p lines is ordered
-          ++shift;
-          ++next_constraint;
-          if (next_constraint == lines.end())
-            // nothing more to do; finish rest
-            // of loop
-            {
-              for (size_type i=row+1; i<n_rows; ++i)
-                new_line.push_back (i-shift);
-              break;
-            };
-        }
-      else
-        new_line.push_back (row-shift);
 
-
-  next_constraint = lines.begin();
-  // note: in this loop we need not check
-  // whether @p next_constraint is a valid
-  // iterator, since @p next_constraint is
-  // only evaluated so often as there are
-  // entries in new_line[*] which tells us
-  // which constraints exist
-  for (size_type row=0; row<uncondensed.size(); ++row)
-    if (new_line[row] != -1)
-      // line not constrained
-      // copy entry
-      condensed(new_line[row]) += uncondensed(row);
-
-    else
-      // line must be distributed
-      {
-        for (size_type q=0; q!=next_constraint->entries.size(); ++q)
-          condensed(new_line[next_constraint->entries[q].first])
-          +=
-            uncondensed(row) * next_constraint->entries[q].second;
-
-        ++next_constraint;
-      };
-}
-
-
-
-template <class VectorType>
-void
-ConstraintMatrix::condense (VectorType &vec) const
-{
-  Assert (sorted == true, ExcMatrixNotClosed());
+  vec = vec_ghosted;
 
   // distribute all entries, and set them to zero. do so in
   // two loops because in the first one we need to add to elements
@@ -166,7 +96,7 @@ ConstraintMatrix::condense (VectorType &vec) const
               ExcMessage ("Inhomogeneous constraint cannot be condensed "
                           "without any matrix specified."));
 
-      const typename VectorType::value_type old_value = vec(constraint_line->line);
+      const typename VectorType::value_type old_value = vec_ghosted(constraint_line->line);
       for (size_type q=0; q!=constraint_line->entries.size(); ++q)
         if (vec.in_local_range(constraint_line->entries[q].first) == true)
           vec(constraint_line->entries[q].first)
@@ -188,6 +118,15 @@ ConstraintMatrix::condense (VectorType &vec) const
 
 
 
+template <class VectorType>
+void
+ConstraintMatrix::condense (VectorType &vec) const
+{
+  condense(vec, vec);
+}
+
+
+
 template<typename number, class VectorType>
 void
 ConstraintMatrix::condense (const SparseMatrix<number> &uncondensed,