From 406d8680bf6dafba525496165ef85288ae475304 Mon Sep 17 00:00:00 2001
From: Timo Heister <timo.heister@gmail.com>
Date: Sun, 7 Jul 2013 15:47:02 +0000
Subject: [PATCH] fix bug in ConstraintMatrix::set_zero

git-svn-id: https://svn.dealii.org/trunk@29940 0785d39b-7218-0410-832d-ea1e28bc413d
---
 .../deal.II/lac/constraint_matrix.templates.h    | 16 +++++++++++++---
 1 file changed, 13 insertions(+), 3 deletions(-)

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 397be1b7cc..a8c1d393c1 100644
--- a/deal.II/include/deal.II/lac/constraint_matrix.templates.h
+++ b/deal.II/include/deal.II/lac/constraint_matrix.templates.h
@@ -738,12 +738,22 @@ namespace internal
       template<class VEC>
       void set_zero_parallel(const std::vector<size_type> &cm, VEC &vec, size_type shift = 0)
       {
-        Assert(!vec.has_ghost_elements(), ExcInternalError());//ExcGhostsPresent());
+        Assert(!vec.has_ghost_elements(), ExcInternalError());
         IndexSet locally_owned = vec.locally_owned_elements();
         for (typename std::vector<size_type>::const_iterator it = cm.begin();
              it != cm.end(); ++it)
-          if (locally_owned.is_element(*it))
-            vec(*it) = 0.;
+          {
+            // If shift>0 then we are working on a part of a BlockVector
+            // so vec(i) is actually the global entry i+shift.
+            // We first make sure the line falls into the range of vec,
+            // then check if is part of the local part of the vector, before
+            // finally setting it to 0.
+            if ((*it)<shift)
+              continue;
+            size_type idx = *it - shift;
+            if (idx<vec.size() && locally_owned.is_element(idx))
+              vec(idx) = 0.;
+          }
       }
 
       template<typename Number>
-- 
2.39.5