diff --git a/dune/solvers/iterationsteps/multigridstep.cc b/dune/solvers/iterationsteps/multigridstep.cc index a8d97ed04dc9656b942711448ef227ec79360c9f..f58eed599052a720ff0dea42a3304755b2f078ed 100644 --- a/dune/solvers/iterationsteps/multigridstep.cc +++ b/dune/solvers/iterationsteps/multigridstep.cc @@ -146,49 +146,52 @@ void MultigridStep<MatrixType, VectorType, BitVectorType>::preprocess() // Set up base solver // ///////////////////////////////////////////// - if (this->basesolver_ == NULL) - DUNE_THROW(SolverError, "You have not provided a base solver!"); + if (basesolver_) + { + if (this->basesolver_ == NULL) + DUNE_THROW(SolverError, "You have not provided a base solver!"); - // If the base solver can ignore dofs give it the ignoreNodes field - if (dynamic_cast<CanIgnore<BitVectorType>*>(this->basesolver_)) - dynamic_cast<CanIgnore<BitVectorType>*>(this->basesolver_)->ignoreNodes_ = ignoreNodesHierarchy_[0]; + // If the base solver can ignore dofs give it the ignoreNodes field + if (dynamic_cast<CanIgnore<BitVectorType>*>(this->basesolver_)) + dynamic_cast<CanIgnore<BitVectorType>*>(this->basesolver_)->ignoreNodes_ = ignoreNodesHierarchy_[0]; - typedef ::LoopSolver<VectorType> DuneSolversLoopSolver; + typedef ::LoopSolver<VectorType> DuneSolversLoopSolver; - if (typeid(*this->basesolver_) == typeid(DuneSolversLoopSolver)) { + if (typeid(*this->basesolver_) == typeid(DuneSolversLoopSolver)) { - DuneSolversLoopSolver* loopBaseSolver = dynamic_cast<DuneSolversLoopSolver*> (this->basesolver_); + DuneSolversLoopSolver* loopBaseSolver = dynamic_cast<DuneSolversLoopSolver*> (this->basesolver_); - typedef LinearIterationStep<MatrixType, VectorType> SmootherType; - assert(dynamic_cast<SmootherType*>(loopBaseSolver->iterationStep_)); + typedef LinearIterationStep<MatrixType, VectorType> SmootherType; + assert(dynamic_cast<SmootherType*>(loopBaseSolver->iterationStep_)); - dynamic_cast<SmootherType*>(loopBaseSolver->iterationStep_)->setProblem(*(this->matrixHierarchy_[0]), *this->xHierarchy_[0], this->rhsHierarchy_[0]); - dynamic_cast<SmootherType*>(loopBaseSolver->iterationStep_)->ignoreNodes_ = ignoreNodesHierarchy_[0]; + dynamic_cast<SmootherType*>(loopBaseSolver->iterationStep_)->setProblem(*(this->matrixHierarchy_[0]), *this->xHierarchy_[0], this->rhsHierarchy_[0]); + dynamic_cast<SmootherType*>(loopBaseSolver->iterationStep_)->ignoreNodes_ = ignoreNodesHierarchy_[0]; - } - // TODO: The following two #if-clauses do the exactly the same thing: they call setProblem - // However, we cannot write the code generically because there is no general abstraction - // for solvers that I can call 'setProblem' for. + } + // TODO: The following two #if-clauses do the exactly the same thing: they call setProblem + // However, we cannot write the code generically because there is no general abstraction + // for solvers that I can call 'setProblem' for. #if HAVE_IPOPT - else if (typeid(*this->basesolver_) == typeid(QuadraticIPOptSolver<MatrixType,VectorType>)) { + else if (typeid(*this->basesolver_) == typeid(QuadraticIPOptSolver<MatrixType,VectorType>)) { - QuadraticIPOptSolver<MatrixType,VectorType>* ipoptBaseSolver = dynamic_cast<QuadraticIPOptSolver<MatrixType,VectorType>*> (this->basesolver_); + QuadraticIPOptSolver<MatrixType,VectorType>* ipoptBaseSolver = dynamic_cast<QuadraticIPOptSolver<MatrixType,VectorType>*> (this->basesolver_); - ipoptBaseSolver->setProblem(*(this->matrixHierarchy_[0]), *this->xHierarchy_[0], this->rhsHierarchy_[0]); - } + ipoptBaseSolver->setProblem(*(this->matrixHierarchy_[0]), *this->xHierarchy_[0], this->rhsHierarchy_[0]); + } #endif #if HAVE_UMFPACK - else if (typeid(*this->basesolver_) == typeid(Dune::Solvers::UMFPackSolver<MatrixType,VectorType>)) { + else if (typeid(*this->basesolver_) == typeid(Dune::Solvers::UMFPackSolver<MatrixType,VectorType>)) { - Dune::Solvers::UMFPackSolver<MatrixType,VectorType>* umfpackBaseSolver - = dynamic_cast<Dune::Solvers::UMFPackSolver<MatrixType,VectorType>*> (this->basesolver_); + Dune::Solvers::UMFPackSolver<MatrixType,VectorType>* umfpackBaseSolver + = dynamic_cast<Dune::Solvers::UMFPackSolver<MatrixType,VectorType>*> (this->basesolver_); - umfpackBaseSolver->setProblem(*(this->matrixHierarchy_[0]), *this->xHierarchy_[0], this->rhsHierarchy_[0]); - } + umfpackBaseSolver->setProblem(*(this->matrixHierarchy_[0]), *this->xHierarchy_[0], this->rhsHierarchy_[0]); + } #endif -else { - DUNE_THROW(SolverError, "You can't use " << typeid(*this->basesolver_).name() - << " as a base solver in a MultigridStep!"); + else { + DUNE_THROW(SolverError, "You can't use " << typeid(*this->basesolver_).name() + << " as a base solver in a MultigridStep!"); + } } } @@ -226,7 +229,7 @@ void MultigridStep<MatrixType, VectorType, BitVectorType>::iterate() std::vector<VectorType>& rhs = this->rhsHierarchy_; // Solve directly if we're looking at the coarse problem - if (level == 0) { + if ((level == 0) and (basesolver_)) { basesolver_->solve(); return; } @@ -239,38 +242,41 @@ void MultigridStep<MatrixType, VectorType, BitVectorType>::iterate() for (int i=0; i<this->nu1_; i++) presmoother_[level]->iterate(); - // ///////////////////////// - // Restriction + if (level>0) + { + // ///////////////////////// + // Restriction - // compute residual - // fineResidual = rhs[level] - mat[level] * x[level]; - VectorType fineResidual = rhs[level]; - mat[level]->mmv(*x[level], fineResidual); + // compute residual + // fineResidual = rhs[level] - mat[level] * x[level]; + VectorType fineResidual = rhs[level]; + mat[level]->mmv(*x[level], fineResidual); - // restrict residual - this->mgTransfer_[level-1]->restrict(fineResidual, rhs[level-1]); + // restrict residual + this->mgTransfer_[level-1]->restrict(fineResidual, rhs[level-1]); - // Set Dirichlet values. - GenericVector::truncate(rhs[level-1], *ignoreNodesHierarchy_[level-1]); + // Set Dirichlet values. + GenericVector::truncate(rhs[level-1], *ignoreNodesHierarchy_[level-1]); - // Choose all zeros as the initial correction - *x[level-1] = 0; + // Choose all zeros as the initial correction + *x[level-1] = 0; - // /////////////////////////////////////// - // Recursively solve the coarser system - level--; - for (int i=0; i<mu_; i++) - iterate(); - level++; + // /////////////////////////////////////// + // Recursively solve the coarser system + level--; + for (int i=0; i<mu_; i++) + iterate(); + level++; - // //////////////////////////////////////// - // Prolong + // //////////////////////////////////////// + // Prolong - // add correction to the presmoothed solution - VectorType tmp; - this->mgTransfer_[level-1]->prolong(*x[level-1], tmp); - *x[level] += tmp; + // add correction to the presmoothed solution + VectorType tmp; + this->mgTransfer_[level-1]->prolong(*x[level-1], tmp); + *x[level] += tmp; + } // Postsmoothing postsmoother_[level]->setProblem(*(mat[level]), *x[level], rhs[level]);