ratioindicator.hh

Go to the documentation of this file.

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_MMESH_REMESHING_RATIOINDICATOR_HH
#define DUNE_MMESH_REMESHING_RATIOINDICATOR_HH
 
#include <dune/common/exceptions.hh>
#include <dune/grid/common/partitionset.hh>
#include <dune/mmesh/remeshing/distance.hh>
#include <memory>
 
namespace Dune {
 
template <class Grid>
class RatioIndicator {
  static constexpr int dim = Grid::dimensionworld;
  using ctype = typename Grid::ctype;
  using GlobalCoordinate = FieldVector<ctype, dim>;
 
 public:
  using DistanceType = Distance<Grid>;
 
  RatioIndicator(ctype h = 0.0, ctype distProportion = 1.0, ctype factor = 1.0)
      : edgeRatio_(4.),  // the edge ratio. Decreases minH!
        K_(2.),  // the factor for the maximal edge length. Defaults to 2,
                 // because we use edge bisection.
        maxH_(K_ * h),
        k_(K_ / (2. * edgeRatio_)),  // ensure that a triangle with two edges
                                     // longer than maxH splits up to a triangle
                                     // where all edges are longer than minH
        minH_(k_ * h),
        radiusRatio_(30.),  // additional radius ratio to avoid very ugly cells
        distProportion_(distProportion),  // cells with distance to interface of
                                          // value greater than distProportion_
                                          // * max(dist) are refined to ...
        factor_(factor)  // ... edge length in [factor * minH_, factor * maxH_]
  {}
 
  void init(const Grid& grid) {
    maxH_ = 0.0;
    minH_ = 1e100;
 
    for (const auto& edge : edges(grid.interfaceGrid().leafGridView())) {
      const ctype h = edge.geometry().volume();
      maxH_ = std::max(maxH_, h);
      minH_ = std::min(minH_, h);
    }
 
    maxH_ = K_ * maxH_;
    minH_ = k_ * minH_;
 
    ctype maxh = 0.0;
    ctype minh = 1e100;
    for (const auto& edge : edges(grid.leafGridView())) {
      const ctype h = edge.geometry().volume();
      maxh = std::max(maxh, h);
      minh = std::min(minh, h);
    }
 
    // fallback if there is no interface
    if (grid.interfaceGrid().size(1) == 0) {
      maxH_ = K_ * maxh;
      minH_ = k_ * minh;
    }
 
    factor_ = maxh / minh;
    distance_ = DistanceType(grid);
  };
 
  void update() {
    distance_.update();
    maxDist_ = distProportion_ * distance_.maximum();
  }
 
  template <class Element>
  int operator()(const Element& element) const {
    static constexpr int edgeCodim = Element::dimension - 1;
 
    int refine = 0;
    int coarse = 0;
 
    ctype minE = 1e100;
    ctype maxE = -1e100;
    ctype sumE = 0.0;
 
    // edge length
    ctype dist = std::min(maxDist_, distance_(element));
    const ctype l = dist / maxDist_;
    const ctype minH = (1. - l) * minH_ + l * factor_ * minH_;
    const ctype maxH = (1. - l) * maxH_ + l * factor_ * maxH_;
 
    for (std::size_t i = 0; i < element.subEntities(edgeCodim); ++i) {
      const auto& edge = element.template subEntity<edgeCodim>(i);
      const ctype len = edge.geometry().volume();
 
      if (len < minH) coarse++;
 
      if (len > maxH) refine++;
 
      minE = std::min(minE, len);
      maxE = std::max(maxE, len);
 
      sumE += len;
    }
 
    // edge ratio criterion
    const ctype edgeRatio = maxE / minE;
    if (edgeRatio > edgeRatio_) coarse++;
 
    // radius ratio criterion
    const auto& geo = element.geometry();
    const ctype innerRadius = dim * geo.volume() / sumE;
    const ctype outerRadius =
        (geo.impl().circumcenter() - geo.corner(0)).two_norm();
    const ctype radiusRatio = outerRadius / (innerRadius * dim);
 
    if (radiusRatio > radiusRatio_) coarse++;
 
    // priority on coarse
    if (coarse > 0)
      if (dim != 3)  // disable coarsening in 3d until it is implemented
        return -1;
 
    // then refine
    if (refine > 0) return 1;
 
    // nothing to do
    return 0;
  }
 
  ctype edgeRatio() const { return edgeRatio_; }
 
  ctype maxH() const { return maxH_; }
 
  ctype& maxH() { return maxH_; }
 
  ctype minH() const { return minH_; }
 
  ctype& minH() { return minH_; }
 
  ctype radiusRatio() const { return radiusRatio_; }
 
  ctype& distProportion() { return distProportion_; }
 
  ctype& factor() { return factor_; }
 
  const DistanceType& distance() const {
    if (!distance_.initialized()) distance_.update();
    return distance_;
  }
 
 private:
  const ctype edgeRatio_;
  const ctype K_;
  ctype maxH_;
  const ctype k_;
  ctype minH_;
  const ctype radiusRatio_;
  ctype distProportion_;
  ctype factor_;
  ctype maxDist_;
  mutable DistanceType distance_;
};
 
}  // end namespace Dune
 
#endif