The hasObstacle hierarchy is now a member that is assembled during the...

The hasObstacle hierarchy is now a member that is assembled during the preprocess step. The only information provided from the outside is a bitfield determining the fine grid obstacles

dune/solvers/iterationsteps/mmgstep.cc

@@ -15,16 +15,24 @@ preprocess()
     // Call preprocess of the base class
     MultigridStep<MatrixType,VectorType>::preprocess();
 
-    // Then specify the subset of entries to be reassembled at each iteration
+    // Then setup the obstacle flags which specify the subset of entries to be reassembled at each iteration
+     for (size_t i=0; i<hasObstacleHierarchy_.size()-1; i++)
+        if (hasObstacleHierarchy_[i])
+            delete(hasObstacleHierarchy_[i]);
+
+    hasObstacleHierarchy_.resize(numLevels());
+    hasObstacleHierarchy_.back() = hasObstacle_;
+
     for (int i=this->mgTransfer_.size()-1; i>=0; i--) {
-        this->mgTransfer_[i]->restrictScalarBitField(hasObstacle_->at(i+1), hasObstacle_->at(i));
-        for (size_t j=0; j<hasObstacle_->at(i).size(); j++)
+        hasObstacleHierarchy_[i] = new Dune::BitSetVector<1>();
+        this->mgTransfer_[i]->restrictScalarBitField(*hasObstacleHierarchy_[i+1], *hasObstacleHierarchy_[i]);
+        for (size_t j=0; j<hasObstacleHierarchy_[i]->size(); j++)
             if ((*this->ignoreNodesHierarchy_[i])[j].any())
-                hasObstacle_->at(i)[j][0]=false;
+                (*hasObstacleHierarchy_[i])[j][0]=false;
     }
 
     for (size_t i=0; i<this->mgTransfer_.size(); i++)
-        dynamic_cast<TruncatedMGTransfer<VectorType>*>(this->mgTransfer_[i])->recompute_ = &hasObstacle_->at(i);
+        dynamic_cast<TruncatedMGTransfer<VectorType>*>(this->mgTransfer_[i])->recompute_ = hasObstacleHierarchy_[i];
     
     // /////////////////////////////////////////////
     //   Set up base solver
@@ -35,7 +43,7 @@ preprocess()
      
         typedef ProjectedBlockGSStep<MatrixType, VectorType> SmootherType;
         
-        dynamic_cast<SmootherType*>(loopBaseSolver->iterationStep_)->hasObstacle_ = &hasObstacle_->at(0);
+        dynamic_cast<SmootherType*>(loopBaseSolver->iterationStep_)->hasObstacle_ = hasObstacleHierarchy_[0];
         dynamic_cast<SmootherType*>(loopBaseSolver->iterationStep_)->obstacles_   = &obstacles_->at(0);
         
 #if HAVE_IPOPT
@@ -63,8 +71,8 @@ void MonotoneMGStep<MatrixType, VectorType>::nestedIteration()
         // //////////////////////////////
         obstacleRestrictor_->restrict((*obstacles_)[i], 
                                     (*obstacles_)[i-1],
-                                    (*hasObstacle_)[i], 
-                                    (*hasObstacle_)[i-1],
+                                    *hasObstacleHierarchy_[i], 
+                                    *hasObstacleHierarchy_[i-1],
                                     *this->mgTransfer_[i-1],
                                     dummy);
 
@@ -136,7 +144,7 @@ void MonotoneMGStep<MatrixType, VectorType>::iterate()
         assert(presmoother);
         presmoother->setProblem(*(mat[level]), *x[level], rhs[level]);
         presmoother->obstacles_ = &(obstacles[level]);
-        presmoother->hasObstacle_ = &((*hasObstacle_)[level]);
+        presmoother->hasObstacle_ = hasObstacleHierarchy_[level];
         presmoother->ignoreNodes_ = this->ignoreNodesHierarchy_[level];
         
         for (int i=0; i<this->nu1_; i++)
@@ -174,7 +182,7 @@ void MonotoneMGStep<MatrixType, VectorType>::iterate()
         
         // Restrict obstacles
         obstacleRestrictor_->restrict(obstacles[level], obstacles[level-1],
-                                    (*hasObstacle_)[level], (*hasObstacle_)[level-1],
+                                    *hasObstacleHierarchy_[level], *hasObstacleHierarchy_[level-1],
                                     *this->mgTransfer_[level-1],
                                     critical);
         
@@ -230,7 +238,7 @@ void MonotoneMGStep<MatrixType, VectorType>::iterate()
         assert(postsmoother);
         postsmoother->setProblem(*(mat[level]), *x[level], rhs[level]);
         postsmoother->obstacles_ = &obstacles[level];
-        postsmoother->hasObstacle_ = &((*hasObstacle_)[level]);
+        postsmoother->hasObstacle_ = hasObstacleHierarchy_[level];
         postsmoother->ignoreNodes_ = this->ignoreNodesHierarchy_[level];
         
         for (int i=0; i<this->nu2_; i++)

dune/solvers/iterationsteps/mmgstep.hh

@@ -47,14 +47,19 @@ public:
     virtual void nestedIteration();
     
     ObstacleRestrictor<VectorType>* obstacleRestrictor_;
-    
-    std::vector<Dune::BitSetVector<1> >* hasObstacle_;
+
+    //! Bitfield determining which fine grid nodes have an obstacle 
+    Dune::BitSetVector<1>* hasObstacle_;
     
     std::vector<std::vector<BoxConstraint<field_type,dim> > >* obstacles_;
     
     // Needed to track changes in the set of critical bits, and allows
     // to check which dofs where critical after the last iteration
     Dune::BitSetVector<dim> oldCritical;
+
+protected:
+    //! Bitfield hierarchy containing the coarse obstacle nodes
+    std::vector<Dune::BitSetVector<1>* > hasObstacleHierarchy_;
 };
 
 #include "mmgstep.cc"