diff --git a/sample/simple3DFilter.cpp b/sample/simple3DFilter.cpp index b3a74d2..345834c 100644 --- a/sample/simple3DFilter.cpp +++ b/sample/simple3DFilter.cpp @@ -7,6 +7,8 @@ #include #include "hdf5_array.hpp" #include +#include +#include using namespace std; using namespace CosmoTool; @@ -17,8 +19,13 @@ struct VCoord{ float v[3]; }; -template -ComputePrecision getVelocity(const VCoord& v) +using boost::format; +using boost::str; +typedef boost::multi_array array_type; +typedef boost::multi_array array3_type; +typedef boost::multi_array array4_type; + +ComputePrecision getVelocity(const VCoord& v, int i) { return v.v[i]; } @@ -32,11 +39,51 @@ typedef SPHSmooth MySmooth; typedef MySmooth::SPHTree MyTree; typedef MyTree::Cell MyCell; +template +void computeInterpolatedField(MyTree *tree1, double boxsize, int Nres, double cx, double cy, double cz, + array3_type& bins, array3_type& arr, FuncT func, double rLimit2) +{ +#pragma omp parallel + { + MySmooth smooth1(tree1, N_SPH); + +#pragma omp for schedule(dynamic) + for (int rz = 0; rz < Nres; rz++) + { + double pz = (rz)*boxsize/Nres-cz; + + cout << format("[%d] %d / %d") % omp_get_thread_num() % rz % Nres << endl; + for (int ry = 0; ry < Nres; ry++) + { + double py = (ry)*boxsize/Nres-cy; + for (int rx = 0; rx < Nres; rx++) + { + double px = (rx)*boxsize/Nres-cx; + + MyTree::coords c = { px, py, pz }; + + double r2 = c[0]*c[0]+c[1]*c[1]+c[2]*c[2]; + if (r2 > rLimit2) + { + arr[rx][ry][rz] = 0; + continue; + } + + uint32_t numInCell = bins[rx][ry][rz]; + if (numInCell > N_SPH) + smooth1.fetchNeighbours(c, numInCell); + else + smooth1.fetchNeighbours(c); + + arr[rx][ry][rz] = smooth1.computeSmoothedValue(c, func); + } + } + } + } +} + int main(int argc, char **argv) { - typedef boost::multi_array array_type; - typedef boost::multi_array array3_type; - typedef boost::multi_array array4_type; char *fname1, *fname2; double rLimit, boxsize, rLimit2, cx, cy, cz; @@ -105,8 +152,6 @@ int main(int argc, char **argv) cout << "Creating smoothing filter..." << endl; - array3_type out_den_1(boost::extents[Nres][Nres][Nres]); - array4_type out_v3d_1(boost::extents[Nres][Nres][Nres][3]); // array3_type out_rad_1(boost::extents[Nres][Nres][Nres]); cout << "Weighing..." << endl; @@ -149,55 +194,18 @@ int main(int argc, char **argv) } cout << "Interpolating..." << endl; -#pragma omp parallel - { - MySmooth smooth1(&tree1, N_SPH); - -#pragma omp for schedule(dynamic) - for (int rz = 0; rz < Nres; rz++) - { - double pz = (rz)*boxsize/Nres-cz; - cout << rz << " / " << Nres << endl; - for (int ry = 0; ry < Nres; ry++) - { - double py = (ry)*boxsize/Nres-cy; - for (int rx = 0; rx < Nres; rx++) - { - double px = (rx)*boxsize/Nres-cx; - - MyTree::coords c = { px, py, pz }; - - double r2 = c[0]*c[0]+c[1]*c[1]+c[2]*c[2]; - if (r2 > rLimit2) - { - out_v3d_1[rx][ry][rz][0] = 0; - out_v3d_1[rx][ry][rz][1] = 0; - out_v3d_1[rx][ry][rz][2] = 0; - out_den_1[rx][ry][rz] = 0; - //out_rad_1[rx][ry][rz] = -1; - continue; - } - - uint32_t numInCell = bins[rx][ry][rz]; - if (numInCell > N_SPH) - smooth1.fetchNeighbours(c, numInCell); - else - smooth1.fetchNeighbours(c); - - //out_rad_1[rx][ry][rz] = smooth1.getSmoothingLen(); - out_v3d_1[rx][ry][rz][0] = smooth1.computeSmoothedValue(c, getVelocity<0>); - out_v3d_1[rx][ry][rz][1] = smooth1.computeSmoothedValue(c, getVelocity<1>); - out_v3d_1[rx][ry][rz][2] = smooth1.computeSmoothedValue(c, getVelocity<2>); - out_den_1[rx][ry][rz] = smooth1.computeSmoothedValue(c, getUnity); - } - } - } + array3_type interpolated(boost::extents[Nres][Nres][Nres]); + + computeInterpolatedField(&tree1, boxsize, Nres, cx, cy, cz, + bins, interpolated, getUnity, rLimit2); + hdf5_write_array(out_f, "density", interpolated); + //out_f.flush(); + for (int i = 0; i < 3; i++) { + computeInterpolatedField(&tree1, boxsize, Nres, cx, cy, cz, + bins, interpolated, boost::bind(getVelocity, _1, i), rLimit2); + hdf5_write_array(out_f, str(format("p%d") % i), interpolated); } - - //hdf5_write_array(out_f, "radii", out_rad_1); - hdf5_write_array(out_f, "velocity", out_v3d_1); - hdf5_write_array(out_f, "density", out_den_1); - + return 0; }; diff --git a/src/sphSmooth.hpp b/src/sphSmooth.hpp index d116dd4..bcd98a5 100644 --- a/src/sphSmooth.hpp +++ b/src/sphSmooth.hpp @@ -82,10 +82,14 @@ namespace CosmoTool return internal.currentCenter; } + template ComputePrecision computeSmoothedValue(const typename SPHTree::coords& c, - computeParticleValue fun, SPHState *state = 0); + FuncT fun, SPHState *state = 0); + + template ComputePrecision computeInterpolatedValue(const typename SPHTree::coords& c, - computeParticleValue fun, SPHState *state = 0); + FuncT fun, SPHState *state = 0); + ComputePrecision getMaxDistance(const typename SPHTree::coords& c, SPHNode *node) const; @@ -101,7 +105,8 @@ namespace CosmoTool } // END - void runForEachNeighbour(runParticleValue fun, SPHState *state = 0); + template + void runForEachNeighbour(FuncT fun, SPHState *state = 0); void addGridSite(const typename SPHTree::coords& c); bool hasNeighbours() const; @@ -127,12 +132,15 @@ namespace CosmoTool uint32_t maxNgb; SPHTree *tree; + template ComputePrecision computeWValue(const typename SPHTree::coords & c, SPHCell& cell, CoordType d, - computeParticleValue fun, SPHState *state); + FuncT fun, SPHState *state); + + template void runUnrollNode(SPHNode *node, - runParticleValue fun); + FuncT fun); }; template diff --git a/src/sphSmooth.tcc b/src/sphSmooth.tcc index 1ae5486..99094cd 100644 --- a/src/sphSmooth.tcc +++ b/src/sphSmooth.tcc @@ -21,10 +21,11 @@ SPHSmooth::~SPHSmooth() } template +template ComputePrecision SPHSmooth::computeWValue(const typename SPHTree::coords& c, SPHCell& cell, CoordType d, - computeParticleValue fun, SPHState *state) + FuncT fun, SPHState *state) { CoordType weight; @@ -117,9 +118,10 @@ SPHSmooth::fetchNeighboursOnVolume(const typename SPHTree::coords } template +template ComputePrecision SPHSmooth::computeSmoothedValue(const typename SPHTree::coords& c, - computeParticleValue fun, SPHState *state) + FuncT fun, SPHState *state) { if (state == 0) state = &internal; @@ -144,8 +146,9 @@ ComputePrecision interpolateOne(const ValType& t) // WARNING ! Cell's weight must be 1 !!! template +template ComputePrecision SPHSmooth::computeInterpolatedValue(const typename SPHTree::coords& c, - computeParticleValue fun, SPHState *state) + FuncT fun, SPHState *state) { if (state == 0) state = &internal; @@ -164,7 +167,8 @@ ComputePrecision SPHSmooth::computeInterpolatedValue(const typena } template -void SPHSmooth::runForEachNeighbour(runParticleValue fun, SPHState *state) +template +void SPHSmooth::runForEachNeighbour(FuncT fun, SPHState *state) { if (state == 0) state = &internal;