Fix a problem where we needed to call compress() in ConstraintMatrix::condense.

git-svn-id: https://svn.dealii.org/trunk@27996 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h
index 0c145932633f01e3ccf7541a05d3339e8152a028..3fae5d6b7867c366b989e90721bc2ad4e5756428 100644 (file)
--- a/deal.II/doc/news/changes.h
+++ b/deal.II/doc/news/changes.h
@@ -172,6 +172,11 @@ DoFHandler, in particular removal of specializations.
 <h3>Specific improvements</h3>
 
 <ol>
+<li> Fixed: The one-argument version of ConstraintMatrix::condense was
+not prepared to deal with parallel vectors. This is now fixed.
+<br>
+(Daniel Arndt, Wolfgang Bangerth, 2013/1/9)
+
 <li> New: Utilities::int_to_string can now deal with any 32-bit
 integer, not just those with up to 6 digits.
 <br>

diff --git a/deal.II/include/deal.II/lac/constraint_matrix.templates.h b/deal.II/include/deal.II/lac/constraint_matrix.templates.h
index 8d7f1b552a77af9e01e6571dd1fdbd974de96305..d1ad4d3565d5b8c4ef7b7c15469c03614fef9a96 100644 (file)
--- a/deal.II/include/deal.II/lac/constraint_matrix.templates.h
+++ b/deal.II/include/deal.II/lac/constraint_matrix.templates.h
@@ -1,7 +1,7 @@
 //---------------------------------------------------------------------------
 //    $Id$
 //
-//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2012 by the deal.II authors
+//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2012, 2013 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -142,24 +142,36 @@ ConstraintMatrix::condense (VectorType &vec) const
 {
   Assert (sorted == true, ExcMatrixNotClosed());
 
-  // distribute all entries, and set them to zero
-  std::vector<ConstraintLine>::const_iterator constraint_line = lines.begin();
-  for (; constraint_line!=lines.end(); ++constraint_line)
+  // distribute all entries, and set them to zero. do so in
+  // two loops because in the first one we need to add to elements
+  // and in the second one we need to set elements to zero. for
+  // parallel vectors, this can only work if we can put a compress()
+  // in between, but we don't want to call compress() twice per entry
+  for (std::vector<ConstraintLine>::const_iterator
+        constraint_line = lines.begin();
+       constraint_line!=lines.end(); ++constraint_line)
     {
-      for (unsigned int q=0; q!=constraint_line->entries.size(); ++q)
-        vec(constraint_line->entries[q].first)
-        += (static_cast<typename VectorType::value_type>
-            (vec(constraint_line->line)) *
-            constraint_line->entries[q].second);
-      vec(constraint_line->line) = 0.;
-
       // in case the constraint is
       // inhomogeneous, this function is not
       // appropriate. Throw an exception.
       Assert (constraint_line->inhomogeneity == 0.,
               ExcMessage ("Inhomogeneous constraint cannot be condensed "
                           "without any matrix specified."));
+
+      const typename VectorType::value_type old_value = vec(constraint_line->line);
+      for (unsigned int q=0; q!=constraint_line->entries.size(); ++q)
+        vec(constraint_line->entries[q].first)
+        += (static_cast<typename VectorType::value_type>
+            (old_value) *
+            constraint_line->entries[q].second);
     }
+
+  vec.compress();
+
+  for (std::vector<ConstraintLine>::const_iterator
+        constraint_line = lines.begin();
+       constraint_line!=lines.end(); ++constraint_line)
+    vec(constraint_line->line) = 0.;
 }