The vector containing the coarse grid obstacles is now a member of the...

The vector containing the coarse grid obstacles is now a member of the mmgstep. The only thing that has to be provided are the fine grid obstacles

dune/solvers/iterationsteps/mmgstep.cc

@@ -37,6 +37,17 @@ preprocess()
     // /////////////////////////////////////////////
     //   Set up base solver
     // /////////////////////////////////////////////
+
+    for (size_t i=0; i<obstacleHierarchy_.size()-1; i++)
+        if (obstacleHierarchy_[i])
+            delete(obstacleHierarchy_[i]);
+
+    obstacleHierarchy_.resize(numLevels());
+    obstacleHierarchy_.back() = obstacles_;
+
+    for (size_t i=0; i<obstacleHierarchy_.size()-1; i++)
+        obstacleHierarchy_[i] = new ObstacleVectorType();
+
     if (typeid(*this->basesolver_) == typeid(LoopSolver<VectorType>)) {
 
         LoopSolver<VectorType>* loopBaseSolver = dynamic_cast<LoopSolver<VectorType>* > (this->basesolver_);
@@ -44,14 +55,14 @@ preprocess()
         typedef ProjectedBlockGSStep<MatrixType, VectorType> SmootherType;
         
         dynamic_cast<SmootherType*>(loopBaseSolver->iterationStep_)->hasObstacle_ = hasObstacleHierarchy_[0];
-        dynamic_cast<SmootherType*>(loopBaseSolver->iterationStep_)->obstacles_   = &obstacles_->at(0);
+        dynamic_cast<SmootherType*>(loopBaseSolver->iterationStep_)->obstacles_   = obstacleHierarchy_[0];
         
 #if HAVE_IPOPT
     } else if (typeid(*this->basesolver_) == typeid(QuadraticIPOptSolver<MatrixType,VectorType>)) {
 
         QuadraticIPOptSolver<MatrixType,VectorType>* ipoptBaseSolver = dynamic_cast<QuadraticIPOptSolver<MatrixType,VectorType>*> (this->basesolver_);
      
-        ipoptBaseSolver->obstacles_ = &(*obstacles_)[0];
+        ipoptBaseSolver->obstacles_ = obstacleHierarchy_[0];
 #endif
     } else {
         DUNE_THROW(SolverError, "You can't use " << typeid(*this->basesolver_).name()
@@ -69,12 +80,9 @@ void MonotoneMGStep<MatrixType, VectorType>::nestedIteration()
         // /////////////////////////////
         //   Restrict obstacles
         // //////////////////////////////
-        obstacleRestrictor_->restrict((*obstacles_)[i], 
-                                    (*obstacles_)[i-1],
-                                    *hasObstacleHierarchy_[i], 
-                                    *hasObstacleHierarchy_[i-1],
-                                    *this->mgTransfer_[i-1],
-                                    dummy);
+        obstacleRestrictor_->restrict(*obstacleHierarchy_[i], *obstacleHierarchy_[i-1],
+                                    *hasObstacleHierarchy_[i], *hasObstacleHierarchy_[i-1],
+                                    *this->mgTransfer_[i-1], dummy);
 
         // Restrict right hand side
         this->mgTransfer_[i-1]->restrict(this->rhsHierarchy_[i], this->rhsHierarchy_[i-1]);
@@ -86,7 +94,7 @@ void MonotoneMGStep<MatrixType, VectorType>::nestedIteration()
         // obstacles may be inconsistent.  We do an ad hoc correction here.
         // The true obstacles on the maxlevel are not touched.
         for (size_t i=0; i<(*obstacles_)[this->level_].size(); i++) {
-            BoxConstraint<field_type,dim>& cO = (*obstacles_)[this->level_][i];
+            BoxConstraint<field_type,dim>& cO = (*obstacleHierarchy_[this->level_])[i];
             for (int j=0; j<dim; j++)
                 if (cO.lower(j) > cO.upper(j))
                     cO.lower(j) = cO.upper(j) = (cO.lower(j) + cO.upper(j)) / 2;
@@ -114,7 +122,7 @@ void MonotoneMGStep<MatrixType, VectorType>::iterate()
     std::vector<std::shared_ptr<const MatrixType> >& mat = this->matrixHierarchy_;
     std::vector<std::shared_ptr<VectorType> >& x   = this->xHierarchy_;
     std::vector<VectorType>& rhs = this->rhsHierarchy_;
-    std::vector<std::vector<BoxConstraint<field_type,dim> > >& obstacles = *this->obstacles_;
+    std::vector<ObstacleVectorType* >& obstacles = obstaclesHierarchy_;
 
     Dune::BitSetVector<dim> critical(x[level]->size(), false);
 
@@ -125,12 +133,12 @@ void MonotoneMGStep<MatrixType, VectorType>::iterate()
 
         // Determine critical nodes  (only for debugging)
         const double eps = 1e-12;
-        for (size_t i=0; i<obstacles[level].size(); i++) {
+        for (size_t i=0; i<obstacles[level]->size(); i++) {
             
             for (int j=0; j<dim; j++) {
                 
-                if (obstacles[level][i].lower(j) >= (*x[level])[i][j] - eps
-                    || obstacles[level][i].upper(j) <= (*x[level])[i][j] + eps) {
+                if ((*obstacles[level])[i].lower(j) >= (*x[level])[i][j] - eps
+                    || (*obstacles[level])[i].upper(j) <= (*x[level])[i][j] + eps) {
                     critical[i][j] = true;
                 }
                 
@@ -143,7 +151,7 @@ void MonotoneMGStep<MatrixType, VectorType>::iterate()
         ProjectedBlockGSStep<MatrixType,VectorType>* presmoother = dynamic_cast<ProjectedBlockGSStep<MatrixType, VectorType>*>(this->presmoother_[level].get());
         assert(presmoother);
         presmoother->setProblem(*(mat[level]), *x[level], rhs[level]);
-        presmoother->obstacles_ = &(obstacles[level]);
+        presmoother->obstacles_ = obstacles[level];
         presmoother->hasObstacle_ = hasObstacleHierarchy_[level];
         presmoother->ignoreNodes_ = this->ignoreNodesHierarchy_[level];
         
@@ -151,11 +159,11 @@ void MonotoneMGStep<MatrixType, VectorType>::iterate()
             presmoother->iterate();
         
         // First, a backup of the obstacles for postsmoothing
-        std::vector<BoxConstraint<field_type,dim> > obstacleBackup = obstacles[level];
+        std::vector<BoxConstraint<field_type,dim> > obstacleBackup = *obstacles[level];
         
         // Compute defect obstacles
-        for (size_t i=0; i<obstacles[level].size(); i++)
-            obstacles[level][i] -= (*x[level])[i];
+        for (size_t i=0; i<obstacles[level]->size(); i++)
+            (*obstacles[level])[i] -= (*x[level])[i];
         
         // ///////////////////////
         //    Truncation
@@ -163,12 +171,12 @@ void MonotoneMGStep<MatrixType, VectorType>::iterate()
         
         // Determine critical nodes
         const double eps = 1e-12;
-        for (size_t i=0; i<obstacles[level].size(); i++) {
+        for (size_t i=0; i<obstacles[level]->size(); i++) {
             
             for (int j=0; j<dim; j++) {
                 
-                if (obstacles[level][i].lower(j) >= -eps
-                    || obstacles[level][i].upper(j) <= eps) {
+                if ((*obstacles[level])[i].lower(j) >= -eps
+                    || (*obstacles[level])[i].upper(j) <= eps) {
                     critical[i][j] = true;
                 }
                 
@@ -181,7 +189,7 @@ void MonotoneMGStep<MatrixType, VectorType>::iterate()
             ->galerkinRestrict(*(mat[level]), *(const_cast<MatrixType*>(mat[level-1].get())), critical);
         
         // Restrict obstacles
-        obstacleRestrictor_->restrict(obstacles[level], obstacles[level-1],
+        obstacleRestrictor_->restrict(*obstacles[level], *obstacles[level-1],
                                     *hasObstacleHierarchy_[level], *hasObstacleHierarchy_[level-1],
                                     *this->mgTransfer_[level-1],
                                     critical);
@@ -216,18 +224,18 @@ void MonotoneMGStep<MatrixType, VectorType>::iterate()
         *x[level] += tmp;
         
         // restore the true (non-defect) obstacles
-        obstacles[level] = obstacleBackup;
+        *obstacles[level] = obstacleBackup;
         
 #ifndef NDEBUG
         // Debug: is the current iterate really admissible?
-        for (size_t i=0; i<obstacles[level].size(); i++)
+        for (size_t i=0; i<obstacles[level]->size(); i++)
             for (int j=0; j<VectorType::block_type::dimension; j++)
-                if (((*x[level])[i][j] - obstacles[level][i].lower(j)<-1e-14
-                    || (*x[level])[i][j] - obstacles[level][i].upper(j) >1e-14 )
+                if (((*x[level])[i][j] - (*obstacles[level])[i].lower(j)<-1e-14
+                    || (*x[level])[i][j] - (*obstacles[level])[i].upper(j) >1e-14 )
                     && (!(*this->ignoreNodesHierarchy_[level])[i][j])) {
                     
                     std::cout << "Obstacle disregarded!\n";
-                    std::cout << (*x[level])[i] << std::endl << obstacles[level][i] << std::endl;
+                    std::cout << (*x[level])[i] << std::endl << (*obstacles[level])[i] << std::endl;
                     std::cout << "level: " << level << "   index: " << i << "   komponent: " << j << std::endl;
                     std::cout << "is " << ((critical[i][j]) ? "" : "not") << "critical" << std::endl;
                 }
@@ -237,7 +245,7 @@ void MonotoneMGStep<MatrixType, VectorType>::iterate()
         ProjectedBlockGSStep<MatrixType,VectorType>* postsmoother = dynamic_cast<ProjectedBlockGSStep<MatrixType, VectorType>*>(this->postsmoother_[level].get());
         assert(postsmoother);
         postsmoother->setProblem(*(mat[level]), *x[level], rhs[level]);
-        postsmoother->obstacles_ = &obstacles[level];
+        postsmoother->obstacles_ = obstacles[level];
         postsmoother->hasObstacle_ = hasObstacleHierarchy_[level];
         postsmoother->ignoreNodes_ = this->ignoreNodesHierarchy_[level];
         

dune/solvers/iterationsteps/mmgstep.hh

@@ -17,7 +17,7 @@ class MonotoneMGStep : public MultigridStep<MatrixType, VectorType>
     static const int dim = VectorType::block_type::dimension;
     
     typedef typename VectorType::field_type field_type;
-    
+    typedef std::vector<BoxConstraint<field_type,dim> > ObstacleVectorType;
 public:
     
     MonotoneMGStep() {}
@@ -50,8 +50,9 @@ public:
 
     //! Bitfield determining which fine grid nodes have an obstacle 
     Dune::BitSetVector<1>* hasObstacle_;
-    
-    std::vector<std::vector<BoxConstraint<field_type,dim> > >* obstacles_;
+
+    //! Vector containing the obstacle values of the fine grid nodes
+    ObstacleVectorType * obstacles_;
     
     // Needed to track changes in the set of critical bits, and allows
     // to check which dofs where critical after the last iteration
@@ -60,6 +61,9 @@ public:
 protected:
     //! Bitfield hierarchy containing the coarse obstacle nodes
     std::vector<Dune::BitSetVector<1>* > hasObstacleHierarchy_;
+
+    //! Hierarchy containing the obstacle values of the coarse obstacles
+    std::vector<ObstacleVectorType*>  obstacleHierarchy_;
 };
 
 #include "mmgstep.cc"