From 6422e08208f4989a915e6d10cacef93ff06e78d1 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Sun, 1 Nov 2009 11:20:50 -0600 Subject: [PATCH 1/4] New octtree tool --- lib/Makefile | 4 +- src/octTree.cpp | 132 ++++++++++++++++++++++++++++++++++++++++++++++++ src/octTree.hpp | 113 +++++++++++++++++++++++++++++++++++++++++ 3 files changed, 247 insertions(+), 2 deletions(-) create mode 100644 src/octTree.cpp create mode 100644 src/octTree.hpp diff --git a/lib/Makefile b/lib/Makefile index 7f39d93..0057b7c 100644 --- a/lib/Makefile +++ b/lib/Makefile @@ -1,5 +1,5 @@ SHLIBS= libCosmoTool.so -SOURCES= loadRamses.cpp yorick.cpp miniargs.cpp fortran.cpp interpolate.cpp load_data.cpp powerSpectrum.cpp +SOURCES= loadRamses.cpp yorick.cpp miniargs.cpp fortran.cpp interpolate.cpp load_data.cpp powerSpectrum.cpp octTree.cpp LIBS= -lnetcdf_c++ -lnetcdf -lgsl -lgslcblas -lm include config.mk @@ -8,7 +8,7 @@ VPATH=../src all: $(SHLIBS) -libCosmoTool.so: loadRamses.o yorick.o miniargs.o fortran.o interpolate.o load_data.o powerSpectrum.o +libCosmoTool.so: loadRamses.o yorick.o miniargs.o fortran.o interpolate.o load_data.o powerSpectrum.o octTree.o depend: $(SOURCES) @echo "[DEPENDS] $^" diff --git a/src/octTree.cpp b/src/octTree.cpp new file mode 100644 index 0000000..30e0175 --- /dev/null +++ b/src/octTree.cpp @@ -0,0 +1,132 @@ +#include +#include +#include +#include "config.hpp" +#include "octTree.hpp" + +using namespace std; +using namespace CosmoTool; + +//#define VERBOSE + +OctTree::OctTree(const FCoordinates *particles, octPtr numParticles, + uint32_t maxMeanTreeDepth, uint32_t maxAbsoluteDepth, + uint32_t threshold) +{ + cout << "MeanTree=" << maxMeanTreeDepth << endl; + numCells = pow(8, maxMeanTreeDepth); + assert(numCells < invalidOctCell); + //#ifdef VERBOSE + cerr << "Allocating " << numCells << " octtree cells" << endl; + //#endif + + cells = new OctCell[numCells]; + Lbox = (float)(octCoordTypeNorm+1); + + cells[0].numberLeaves = 0; + for (int i = 0; i < 8; i++) + cells[0].children[i] = emptyOctCell; + + lastNode = 1; + this->particles = particles; + this->numParticles = numParticles; + buildTree(maxAbsoluteDepth); + //#ifdef VERBOSE + cerr << "Used " << lastNode << " cells" << endl; + //#endif +} + +OctTree::~OctTree() +{ + delete cells; +} + +void OctTree::buildTree(uint32_t maxAbsoluteDepth) +{ + for (octPtr i = 0; i < numParticles; i++) + { + OctCoords rootCenter = { octCoordCenter, octCoordCenter, octCoordCenter }; + insertParticle(0, // root node + rootCenter, + octCoordCenter, + i, + maxAbsoluteDepth); + } +} + + +void OctTree::insertParticle(octPtr node, + const OctCoords& icoord, + octCoordType halfNodeLength, + octPtr particleId, + uint32_t maxAbsoluteDepth) +{ + +#ifdef VERBOSE + cout << "Entering " << node << " (" << icoord[0] << "," << icoord[1] << "," << icoord[2] << ")" << endl; +#endif + int octPos = 0; + int ipos[3] = { 0,0,0}; + octPtr newNode; + OctCoords newCoord; + + cells[node].numberLeaves++; + if (maxAbsoluteDepth == 0) + { + // All children must be invalid. + for (int i = 0 ; i < 8; i++) + cells[node].children[i] = invalidOctCell; + + return; + } + + for (int j = 0; j < 3; j++) + if ((octPtr)(particles[particleId][j]*Lbox) > icoord[j]) + { + octPos |= (1 << j); + ipos[j] = 1; + } + + if (cells[node].children[octPos] == emptyOctCell) + { + // Put the particle there. + cells[node].children[octPos] = particleId | octParticleMarker; + return; + } + + // If it is a node, explores it. + if (!(cells[node].children[octPos] & octParticleMarker)) + { + assert(halfNodeLength >= 2); + // Compute coordinates + for (int j = 0; j < 3; j++) + newCoord[j] = icoord[j]+(2*ipos[j]-1)*halfNodeLength/2; + insertParticle(cells[node].children[octPos], newCoord, halfNodeLength/2, + particleId, maxAbsoluteDepth-1); + return; + } + + // We have a particle there. + // Make a new node and insert the old particle into this node. + // Insert the new particle into the node also + // Finally put the node in place + + newNode = lastNode++; + assert(lastNode != numCells); + + for (int j = 0; j < 8; j++) + cells[newNode].children[j] = emptyOctCell; + cells[newNode].numberLeaves = 0; + + // Compute coordinates + for (int j = 0; j < 3; j++) + newCoord[j] = icoord[j]+(2*ipos[j]-1)*halfNodeLength/2; + + octPtr oldPartId = cells[node].children[octPos] & octParticleMask; + + insertParticle(newNode, newCoord, halfNodeLength/2, + oldPartId, maxAbsoluteDepth-1); + insertParticle(newNode, newCoord, halfNodeLength/2, + particleId, maxAbsoluteDepth-1); + cells[node].children[octPos] = newNode; +} diff --git a/src/octTree.hpp b/src/octTree.hpp new file mode 100644 index 0000000..3aee37c --- /dev/null +++ b/src/octTree.hpp @@ -0,0 +1,113 @@ +#ifndef __COSMOTOOL_AMR_HPP +#define __COSMOTOOL_AMR_HPP + +#include "config.hpp" + +namespace CosmoTool +{ + + typedef uint32_t octPtr; + typedef uint16_t octCoordType; + static const uint16_t octCoordTypeNorm = 0xffff; + static const uint16_t octCoordCenter = 0x8000; + + // This is also the root cell, but this one + // is never referenced (really ??). + static const octPtr invalidOctCell = 0xffffffff; + static const octPtr emptyOctCell = 0; + static const octPtr octParticleMarker = 0x80000000; + static const octPtr octParticleMask = 0x7fffffff; + + typedef octCoordType OctCoords[3]; + + struct OctCell + { + octPtr numberLeaves; + octPtr children[8]; + }; + + class OctTree + { + public: + //Coordinates of particles must be in the [0:1] range + OctTree(const FCoordinates *particles, octPtr numParticles, + uint32_t maxTreeDepth, uint32_t maxAbsoluteDepth, + uint32_t threshold = 1); + ~OctTree(); + + void buildTree(uint32_t maxAbsoluteDepth); + void insertParticle(octPtr node, + const OctCoords& icoord, + octCoordType halfNodeLength, + octPtr particleId, + uint32_t maxAbsoluteDepth); + + + octPtr getNumberLeaves() const { + return cells[0].numberLeaves; + } + + template + void walkTree(FunT f) + { + OctCoords rootCenter = { octCoordCenter, octCoordCenter, octCoordCenter }; + + walkTreeElements(f, 0, rootCenter, octCoordCenter); + } + + protected: + const FCoordinates *particles; + octPtr numParticles; + OctCell *cells; + float Lbox; + octPtr lastNode; + octPtr numCells; + + + template + void walkTreeElements(FunT f, octPtr node, + const OctCoords& icoord, + octCoordType halfNodeLength) + { + OctCoords newCoord; + FCoordinates realCenter; + + for (int j = 0; j < 3; j++) + realCenter[j] = icoord[j]/(2.*octCoordCenter); + + f(realCenter, cells[node].numberLeaves, halfNodeLength/(float)octCoordCenter, + cells[node].children[0] == invalidOctCell, // True if this is a meta-node + false); + + for (int i = 0; i < 8; i++) + { + octPtr newNode = cells[node].children[i]; + int ipos[3] = { (i&1), (i&2)>>1, (i&4)>>2 }; + + if (newNode == emptyOctCell || newNode == invalidOctCell) + continue; + + for (int j = 0; j < 3; j++) + newCoord[j] = icoord[j]+(2*ipos[j]-1)*halfNodeLength/2; + + if (newNode & octParticleMarker) + { + for (int j = 0; j < 3; j++) + realCenter[j] = newCoord[j]/(2.*octCoordCenter); + + f(realCenter, + 1, halfNodeLength/(2.*octCoordCenter), + false, true); + continue; + } + walkTreeElements(f, cells[node].children[i], newCoord, halfNodeLength/2); + } + + } + + }; + +}; + + +#endif From bab3c3b1114669c86b0328a621fc78b7dbf9ad4c Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Sun, 1 Nov 2009 11:21:14 -0600 Subject: [PATCH 2/4] Small compilation fixes --- src/miniargs.cpp | 1 + src/mykdtree.hpp | 6 +++--- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/miniargs.cpp b/src/miniargs.cpp index 95a832a..4ad71c6 100644 --- a/src/miniargs.cpp +++ b/src/miniargs.cpp @@ -1,4 +1,5 @@ #include +#include #include #include "miniargs.hpp" #include diff --git a/src/mykdtree.hpp b/src/mykdtree.hpp index 9d83d5b..10c96b4 100644 --- a/src/mykdtree.hpp +++ b/src/mykdtree.hpp @@ -90,9 +90,9 @@ namespace CosmoTool { Cell *getNearestNeighbour(const coords& x); - void getNearestNeighbours(const coords& x, uint32_t N, + void getNearestNeighbours(const coords& x, uint32_t NumCells, Cell **cells); - void getNearestNeighbours(const coords& x, uint32_t N, + void getNearestNeighbours(const coords& x, uint32_t NumCells, Cell **cells, CoordType *distances); @@ -114,7 +114,7 @@ namespace CosmoTool { Cell **sortingHelper; Node *buildTree(Cell **cell0, - uint32_t N, + uint32_t NumCells, uint32_t depth, coords minBound, coords maxBound); From bdb98252138733fc693129d7f0f3dbc4e2c7f40d Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Sun, 1 Nov 2009 18:56:10 -0600 Subject: [PATCH 3/4] Adapt for a more general octTree --- src/octTree.cpp | 13 +++++++++++++ src/octTree.hpp | 14 +++++++++++--- 2 files changed, 24 insertions(+), 3 deletions(-) diff --git a/src/octTree.cpp b/src/octTree.cpp index 30e0175..91648cf 100644 --- a/src/octTree.cpp +++ b/src/octTree.cpp @@ -19,6 +19,19 @@ OctTree::OctTree(const FCoordinates *particles, octPtr numParticles, //#ifdef VERBOSE cerr << "Allocating " << numCells << " octtree cells" << endl; //#endif + + for (int j = 0; j < 3; j++) + xMin[j] = particles[0][j]; + + for (octPtr i = 1; i < numParticles; i++) + { + for (int j = 0; j < 3; j++) + { + if (particles[i][j] < xMin[j]) + xMin[j] = particles[i][j]; + } + } + cells = new OctCell[numCells]; Lbox = (float)(octCoordTypeNorm+1); diff --git a/src/octTree.hpp b/src/octTree.hpp index 3aee37c..2402a5e 100644 --- a/src/octTree.hpp +++ b/src/octTree.hpp @@ -62,6 +62,8 @@ namespace CosmoTool float Lbox; octPtr lastNode; octPtr numCells; + float lenNorm; + float xMin[3]; template @@ -70,10 +72,13 @@ namespace CosmoTool octCoordType halfNodeLength) { OctCoords newCoord; - FCoordinates realCenter; + FCoordinates center, realCenter; for (int j = 0; j < 3; j++) - realCenter[j] = icoord[j]/(2.*octCoordCenter); + { + center[j] = icoord[j]/(2.*octCoordCenter); + realCenter[j] = center[j]*lenNorm+xMin[j]; + } f(realCenter, cells[node].numberLeaves, halfNodeLength/(float)octCoordCenter, cells[node].children[0] == invalidOctCell, // True if this is a meta-node @@ -93,7 +98,10 @@ namespace CosmoTool if (newNode & octParticleMarker) { for (int j = 0; j < 3; j++) - realCenter[j] = newCoord[j]/(2.*octCoordCenter); + { + center[j] = newCoord[j]/(2.*octCoordCenter); + realCenter[j] = xMin[j] + lenNorm*center[j]; + } f(realCenter, 1, halfNodeLength/(2.*octCoordCenter), From a4ffa66b8da08e1bc1164a502ad3e434933fd234 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Mon, 2 Nov 2009 07:27:00 -0600 Subject: [PATCH 4/4] Fixed representation and walking of arbitrary sized OctTree --- src/octTree.cpp | 24 +++++++++++++++++++----- src/octTree.hpp | 6 +++--- 2 files changed, 22 insertions(+), 8 deletions(-) diff --git a/src/octTree.cpp b/src/octTree.cpp index 91648cf..5435b25 100644 --- a/src/octTree.cpp +++ b/src/octTree.cpp @@ -32,6 +32,17 @@ OctTree::OctTree(const FCoordinates *particles, octPtr numParticles, } } + lenNorm = 0; + for (octPtr i = 0; i < numParticles; i++) + { + for (int j = 0; j < 3; j++) + { + float delta = particles[i][j]-xMin[j]; + if (delta > lenNorm) + lenNorm = delta; + } + } + cout << xMin[0] << " " << xMin[1] << " " << xMin[2] << " lNorm=" << lenNorm << endl; cells = new OctCell[numCells]; Lbox = (float)(octCoordTypeNorm+1); @@ -94,11 +105,14 @@ void OctTree::insertParticle(octPtr node, } for (int j = 0; j < 3; j++) - if ((octPtr)(particles[particleId][j]*Lbox) > icoord[j]) - { - octPos |= (1 << j); - ipos[j] = 1; - } + { + float treePos = (particles[particleId][j]-xMin[j])*Lbox/lenNorm; + if ((octPtr)(treePos) > icoord[j]) + { + octPos |= (1 << j); + ipos[j] = 1; + } + } if (cells[node].children[octPos] == emptyOctCell) { diff --git a/src/octTree.hpp b/src/octTree.hpp index 2402a5e..6dd846d 100644 --- a/src/octTree.hpp +++ b/src/octTree.hpp @@ -77,10 +77,10 @@ namespace CosmoTool for (int j = 0; j < 3; j++) { center[j] = icoord[j]/(2.*octCoordCenter); - realCenter[j] = center[j]*lenNorm+xMin[j]; + realCenter[j] = xMin[j] + center[j]*lenNorm; } - f(realCenter, cells[node].numberLeaves, halfNodeLength/(float)octCoordCenter, + f(realCenter, cells[node].numberLeaves, lenNorm*halfNodeLength/(float)octCoordCenter, cells[node].children[0] == invalidOctCell, // True if this is a meta-node false); @@ -104,7 +104,7 @@ namespace CosmoTool } f(realCenter, - 1, halfNodeLength/(2.*octCoordCenter), + 1, lenNorm*halfNodeLength/(2.*octCoordCenter), false, true); continue; }