From 66ac772dcd0fc2755c4d0795548f7ccc46a6e04b Mon Sep 17 00:00:00 2001
From: Jonathan Youett <youett@math.fu-berlin.de>
Date: Fri, 27 Jun 2014 16:00:02 +0200
Subject: [PATCH] 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 | 62 +++++++++++++++-----------
dune/solvers/iterationsteps/mmgstep.hh | 10 +++--
2 files changed, 42 insertions(+), 30 deletions(-)
diff --git a/dune/solvers/iterationsteps/mmgstep.cc b/dune/solvers/iterationsteps/mmgstep.cc
index 11ae553e..68dfcfdc 100644
--- a/dune/solvers/iterationsteps/mmgstep.cc
+++ b/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];
diff --git a/dune/solvers/iterationsteps/mmgstep.hh b/dune/solvers/iterationsteps/mmgstep.hh
index 01031cdf..b7917a8f 100644
--- a/dune/solvers/iterationsteps/mmgstep.hh
+++ b/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"
--
GitLab