Commit 7fccb4cf authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

kdtree implementation using entity-seeds

parent ee990734
...@@ -6,27 +6,9 @@ ...@@ -6,27 +6,9 @@
#include <iostream> #include <iostream>
namespace AMDiS { namespace extensions { namespace AMDiS
{
using namespace nanoflann; template <class FeSpace>
// =====================================================================================
typedef WorldVector<double> PointType;
typedef boost::tuple<const MacroElement*, int, unsigned long> DataType;
typedef std::vector<PointType> VectorOfPointsType;
typedef std::vector<DataType> VectorOfDataType;
// =====================================================================================
/** A simple vector-of-vectors adaptor for nanoflann, without duplicating the storage.
* The i'th vector represents a point in the state space.
*
* \tparam dim If set to >0, it specifies a compile-time fixed dimensionality for the points in the data set, allowing more compiler optimizations.
* \tparam value_type The type of the point coordinates (typically, double or float).
* \tparam Distance The distance metric to use: nanoflann::metric_L1, nanoflann::metric_L2, nanoflann::metric_L2_Simple, etc.
* \tparam index_type The type for indices in the KD-tree index (typically, size_t of int)
*/
template < class FeSpace >
struct KDTreeMeshAdaptor struct KDTreeMeshAdaptor
{ {
using MeshView = typename FeSpace::GridView; using MeshView = typename FeSpace::GridView;
...@@ -40,17 +22,17 @@ typedef std::vector<DataType> VectorOfDataType; ...@@ -40,17 +22,17 @@ typedef std::vector<DataType> VectorOfDataType;
using PointType = typename Codim0::Geometry::GlobalCoordinate; using PointType = typename Codim0::Geometry::GlobalCoordinate;
using VectorOfPointsType = std::vector<PointType> ; using VectorOfPointsType = std::vector<PointType> ;
using DataType = std::vector<Dune::FieldVector<double,1> >; using DataType = typename MeshView::Grid::EntitySeed;
using VectorOfDataType = std::vector<std::pair<std::hash<DataType>, DataType>>>; using VectorOfDataType = std::vector<DataType>;
using Distance = nanoflann::metric_L2_Simple; using Distance = nanoflann::metric_L2_Simple;
using Self = KDTreeMeshAdaptor<VectorOfPointsType, VectorOfDataType, ValueType, dimensionworld, Distance, SizeType>; using Self = KDTreeMeshAdaptor<FeSpace>;
using metric_t = typename Distance::template traits<ValueType, Self>::distance_t; using Metric = typename Distance::template traits<ValueType, Self>::distance_t;
using index_t = KDTreeSingleIndexAdaptor<metric_t, Self, dimensionworld, SizeType>; using IndexType = nanoflann::KDTreeSingleIndexAdaptor<Metric, Self, dimensionworld, SizeType>;
public: public:
unique_ptr<index_t> index; //! The kd-tree index for the user to call its methods as usual with any other FLANN index. unique_ptr<IndexType> index; //! The kd-tree index for the user to call its methods as usual with any other FLANN index.
FeSpace const& feSpace; FeSpace const& feSpace;
VectorOfPointsType_ m_points; VectorOfPointsType_ m_points;
...@@ -63,11 +45,13 @@ typedef std::vector<DataType> VectorOfDataType; ...@@ -63,11 +45,13 @@ typedef std::vector<DataType> VectorOfDataType;
KDTreeMeshAdaptor(FeSpace const& feSpace_, KDTreeMeshAdaptor(FeSpace const& feSpace_,
const int leaf_max_size = 10) const int leaf_max_size = 10)
: feSpace(feSpace_) : feSpace(feSpace_)
, localView(feSpace.localView())
, localIndexSet(feSpace.localIndexSet())
{ {
init(); init();
index = make_unique<index_t>(dimensionworld, *this /* adaptor */, index = make_unique<IndexType>(dimensionworld, *this /* adaptor */,
nanoflann::KDTreeSingleIndexAdaptorParams(leaf_max_size, dimensionworld)); nanoflann::KDTreeSingleIndexAdaptorParams(leaf_max_size, dimensionworld));
index->buildIndex(); index->buildIndex();
} }
...@@ -76,21 +60,9 @@ typedef std::vector<DataType> VectorOfDataType; ...@@ -76,21 +60,9 @@ typedef std::vector<DataType> VectorOfDataType;
void init() void init()
{ {
auto const& localView = feSpace.localView();
auto const& indexSet = feSpace.localIndexSet();
for (auto const& element : elements(feSpace.gridView())) { for (auto const& element : elements(feSpace.gridView())) {
localView.bind(element);
localIndexSet.bind(localView);
m_points.push_back(element.geometry().center()); m_points.push_back(element.geometry().center());
m_data.push_back(element.seed());
auto const& localBasis = localView.tree().finiteElement().localBasis();
std::vector<SizeType> indices(localBasis.size());
for (size_t j = 0; j < localBasis.size(); ++j)
indices[i] = localIndexSet.index(j);
m_data.push_back(indices);
} }
} }
...@@ -99,8 +71,8 @@ typedef std::vector<DataType> VectorOfDataType; ...@@ -99,8 +71,8 @@ typedef std::vector<DataType> VectorOfDataType;
m_points.clear(); m_points.clear();
m_data.clear(); m_data.clear();
init(); init();
index.reset(new index_t(dimensionworld, *this /* adaptor */, index.reset(new IndexType(dimensionworld, *this /* adaptor */,
nanoflann::KDTreeSingleIndexAdaptorParams(leaf_max_size, dimensionworld)); nanoflann::KDTreeSingleIndexAdaptorParams(leaf_max_size, dimensionworld));
index->buildIndex(); index->buildIndex();
} }
...@@ -139,20 +111,20 @@ typedef std::vector<DataType> VectorOfDataType; ...@@ -139,20 +111,20 @@ typedef std::vector<DataType> VectorOfDataType;
} }
// Returns the distance between the vector "p1[0:size-1]" and the data point with index "idx_p2" stored in the class: // Returns the distance between the vector "p1[0:size-1]" and the data point with index "idx_p2" stored in the class:
inline value_type kdtree_distance(const double *p1, inline ValueType kdtree_distance(const double *p1,
const SizeType idx_p2, const SizeType idx_p2,
size_t size) const size_t size) const
{ {
double s = 0.0; ValueType s = 0.0;
for (size_t i = 0; i < size; i++) { for (size_t i = 0; i < size; i++) {
const double d = p1[i] - m_points[idx_p2][i]; const ValueType d = p1[i] - m_points[idx_p2][i];
s += d*d; s += d*d;
} }
return s; return s;
} }
// Returns the dim'th component of the idx'th point in the class: // Returns the dim'th component of the idx'th point in the class:
inline value_type kdtree_get_pt(const SizeType idx, size_t dim) const inline ValueType kdtree_get_pt(const SizeType idx, size_t dim) const
{ {
return m_points[idx][dim]; return m_points[idx][dim];
} }
...@@ -174,72 +146,75 @@ typedef std::vector<DataType> VectorOfDataType; ...@@ -174,72 +146,75 @@ typedef std::vector<DataType> VectorOfDataType;
* == evalAtPoint_simple * == evalAtPoint_simple
**/ **/
template <class DOFVector> template <class DOFVector>
typename DOFVector::value_type Value_t<DOFVector> evalAtPoint(DOFVector const& vec, PointType const& x)
evalAtPoint(DOFVector const& vec, PointType const& x, int strategy = 0, int nrOfPoints = -1)
{ {
using T = typename DOFVector::value_type; Value_t<DOFVector> value = 0;
T value = 0;
std::vector<SizeType> indices; auto const& grid = vec.getFeSpace().gridView().grid();
if (getNearestElement(x, indices)) { DataType elementSeed;
// get lambda of point... if (getNearestElement(x, elementSeed)) {
auto element = grid.entity(elementSeed);
return localEval(vec, indices, lambda); auto lambda = element.geometry().local(x);
value = localEval(vec, element, lambda);
} }
return value; return value;
} }
private: private:
/// evaluate DOFVector at barycentric coordinates \p lambda in element that /// evaluate DOFVector at barycentric coordinates \p lambda in element
/// is bind to instance in \ref bind(). template <class FE, class Value, class Element, class LocalCoords>
template <class FE, class ValueType, class Indices, class LocalCoord> Value localEval(DOFVector<FE, Value> const& vec,
ValueType localEval(DOFVector<FE, ValueType> const& vec, Element const& element,
Indices const& indices, LocalCoord const& lambda) LocalCoords const& lambda)
{ {
auto const& localBasis = localView.tree().finiteElement().localBasis(); localView.bind(element);
localIndexSet.bind(localView);
// find a way to evaluate basisFcts at lambda!... auto const& localFiniteElem = localView.tree().finiteElement();
auto const& shape = shapeValues[indices.first]; auto const& localBasis = localFiniteElem.localBasis();
const size_t nBasisFct = localFiniteElem.size();
std::vector<Dune::FieldVector<double,1> > shapeValues;
localBasis.evaluateFunction(lambda, shapeValues); localBasis.evaluateFunction(lambda, shapeValues);
assert( shape.size() == indices.size() ); Value data = 0;
for (size_t j = 0; j < nBasisFct; ++j) {
ValueType data = 0; auto const global_idx = localIndexSet.index(j);
for (size_t j = 0; j < shape.size(); ++j) data += vector[global_idx] * shapeValues[j];
data += vec[indices.second[j]] * shape[j]; }
return data; return data;
} }
template <class Indices>
bool getNearestElement(PointType& x, Indices& indices) bool getNearestElement(PointType const& x, DataType& elementSeed)
{ {
const size_t nnp = 1; const size_t nrOfPoints = 1;
std::vector<DegreeOfFreedom> ret_indexes(nnp); std::vector<DegreeOfFreedom> ret_indexes(nrOfPoints);
std::vector<double> out_dists_sqr(nnp); std::vector<double> out_dists_sqr(nrOfPoints);
nanoflann::KNNResultSet<double, DegreeOfFreedom> resultSet(nnp); nanoflann::KNNResultSet<double, DegreeOfFreedom> resultSet(nrOfPoints);
resultSet.init(&ret_indexes[0], &out_dists_sqr[0]); resultSet.init(&ret_indexes[0], &out_dists_sqr[0]);
index->findNeighbors(resultSet, x.begin(), nanoflann::SearchParams()); index->findNeighbors(resultSet, x.begin(), nanoflann::SearchParams());
indices = m_data[ret_indexes[0]]; elementSeed = m_data[ret_indexes[0]];
return true; return true;
} }
private:
typename FeSpace::LocalView localView;
typename FeSpace::LocalIndexSet localIndexSet;
}; // end of KDTreeMeshAdaptor }; // end of KDTreeMeshAdaptor
template <class FeSpace> template <class FeSpace>
auto provideKDTree(FeSpace const& gridView) auto provideKDTree(FeSpace const& feSpace)
{ {
using GridView = typename FeSpace::GridView; using KD_Tree = KDTreeMeshAdaptor<FeSpace>;
static constexpr int dow = GridView::dimensionworld;
using KD_Tree = KDTreeMeshAdaptor<FeSpace, VectorOfPointsType, VectorOfDataType, double, dow>;
return make_unique<KD_Tree_M>(dow, gridView); return make_unique<KD_Tree>(feSpace);
} }
} // end namespace AMDiS
} }
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment