diff --git a/dune/solvers/iterationsteps/mmgstep.cc b/dune/solvers/iterationsteps/mmgstep.cc index 4d543d876e508c8995a1d10a38e8f9e023040512..11ae553e4163d28b76d4c611c7e62bccbd05e57b 100644 --- a/dune/solvers/iterationsteps/mmgstep.cc +++ b/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++) diff --git a/dune/solvers/iterationsteps/mmgstep.hh b/dune/solvers/iterationsteps/mmgstep.hh index de32f18c51cd36c260dc0e25eab72763f3a8c3dc..01031cdf695a7b782dd6e3a53732dfe2af8f3967 100644 --- a/dune/solvers/iterationsteps/mmgstep.hh +++ b/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"