[cig-commits] r14073 - cs/cigma/trunk/src

[luis at geodynamics.org]

Author: luis
Date: 2009-02-18 08:14:18 -0800 (Wed, 18 Feb 2009)
New Revision: 14073

Added:
   cs/cigma/trunk/src/loc_kdtree.cpp
   cs/cigma/trunk/src/loc_kdtree.h
Log:
Added KdtreeLocator, which uses a different search strategy more suited for both 2d and 3d cells.

Added: cs/cigma/trunk/src/loc_kdtree.cpp
===================================================================
--- cs/cigma/trunk/src/loc_kdtree.cpp	                        (rev 0)
+++ cs/cigma/trunk/src/loc_kdtree.cpp	2009-02-18 16:14:18 UTC (rev 14073)
@@ -0,0 +1,103 @@
+#include "loc_kdtree.h"
+#include "tri_logger.hpp"
+#include <cassert>
+#include <iostream>
+#include <boost/shared_ptr.hpp>
+
+using namespace cigma;
+using boost::shared_ptr;
+
+int KdtreeLocator::nnk = 10;
+
+KdtreeLocator::KdtreeLocator()
+{
+    TRI_LOG_STR("KdtreeLocator::KdtreeLocator()");
+
+    dataPoints = 0;
+    npts = 0;
+    ndim = 0;
+    kdtree = 0;
+
+    queryPoint = 0;
+    searchRadius = 0;
+
+    //nnk = 8;
+    nnIdx = 0;
+    nnDists = 0;
+}
+
+KdtreeLocator::~KdtreeLocator()
+{
+    TRI_LOG_STR("KdtreeLocator::~KdtreeLocator()");
+
+    if (kdtree != 0)
+    {
+        delete kdtree;
+        kdtree = 0;
+    }
+
+    if (dataPoints != 0)
+    {
+        annDeallocPts(dataPoints);
+        dataPoints = 0;
+    }
+
+    if (queryPoint != 0)
+    {
+        annDeallocPt(queryPoint);
+        queryPoint = 0;
+    }
+}
+
+void KdtreeLocator::init(MeshPart& mesh)
+{
+    TRI_LOG_STR("KdtreeLocator::init()");
+
+    npts = mesh.connect->n_cells();
+    ndim = mesh.coords->n_dim();
+
+    assert(nnk > 0);
+    assert(npts > 0);
+    assert((ndim == 2) || (ndim == 3));
+
+    dataPoints = annAllocPts(npts, ndim);
+    queryPoint = annAllocPt(ndim);
+
+    if (mesh.centroid.size() == 0)
+    {
+        mesh.computeCellCentroids();
+    }
+
+    for (int i = 0; i < npts; i++)
+    {
+        ANNpoint pt = dataPoints[i];
+        for (int j = 0; j < ndim; j++)
+        {
+            pt[j] = mesh.centroid[ndim*i + j];
+        }
+    }
+
+    kdtree = new ANNkd_tree(dataPoints, npts, ndim);
+
+    nnIdx.resize(nnk);
+    nnDists.resize(nnk);
+
+    searchRadius = 3 * mesh.getMaxCellDiameter();
+
+    TRI_LOG(npts);
+    TRI_LOG(ndim);
+    TRI_LOG(nnk);
+    TRI_LOG(searchRadius);
+}
+
+void KdtreeLocator::search(double *point)
+{
+    for (int j = 0; j < ndim; j++)
+    {
+        queryPoint[j] = point[j];
+    }
+
+    kdtree->annkSearch(queryPoint, nnk, &nnIdx[0], &nnDists[0], searchRadius);
+}
+
+

Added: cs/cigma/trunk/src/loc_kdtree.h
===================================================================
--- cs/cigma/trunk/src/loc_kdtree.h	                        (rev 0)
+++ cs/cigma/trunk/src/loc_kdtree.h	2009-02-18 16:14:18 UTC (rev 14073)
@@ -0,0 +1,64 @@
+#ifndef CIGMA_KDTREE_LOCATOR_H
+#define CIGMA_KDTREE_LOCATOR_H
+
+#include "Locator.h"
+#include "MeshPart.h"
+#include "ANN/ANN.h"
+#include <valarray>
+
+namespace cigma
+{
+    class KdtreeLocator;
+}
+
+class cigma::KdtreeLocator : public cigma::Locator
+{
+public:
+    KdtreeLocator();
+    ~KdtreeLocator();
+    LocatorType getType() const { return Locator::KdtreeType; }
+
+    void init(MeshPart& mesh);
+    
+    int n_dim() const;
+    void search(double *point);
+
+    int n_idx() const;
+    int idx(int i) const;
+
+public:
+    static int nnk;         /// number of nearest neighbors
+
+private:
+    int npts;               /// number of points
+    int ndim;               /// dimension of single point
+    double searchRadius;    /// search radius
+
+private:
+    // XXX: reusing data in nc_array?
+    ANNpointArray dataPoints;
+    ANNkd_tree *kdtree;
+    ANNpoint queryPoint;                /// query point for search routine
+
+public:
+    std::valarray<ANNidx> nnIdx;        /// near neighbor indices
+    std::valarray<ANNdist> nnDists;     /// near neighbor distances
+};
+
+inline int cigma::KdtreeLocator::n_dim() const
+{
+    return ndim;
+}
+
+inline int cigma::KdtreeLocator::n_idx() const
+{
+    return nnk;
+}
+
+inline int cigma::KdtreeLocator::idx(int i) const
+{
+    //cigma_assert((0 <= i) && (i < nnk));
+    return nnIdx[i];
+}
+
+#endif