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/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); diff --git a/src/octTree.cpp b/src/octTree.cpp new file mode 100644 index 0000000..5435b25 --- /dev/null +++ b/src/octTree.cpp @@ -0,0 +1,159 @@ +#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 + + 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]; + } + } + + 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); + + 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++) + { + 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) + { + // 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..6dd846d --- /dev/null +++ b/src/octTree.hpp @@ -0,0 +1,121 @@ +#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; + float lenNorm; + float xMin[3]; + + + template + void walkTreeElements(FunT f, octPtr node, + const OctCoords& icoord, + octCoordType halfNodeLength) + { + OctCoords newCoord; + FCoordinates center, realCenter; + + for (int j = 0; j < 3; j++) + { + center[j] = icoord[j]/(2.*octCoordCenter); + realCenter[j] = xMin[j] + center[j]*lenNorm; + } + + f(realCenter, cells[node].numberLeaves, lenNorm*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++) + { + center[j] = newCoord[j]/(2.*octCoordCenter); + realCenter[j] = xMin[j] + lenNorm*center[j]; + } + + f(realCenter, + 1, lenNorm*halfNodeLength/(2.*octCoordCenter), + false, true); + continue; + } + walkTreeElements(f, cells[node].children[i], newCoord, halfNodeLength/2); + } + + } + + }; + +}; + + +#endif