Marker.inc.hpp

#include <cmath>

namespace AMDiS {

template <class Grid>
void Marker<Grid>::mark(Element const& elem, int newMark)
{
  int oldMark = grid_->getMark(elem);

  if (!maximumMarking_ || (newMark > oldMark)) {
    bool marked = grid_->mark(newMark, elem);
    if (marked) {
      if (oldMark > 0) {
        elMarkRefine_--;
      } else if (oldMark < 0) {
        elMarkCoarsen_--;
      }

      if (newMark > 0) {
        elMarkRefine_++;
      } else if (newMark < 0) {
        elMarkCoarsen_++;
      }
    } else {
      msg("Marking failed");
    }
  }
}


template <class Grid>
void Marker<Grid>::initMarking(AdaptInfo& /*adaptInfo*/)
{
  this->elMarkRefine_ = 0;
  this->elMarkCoarsen_ = 0;
}


template <class Grid>
void Marker<Grid>::finishMarking(AdaptInfo& /*adaptInfo*/)
{
  if (info_ > 0) {
    msg("{} elements marked for refinement", elMarkRefine_);
    msg("{} elements marked for coarsening", elMarkCoarsen_);
  }
}


template <class Grid>
Flag Marker<Grid>::markGrid(AdaptInfo& adaptInfo)
{
  test_exit(bool(this->grid_), "No grid!");

  initMarking(adaptInfo);

  if (!this->coarsenAllowed_ && !this->refineAllowed_)
    return 0;

  for (const auto& elem : Dune::elements(this->grid_->leafGridView()))
    markElement(adaptInfo, elem);

  finishMarking(adaptInfo);

  Flag markFlag;
  if (this->elMarkRefine_)
    markFlag = 1;
  if (this->elMarkCoarsen_)
    markFlag |= 2;

  return markFlag;
}

template <class Grid>
void EstimatorMarker<Grid>::initMarking(AdaptInfo& adaptInfo)
{
  Super::initMarking(adaptInfo);
  estSum_ = std::pow(adaptInfo.estSum(component_), p_);
  estMax_ = adaptInfo.estMax(component_);
  this->refineAllowed_ = adaptInfo.isRefinementAllowed(component_);
  this->coarsenAllowed_ = adaptInfo.isCoarseningAllowed(component_);
}


template <class Grid>
void EstimatorMarker<Grid>::markElement(AdaptInfo& /*adaptInfo*/, const Element& elem)
{
  const auto& index = this->grid_->leafIndexSet().index(elem);
  double lError = est_[index];

  if (lError > markRLimit_ && this->refineAllowed_
      && elem.level() < this->maxRefineLevel_) {
    this->mark(elem, 1);
  } else  if (lError <= markCLimit_ && this->coarsenAllowed_
              && elem.level() > this->minRefineLevel_) {
    this->mark(elem, -1);
  }
}


template <class Grid>
std::unique_ptr<EstimatorMarker<Grid>> EstimatorMarker<Grid>::
createMarker(std::string const& name, std::string const& component,
             Estimates const& est, std::shared_ptr<Grid> const& grid)
{
  int strategy = 0;
  Parameters::get(name + "->strategy", strategy);

  switch (strategy) {
  case 0:  // no refinement/coarsening
    break;
  case 1:
    return std::make_unique<GRMarker<Grid> >(name, component, est, grid);
    break;
  case 2:
    return std::make_unique<MSMarker<Grid> >(name, component, est, grid);
    break;
  case 3:
    return std::make_unique<ESMarker<Grid> >(name, component, est, grid);
    break;
  case 4:
    return std::make_unique<GERSMarker<Grid> >(name, component, est, grid);
    break;
  default:
    error_exit("invalid strategy");
  }

  return {};
}


template <class Grid>
void MSMarker<Grid>::initMarking(AdaptInfo& adaptInfo)
{
  Super::initMarking(adaptInfo);

  double msGammaP = std::pow(msGamma_, this->p_);
  double msGammaCP = std::pow(msGammaC_, this->p_);

  this->markRLimit_ = msGammaP * adaptInfo.estMax(this->component_);
  this->markCLimit_ = msGammaCP * adaptInfo.estMax(this->component_);

  msg("start max_est: {}, mark_limits: {}, {}",
    adaptInfo.estMax(this->component_), this->markRLimit_ , this->markCLimit_);
}


template <class Grid>
void ESMarker<Grid>::initMarking(AdaptInfo& adaptInfo)
{
  Super::initMarking(adaptInfo);

  double esThetaP = std::pow(esTheta_, this->p_);
  double esThetaCP = std::pow(esThetaC_, this->p_);
  double epsP = std::pow(adaptInfo.spaceTolerance(this->component_), this->p_);

  int nLeaves = (this->grid_->leafGridView()).size(0);
#if AMDIS_HAS_PARALLEL
  Dune::CollectiveCommunication<>{}.sum(nLeaves, 1);
#endif

  this->markRLimit_ = esThetaP * epsP / nLeaves;
  this->markCLimit_ = esThetaCP * epsP / nLeaves;

  msg("start mark_limits: {}, {}; nt = {}", this->markRLimit_, this->markCLimit_, nLeaves);
}


template <class Grid>
Flag GERSMarker<Grid>::markGrid(AdaptInfo& adaptInfo)
{
  Super::initMarking(adaptInfo);

  if (!this->coarsenAllowed_ && !this->refineAllowed_)
    return 0;

  gersSum_ = 0.0;

  double LTheta = std::pow(1.0 - gersThetaStar_, this->p_);
  double epsP = std::pow(adaptInfo.spaceTolerance(this->component_), this->p_);

  if (this->estSum_ < oldErrSum_) {
    double improv = this->estSum_ / oldErrSum_;
    double wanted = 0.8 * epsP / this->estSum_;
    double redfac = std::min((1.0 - wanted) / (1.0 - improv), 1.0);
    redfac = std::max(redfac, 0.0);

    if (redfac < 1.0) {
      LTheta *= redfac;
      msg("GERS: use extrapolated theta_star = {}", std::pow(LTheta, 1.0 / this->p_));
    }
  }

  oldErrSum_ = this->estSum_;
  double GERSGamma = 1.0;

  if (this->refineAllowed_) {
    if (LTheta > 0) {
      do {
        gersSum_ = 0.0;
        GERSGamma -= gersNu_;
        this->markRLimit_ = GERSGamma * this->estMax_;

        for (const auto& elem : Dune::elements(this->grid_->leafGridView()))
          markElementForRefinement(adaptInfo, elem);

      } while((GERSGamma > 0) && (gersSum_ < LTheta * this->estSum_));
    }

    msg("GERS refinement with gamma = {}", GERSGamma);
  }

  if (this->coarsenAllowed_) {
    GERSGamma = 0.3;
    LTheta = gersThetaC_ * epsP;

    do {
      gersSum_ = 0.0;
      GERSGamma -= gersNu_;
      this->markCLimit_ = GERSGamma * this->estMax_;

      for (const auto& elem : Dune::elements(this->grid_->leafGridView()))
        markElementForCoarsening(adaptInfo, elem);

      msg("coarse loop: gamma = {}, sum = {}, limit = {}", GERSGamma, gersSum_, LTheta);
    } while(gersSum_ > LTheta);

    msg("GERS coarsening with gamma = {}", GERSGamma);
  }

  Super::finishMarking(adaptInfo);

  Flag markFlag;
  if (this->elMarkRefine_)
    markFlag = 1;
  if (this->elMarkCoarsen_)
    markFlag |= 2;

  return markFlag;
}


template <class Grid>
void GERSMarker<Grid>::markElementForRefinement(AdaptInfo& /*adaptInfo*/, const Element& elem)
{
  double lError = this->est_[this->grid_->leafIndexSet().index(elem)];

  if (lError > this->markRLimit_) {
    gersSum_ += lError;
    this->mark(elem, 1);
  }
}


template <class Grid>
void GERSMarker<Grid>::markElementForCoarsening(AdaptInfo& /*adaptInfo*/, const Element& elem)
{
  double lError = this->est_[this->grid_->leafIndexSet().index(elem)];

  if (this->grid_->getMark(elem) <= 0) {
    if (lError <= this->markCLimit_) {
      gersSum_ += lError;
      this->mark(elem, -1);
    } else {
      this->mark(elem, 0);
    }
  }
}


template <class Grid, class PreGridFct>
Flag GridFunctionMarker<Grid, PreGridFct>::markGrid(AdaptInfo& adaptInfo)
{
  test_exit(bool(this->grid_), "No grid!");

  Super::initMarking(adaptInfo);

  if (!this->coarsenAllowed_ && !this->refineAllowed_)
    return 0;

  auto localFct = localFunction(gridFct_);

  for (auto const& e : Dune::elements(this->grid_->leafGridView())) {
    localFct.bind(e);
    int currentLevel = e.level();
    auto refElem = Dune::referenceElement<typename Grid::ctype,Grid::dimension>(e.type());

    // evaluate in the center of the element
    int targetLevel = int(std::round(localFct(refElem.position(0,0))));

    int m =   ((((targetLevel > currentLevel) && (currentLevel < this->maxRefineLevel_))
                  || (currentLevel < this->minRefineLevel_))
                && this->refineAllowed_)
            - ((((targetLevel < currentLevel) && (currentLevel > this->minRefineLevel_))
                  || (currentLevel > this->maxRefineLevel_))
                && this->coarsenAllowed_);
    this->mark(e, m);
    localFct.unbind();
  }

  Super::finishMarking(adaptInfo);

  Flag markFlag;
  if (this->elMarkRefine_)
    markFlag = 1;
  if (this->elMarkCoarsen_)
    markFlag |= 2;

  return markFlag;
}

} // end namespace AMDiS