Add hierarchic search for a/the closest element

dune/fufem/CMakeLists.txt

@@ -34,6 +34,7 @@ install(FILES
     facelocalfiniteelement.hh
     formatstring.hh
     globalintersectioniterator.hh
+    hierarchic-approximation.hh
     improvegrid.hh
     interval.hh
     lumpmatrix.hh

dune/fufem/hierarchic-approximation.hh

+#ifndef DUNE_FUFEM_HIERARCHIC_APPROXIMATION_HH
+#define DUNE_FUFEM_HIERARCHIC_APPROXIMATION_HH
+
+// Based on
+// dune/grid/utility/hierarchicsearch.hh
+
+/**
+   @file
+   @brief Utility class for hierarchically searching for the Entity
+   closest to a given point.
+ */
+
+#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/polyhedrondistance.hh>
+
+/**
+   @brief Search an IndexSet for the Entity closest to a given point.
+ */
+template <class Grid, class IS> class HierarchicApproximation {
+  //! get dimension from the grid
+  enum { dim = Grid::dimension };
+
+  //! get world dimension from the grid
+  enum { dimw = Grid::dimensionworld };
+
+  //! get coord type from the grid
+  typedef typename Grid::ctype ct;
+
+  //! get entity from the grid
+  typedef typename Grid::template Codim<0>::Entity Entity;
+
+  //! type of EntityPointer
+  typedef typename Grid::template Codim<0>::EntityPointer EntityPointer;
+
+  //! type of HierarchicIterator
+  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;
+  }
+
+  /**
+     internal helper method
+
+     @param[in] entity Entity whos children should be searched
+     @param[in] global Point you are searching for
+
+     Search the child entity closest to the point global. Recursively
+     continue until we found an entity that is part of
+     the IndexSet.
+   */
+  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:
+  /**
+     @brief Construct a HierarchicApproximation object from a Grid, an IndexSet,
+     and an approximation tolerance
+   */
+  HierarchicApproximation(const Grid &g, const IS &is, double tolerance)
+      : grid_(g), indexSet_(is), tolerance_(tolerance) {}
+
+  /**
+     @brief Search the IndexSet of this HierarchicSearch for the Entity
+     closest to the point global.
+   */
+  Entity findEntity(const Dune::FieldVector<ct, dimw> &global) const {
+    return findEntity<Dune::All_Partition>(global);
+  }
+
+  /**
+     @brief Search the IndexSet of this HierarchicSearch for the Entity
+     closest to the point 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