marking.hh

#ifndef MARKING_HH
#define MARKING_HH
 
#include "../estimators/estimatorinterface.hh"
 
namespace Dune {
 
namespace ACFem {
 
template<class Estimator>
class MarkingStrategy
{
 public:
  typedef EstimatorInterface<Estimator> EstimatorType;
 protected:
  typedef typename EstimatorType::DiscreteFunctionSpaceType
    DiscreteFunctionSpaceType;
  typedef typename EstimatorType::ElementType ElementType;
  typedef typename DiscreteFunctionSpaceType::IteratorType IteratorType;
 
  const EstimatorType& estimator_;
  const std::string name_;
 public:
  MarkingStrategy(const EstimatorType& estimator, const std::string& name = "adaptation")
    : estimator_(estimator), name_(name)
  {}
 
  bool mark (const double tolerance) const
  {
    // possible strategies
    enum Strategy { none = 0, maximum = 1, equitol = 2, equiv = 3, uniform = 4 };
    static const std::string strategyNames []
      = { "none", "maximum", "equitolerance", "equidistribution", "uniform" };
    Strategy strategy = (Strategy) Dune::Fem::Parameter::getEnum(name_ + "." + "markingStrategy", strategyNames);
 
    double localTol2 = Dune::Fem::Parameter::getValue<double>(name_ + ".refineTolerance", 0);
    double coarseTol2 =  Dune::Fem::Parameter::getValue<double>(name_ + ".coarsenTolerance", std::numeric_limits<double>::max());
    double coarseTolUnif =  Dune::Fem::Parameter::getValue<double>(name_ + ".coarsenToleranceUniform", 0.05);
 
    int refineBisections = 1;
    int coarseBisections = 1;
 
    localTol2 *= localTol2;
    coarseTol2 *= coarseTol2;
 
    switch(strategy) {
    case maximum: {
      double maxError2 = 0;
      // sum up local estimators
      const typename EstimatorType::IteratorType endEstimator = estimator_.end();
      for (typename EstimatorType::IteratorType it = estimator_.begin();
           it != endEstimator; ++it) {
        maxError2 = std::max(*it, maxError2);
      }
 
      // get global maxError
      maxError2 = estimator_.grid().comm().max(maxError2);
      localTol2 *= maxError2;
      coarseTol2 *= maxError2;
 
      if (Dune::Fem::Parameter::verbose()) {
        std::cerr << "max: " << maxError2 << std::endl;
        std::cerr << "local: " << localTol2 << std::endl;
        std::cerr << "coarse: " << coarseTol2 << std::endl;
      }
      break;
    }
 
    case equitol: { // try to equidistribute the given tolerance
      // sum up local estimators
      size_t indSize = estimator_.size();
 
      // get global number of elements
      estimator_.grid().comm().sum(&indSize, 1);
 
      double avg2 = tolerance * tolerance / indSize;
 
      localTol2 *= avg2;
      coarseTol2 *= avg2;
 
      if (Dune::Fem::Parameter::verbose()) {
        std::cerr << "squared mean tolerance: " << avg2 << std::endl;
        std::cerr << "local: " << localTol2 << std::endl;
        std::cerr << "coarse: " << coarseTol2 << std::endl;
      }
      break;
    }
 
    case equiv: { //independent form given tolerance
      double sumError2 = 0;
      // sum up local estimators
      size_t indSize = estimator_.size();
      for (size_t i = 0; i < indSize; ++i)  {
       sumError2 += estimator_[i];
      }
      // get global sum of number of elements and local error sum
      double buffer[2] = { (double)indSize , sumError2 };
      estimator_.grid().comm().sum(buffer, 2);
 
      sumError2 = buffer[1];
      indSize   = buffer[0];
 
      double avg2 = sumError2 / indSize;
 
      localTol2 *= avg2;
      coarseTol2 *= avg2;
 
      if (Dune::Fem::Parameter::verbose()) {
        std::cerr << "estimate: " << std::sqrt(sumError2) << std::endl;
        std::cerr << "squared mean: " << avg2 << std::endl;
        std::cerr << "local: " << localTol2 << std::endl;
        std::cerr << "coarse: " << coarseTol2 << std::endl;
      }
      break;
    }
 
    case uniform: {
      double sumError2 = 0;
      // sum up local estimators
      size_t indSize = estimator_.size();
      for (size_t i=0; i < indSize; ++i)  {
       sumError2 += estimator_[i];
      }
      // get global sum of number of elements and local error sum
      double buffer[2] = { (double)indSize , sumError2 };
      estimator_.grid().comm().sum(buffer, 2);
 
      sumError2 = buffer[1];
      indSize   = buffer[0];
 
      coarseTolUnif *= coarseTolUnif;
 
      refineBisections = coarseBisections = Dune::DGFGridInfo<typename Estimator::GridType>::refineStepsForHalf();
 
      // if global sum > tol, then refine all
      if (sumError2 > tolerance*tolerance ) {
        localTol2 = -1;
        coarseTol2 = -1;
      }
      // if global sum < tol*coarsenToleranceUniform (default 0.05), then coarse all
      else if (sumError2 < tolerance*tolerance*coarseTolUnif ) {
        localTol2 = std::numeric_limits<double>::max();
        coarseTol2 = std::numeric_limits<double>::max();
      }
      // else do nothing
      else {
        localTol2 = std::numeric_limits<double>::max();
        coarseTol2 = -1;
      }
      break;
    }
 
    case none: {
      // settings to do nothing
      localTol2 = std::numeric_limits<double>::max();
      coarseTol2 = -1;
      break;
    }
    default:
      break;
    }
 
    int marked = 0;
    int notmarked = 0;
    int coarseMarked = 0;
    // loop over all elements
    const IteratorType end = estimator_.space().end();
    for (IteratorType it = estimator_.space().begin(); it != end; ++it) {
      const ElementType &entity = *it;
      // check local error indicator
      if (estimator_[entity] > localTol2) {
        // mark entity for refinement
        estimator_.grid().mark(refineBisections, entity);
        //std::cerr << "refineBisections: " << refineBisections << std::endl;
        // grid was marked
        ++marked;
      } else if (estimator_[entity] < coarseTol2) {
        estimator_.grid().mark(-coarseBisections, entity);
        ++coarseMarked;
      } else {
        ++notmarked;
      }
    }
 
    // get globally marked
    int comarray[] = { marked, coarseMarked, notmarked };
    estimator_.grid().comm().sum(comarray, 3);
    marked = comarray[0];
    coarseMarked = comarray[1];
    notmarked = comarray[2];
 
    if (Dune::Fem::Parameter::verbose()) {
      std::cout << "Marked   Elements:   " << marked << std::endl;
      std::cout << "Coarsen Marked:      " << coarseMarked << std::endl;
      std::cout << "Elements not marked: " << notmarked << std::endl;
    }
 
    return bool(marked + coarseMarked);
  }
 
};
 
} // ACFem
} // Dune
 
#endif