Store the hasObstacle information component-wise and not node-wise

dune/solvers/iterationsteps/projectedblockgsstep.cc

@@ -40,7 +40,7 @@ void ProjectedBlockGSStep<OperatorType, DiscFuncType>::iterate()
         VectorBlock& x = (*this->x_)[i];
 
 
-        if ((*hasObstacle_)[i] == false) {
+        if ((*hasObstacle_)[i].none()) {
 
             // No obstacle --> solve linear problem
             mat[i][i].solve(v, r);
@@ -68,10 +68,12 @@ void ProjectedBlockGSStep<OperatorType, DiscFuncType>::iterate()
                     v[k] += sr / mat[i][i][k][k];
 
                     // Project
-                    if (v[k] < defectObstacle.lower(k))
-                        v[k] = defectObstacle.lower(k);
-                    else if (v[k] > defectObstacle.upper(k))
-                        v[k] = defectObstacle.upper(k);
+                    if (hasObstacle_[i][j]==true) {
+                        if (v[k] < defectObstacle.lower(k))
+                            v[k] = defectObstacle.lower(k);
+                        else if (v[k] > defectObstacle.upper(k))
+                            v[k] = defectObstacle.upper(k);
+                    }
 
                 }
 

dune/solvers/iterationsteps/projectedblockgsstep.hh

@@ -31,7 +31,7 @@ public:
     
     virtual DiscFuncType getSol();
     
-    Dune::BitSetVector<1>* hasObstacle_;
+    const Dune::BitSetVector<BlockSize>* hasObstacle_;
     
     const std::vector<BoxConstraint<field_type,BlockSize> >* obstacles_;
 };