distance.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_DISTANCE_HH
#define DUNE_MMESH_REMESHING_DISTANCE_HH
 
#include <dune/common/exceptions.hh>
#include <dune/grid/common/partitionset.hh>
#include <memory>
 
namespace Dune {
 
template <class Grid>
class Distance {
  using ThisType = Distance<Grid>;
  static constexpr int dim = Grid::dimensionworld;
  using ctype = typename Grid::ctype;
  using GlobalCoordinate = FieldVector<ctype, dim>;
  using Vertex = typename Grid::Vertex;
  using Element = typename Grid::template Codim<0>::Entity;
  using Facet = typename Grid::template Codim<1>::Entity;
  using InterfaceElement =
      typename Grid::InterfaceGrid::template Codim<0>::Entity;
 
 public:
  Distance() {}
 
  Distance(const Grid& grid) : grid_(&grid), initialized_(false) {}
 
  void update() {
    // Resize distance_ and set to high default value
    distances_.resize(indexSet().size(dim));
    std::fill(distances_.begin(), distances_.end(), 1e100);
 
    // Set all interface vertices to zero and initialize queue
 
    for (const InterfaceElement& ielement : elements(
             grid_->interfaceGrid().leafGridView(), Partitions::interior)) {
      // Convert to bulk facet
      const Facet facet = grid_->entity(ielement.impl().hostEntity());
 
      // Compute vertex distances to this facet
      handleFacet(facet);
    }
 
    initialized_ = true;
  };
 
  bool initialized() const { return initialized_; }
 
  ctype operator()(const Vertex& vertex) const {
    assert(initialized_);
    assert(indexSet().index(vertex) < size());
    return distances_[indexSet().index(vertex)];
  }
 
  void set(const Vertex& vertex, ctype value) {
    assert(indexSet().index(vertex) < size());
    distances_[indexSet().index(vertex)] = value;
  }
 
  ctype operator()(const Element& element) const {
    assert(initialized_);
    ctype dist = 0.0;
    for (std::size_t i = 0; i < dim + 1; ++i) {
      const auto& vertex = element.template subEntity<dim>(i);
      dist += distances_[indexSet().index(vertex)];
    }
    dist /= dim + 1;
    return dist;
  }
 
  ctype operator()(const InterfaceElement& element) const { return 0.0; }
 
  template <class Index>
  ctype operator[](const Index& index) const {
    assert(initialized_);
    return distances_[index];
  }
 
  ctype maximum() const {
    assert(initialized_);
    double maximum = 0.0;
    for (const auto& d : distances_) maximum = std::max(d, maximum);
    return maximum;
  }
 
  std::size_t size() const {
    assert(initialized_);
    return distances_.size();
  }
 
 private:
  void handleFacet(const Facet& facet) {
    for (const auto& v :
         vertices(grid_->leafGridView(), Partitions::interior)) {
      ctype dist = computeDistance(v, facet);
      ctype& currentDist = distances_[indexSet().index(v)];
      if (dist < currentDist) currentDist = dist;
    }
  }
 
  ctype computeDistance(const Vertex& v, const Facet& facet) const {
    GlobalCoordinate vp = v.geometry().center();
    const auto& geo = facet.geometry();
 
    if constexpr (dim == 2) {
      for (std::size_t i = 0; i < dim; ++i) {
        const GlobalCoordinate& vi = geo.corner(i);
        const GlobalCoordinate& vj = geo.corner(1 - i);
        ctype sgn = (vi - vj) * (vp - vi);
        if (sgn >= 0) return (vp - vi).two_norm();
      }
 
      Dune::Plane<GlobalCoordinate> plane(geo.corner(0), geo.corner(1));
      return std::abs(plane.signedDistance(vp));
    } else  // dim == 3
    {
      // We only use an estimate for 3d.
      ctype dist = 1e100;
      for (std::size_t i = 0; i < dim; ++i) {
        const GlobalCoordinate& vi = geo.corner(i);
        dist = std::min(dist, (vp - vi).two_norm2());
 
        const GlobalCoordinate& vj = geo.corner((i + 1) % dim);
        dist = std::min(dist, (vp - 0.5 * (vi + vj)).two_norm2());
      }
 
      const GlobalCoordinate& vc = geo.center();
      dist = std::min(dist, (vp - vc).two_norm2());
 
      return std::sqrt(dist);
    }
 
    DUNE_THROW(InvalidStateException, "Compute distance failed");
    return -1.;
  }
 
  const typename Grid::LeafIndexSet& indexSet() const {
    return grid_->leafIndexSet();
  }
 
  std::vector<ctype> distances_;
  const Grid* grid_;
  bool initialized_;
};
 
}  // end namespace Dune
 
#endif