hierarchic-approximation.hh

Go to the documentation of this file.

// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set ts=8 sw=4 et sts=4:
 
// SPDX-FileCopyrightText: Copyright © DUNE-FUFEM Project contributors, see file AUTHORS.md
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
 
#ifndef DUNE_FUFEM_GRID_HIERARCHIC_APPROXIMATION_HH
#define DUNE_FUFEM_GRID_HIERARCHIC_APPROXIMATION_HH
 
// Based on
// dune/grid/utility/hierarchicsearch.hh
 
#include <cstddef>
#include <sstream>
#include <string>
#include <utility>
 
#include <dune/common/classname.hh>
#include <dune/common/fvector.hh>
 
#include <dune/grid/common/grid.hh>
#include <dune/grid/common/gridenums.hh>
 
#include <dune/fufem/geometry/polyhedrondistance.hh>
 
template <class Grid, class IS> class HierarchicApproximation {
  enum { dim = Grid::dimension };
 
  enum { dimw = Grid::dimensionworld };
 
  typedef typename Grid::ctype ct;
 
  typedef typename Grid::template Codim<0>::Entity Entity;
 
  typedef typename Grid::HierarchicIterator HierarchicIterator;
 
  template <class Iterator>
  Entity findBestEntity(Iterator first, Iterator last,
                        const Dune::FieldVector<ct, dimw> &global) const {
    double minDistance = std::numeric_limits<double>::infinity();
    Entity ret;
 
    for (Iterator it = first; it != last; ++it) {
      const Entity child = *it;
      const double d = distance(global, child.geometry(), tolerance_);
      if (d < minDistance) {
        ret = child;
        minDistance = d;
      }
    }
    return ret;
  }
 
  Entity hFindEntity(const Entity &entity,
                     const Dune::FieldVector<ct, dimw> &global) const {
    const int childLevel = entity.level() + 1;
    Entity const best = findBestEntity(entity.hbegin(childLevel),
                                       entity.hend(childLevel), global);
    if (indexSet_.contains(best))
      return best;
    else
      return hFindEntity(best, global);
  }
 
public:
  HierarchicApproximation(const Grid &g, const IS &is, double tolerance)
      : grid_(g), indexSet_(is), tolerance_(tolerance) {}
 
  Entity findEntity(const Dune::FieldVector<ct, dimw> &global) const {
    return findEntity<Dune::All_Partition>(global);
  }
 
  template <Dune::PartitionIteratorType partition>
  Entity findEntity(const Dune::FieldVector<ct, dimw> &global) const {
    typedef typename Grid::LevelGridView LevelGV;
    const LevelGV &gv = grid_.template levelGridView(0);
 
    Entity const best = findBestEntity(gv.template begin<0, partition>(),
                                       gv.template end<0, partition>(), global);
    if (indexSet_.contains(best))
      return best;
    else
      return hFindEntity(best, global);
  }
 
private:
  const Grid &grid_;
  const IS &indexSet_;
  const double tolerance_;
};
#endif