[cig-commits] r9132 - in cs/benchmark/cigma/trunk/src: . tests
luis at geodynamics.org
luis at geodynamics.org
Fri Jan 25 07:37:17 PST 2008
Author: luis
Date: 2008-01-25 07:37:17 -0800 (Fri, 25 Jan 2008)
New Revision: 9132
Added:
cs/benchmark/cigma/trunk/src/AnnLocator.cpp
cs/benchmark/cigma/trunk/src/AnnLocator.h
cs/benchmark/cigma/trunk/src/tests/TestAnnLocator.cpp
Log:
Added class for finding candidate cells from a point.
* Uses ANN library to create a kd-tree over list of bounding boxes
Added: cs/benchmark/cigma/trunk/src/AnnLocator.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/AnnLocator.cpp (rev 0)
+++ cs/benchmark/cigma/trunk/src/AnnLocator.cpp 2008-01-25 15:37:17 UTC (rev 9132)
@@ -0,0 +1,90 @@
+#include <cassert>
+#include "AnnLocator.h"
+
+
+// ---------------------------------------------------------------------------
+
+cigma::AnnLocator::AnnLocator()
+{
+ nnk = 8;
+ epsilon = 0;
+
+ npts = 0;
+ dim = 0;
+
+ dataPoints = 0;
+ kdtree = 0;
+
+ nnIdx = 0;
+ nnDists = 0;
+}
+
+cigma::AnnLocator::~AnnLocator()
+{
+ if (kdtree != 0) delete kdtree;
+ if (dataPoints != 0) annDeallocPts(dataPoints);
+ if (nnIdx != 0) delete [] nnIdx;
+ if (nnDists != 0) delete [] nnDists;
+}
+
+
+// ---------------------------------------------------------------------------
+
+void cigma::AnnLocator::initialize(MeshPart *meshPart)
+{
+ npts = meshPart->nel;
+ nsd = meshPart->nsd;
+ dim = 2 * nsd;
+
+ assert(npts > 0);
+ assert(nsd > 0);
+
+ dataPoints = annAllocPts(npts, dim);
+ queryPoint = annAllocPt(dim);
+
+ nnIdx = new ANNidx[nnk];
+ nnDists = new ANNdist[nnk];
+
+ int i,j;
+ double minpt[nsd];
+ double maxpt[nsd];
+
+ for (i = 0; i < npts; i++)
+ {
+ meshPart->set_cell(i);
+ meshPart->cell->bbox(minpt, maxpt);
+
+ ANNpoint pt = dataPoints[i];
+ for (j = 0; j < nsd; j++)
+ {
+ pt[nsd*0 + j] = minpt[j];
+ pt[nsd*1 + j] = maxpt[j];
+ }
+ }
+
+ kdtree = new ANNkd_tree(dataPoints, npts, dim);
+}
+
+
+// ---------------------------------------------------------------------------
+
+void cigma::AnnLocator::search(double *globalPoint, int *cellIndices, int k)
+{
+ int i;
+
+ for (i = 0; i < nsd; i++)
+ {
+ queryPoint[nsd*0 + i] = globalPoint[i];
+ queryPoint[nsd*1 + i] = globalPoint[i];
+ }
+
+ kdtree->annkSearch(queryPoint, k, nnIdx, nnDists, epsilon);
+
+ for (i = 0; i < k; i++)
+ {
+ cellIndices[i] = nnIdx[i];
+ }
+}
+
+
+// ---------------------------------------------------------------------------
Added: cs/benchmark/cigma/trunk/src/AnnLocator.h
===================================================================
--- cs/benchmark/cigma/trunk/src/AnnLocator.h (rev 0)
+++ cs/benchmark/cigma/trunk/src/AnnLocator.h 2008-01-25 15:37:17 UTC (rev 9132)
@@ -0,0 +1,47 @@
+#ifndef __ANN_LOCATOR_H__
+#define __ANN_LOCATOR_H__
+
+#include "ANN/ANN.h"
+#include "MeshPart.h"
+
+
+namespace cigma
+{
+ class AnnLocator;
+}
+
+
+
+class cigma::AnnLocator
+{
+public:
+ AnnLocator();
+ ~AnnLocator();
+
+public:
+ void initialize(MeshPart *meshPart);
+
+public:
+ void search(double *globalPoint, int *cellIndices, int k);
+
+public:
+ int nnk; // number of nearest neighbors
+ int npts; // number of points (bounding boxes)
+ int nsd; // spatial dimensions
+ int dim; // dimension of bounding box (2 * nsd)
+ double epsilon; // tolerance (in bbox space)
+
+public:
+ ANNpointArray dataPoints;
+ ANNkd_tree *kdtree;
+
+public:
+ ANNpoint queryPoint; // query point for find()
+ ANNidxArray nnIdx; // near neighbor indices
+ ANNdistArray nnDists; // near neighbord distances
+
+};
+
+
+
+#endif
Added: cs/benchmark/cigma/trunk/src/tests/TestAnnLocator.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/tests/TestAnnLocator.cpp (rev 0)
+++ cs/benchmark/cigma/trunk/src/tests/TestAnnLocator.cpp 2008-01-25 15:37:17 UTC (rev 9132)
@@ -0,0 +1,172 @@
+#include <iostream>
+#include <ctime>
+#include <cassert>
+
+#include "AnnLocator.h"
+#include "MeshPart.h"
+#include "VtkUgMeshPart.h"
+#include "VtkUgReader.h"
+#include "Misc.h"
+#include "StringUtils.h"
+
+#include "Numeric.h"
+#include "Cell.h"
+#include "Tet.h"
+#include "Hex.h"
+
+
+using namespace std;
+using namespace cigma;
+
+int main(int argc, char *argv[])
+{
+ int npts = 100;
+
+ string filename = "strikeslip_tet4_1000m_t0.vtk";
+ string inputstr;
+
+ if (argc > 1)
+ {
+ filename = argv[1];
+ }
+ if (argc > 2)
+ {
+ inputstr = argv[2];
+ string_to_int(inputstr, npts);
+ assert(npts > 0);
+ }
+
+
+ int i,j;
+
+ int nno, nsd;
+ double *coords;
+ int nel, ndofs;
+ int *connect;
+
+
+ VtkUgReader reader;
+ reader.open(filename);
+ reader.get_coordinates(&coords, &nno, &nsd);
+ reader.get_connectivity(&connect, &nel, &ndofs);
+
+
+ MeshPart *meshPart;
+ meshPart = new VtkUgMeshPart();
+ meshPart->set_coordinates(coords, nno, nsd);
+ meshPart->set_connectivity(connect, nel, ndofs);
+ meshPart->cell = 0;
+ switch (ndofs)
+ {
+ case 4:
+ meshPart->cell = new Tet();
+ break;
+ case 8:
+ meshPart->cell = new Hex();
+ break;
+ }
+ assert(meshPart->cell != 0);
+
+
+
+ AnnLocator *locator = new AnnLocator();
+ locator->initialize(meshPart);
+
+
+
+
+ double minpt[nsd], maxpt[nsd];
+ meshPart->get_bbox(minpt, maxpt);
+
+
+
+ double *points = new double[npts * nsd];
+ for (i = 0; i < npts; i++)
+ {
+ double *pt = &points[nsd*i];
+ for (j = 0; j < nsd; j++)
+ {
+ pt[j] = pick_from_interval(minpt[j],maxpt[j]);
+ }
+ }
+
+
+
+ int k = locator->nnk;
+ int *candidateCells= new int[k];
+ int *parentCells = new int[npts];
+
+ for (i = 0; i < npts; i++)
+ {
+ parentCells[i] = -1;
+ }
+
+
+ time_t t0, t1;
+
+ cout << "Starting...";
+ time(&t0);
+ for (i = 0; i < npts; i++)
+ {
+ double *pt = &points[nsd*i];
+ locator->search(pt, candidateCells, k);
+
+ bool found = false;
+ for (j = 0; j < k; j++)
+ {
+ int e = candidateCells[j];
+ meshPart->set_cell(e);
+
+ double uvw[3];
+ meshPart->cell->xyz2uvw(pt, uvw);
+ found = meshPart->cell->interior(uvw[0], uvw[1], uvw[2]);
+ if (found)
+ {
+ parentCells[i] = e;
+ break;
+ }
+ }
+ assert(found);
+ }
+ time(&t1);
+ cout << "done!" << endl;
+
+
+
+ double total_time = (t1 - t0);
+ cout << "time: " << total_time << endl;
+
+ cout << "bbox: ";
+ for (j = 0; j < nsd; j++)
+ {
+ cout << minpt[j] << " ";
+ }
+ for (j = 0; j < nsd; j++)
+ {
+ cout << maxpt[j] << " ";
+ }
+ cout << endl;
+
+
+ cout << "points: " << npts << endl;
+ for (i = 0; i < npts; i++)
+ {
+ double *pt = &points[nsd*i];
+ cout << i << " (";
+ for (j = 0; j < nsd; j++)
+ {
+ cout << pt[j];
+ if (j != nsd-1)
+ cout << " ";
+ }
+ cout << ") -> " << parentCells[i] << endl;
+ }
+
+
+ delete [] points;
+ delete [] parentCells;
+ delete [] candidateCells;
+
+
+ return 0;
+}
More information about the cig-commits
mailing list