Fixed conflict

CMakeLists.txt

@@ -31,6 +31,8 @@ find_library(NETCDF_LIBRARY netcdf)
 find_library(GSL_LIBRARY gsl)
 find_library(GSLCBLAS_LIBRARY gslcblas)
 
+find_package(Boost 1.53)
+
 if (ENABLE_SHARP)
   SET(SHARP_SOURCE ${CMAKE_SOURCE_DIR}/external/sharp)
   SET(DEP_BUILD ${CMAKE_SOURCE_DIR}/external/sharp/auto)

sample/CMakeLists.txt

@@ -1,5 +1,7 @@
 SET(tolink ${GSL_LIBRARIES} CosmoTool ${CosmoTool_LIBS})
-include_directories(${CMAKE_SOURCE_DIR}/src ${FFTW3_INCLUDE_DIRS} ${EIGEN3_INCLUDE_DIRS} ${NETCDF_INCLUDE_PATH} ${GSL_INCLUDE_PATH})
+include_directories(${CMAKE_SOURCE_DIR}/src)
+include_directories(${FFTW3_INCLUDE_DIRS}
+${EIGEN3_INCLUDE_DIRS} ${NETCDF_INCLUDE_PATH} ${GSL_INCLUDE_PATH})
 
 IF(SHARP_INCLUDE_PATH)
   include_directories(BEFORE ${SHARP_INCLUDE_PATH})
@@ -11,9 +13,6 @@ target_link_libraries(testBQueue ${tolink})
 add_executable(testInterpolate testInterpolate.cpp)
 target_link_libraries(testInterpolate ${tolink})
 
-#add_executable(testSmooth testSmooth.cpp)
-#target_link_libraries(testSmooth ${tolink})
-
 add_executable(testkd testkd.cpp)
 target_link_libraries(testkd ${tolink})
 
@@ -34,7 +33,7 @@ target_link_libraries(testPool ${tolink})
 
 if (HDF5_FOUND)
 	  add_executable(testReadFlash testReadFlash.cpp)
-	  target_link_libraries(testReadFlash ${tolink})
+    target_link_libraries(testReadFlash ${tolink})
 endif (HDF5_FOUND)
 
 
@@ -69,3 +68,15 @@ endif (ENABLE_SHARP AND SHARP_LIBRARY AND SHARP_INCLUDE_PATH AND EIGEN3_FOUND)
 add_executable(test_cosmopower test_cosmopower.cpp)
 target_link_libraries(test_cosmopower ${tolink})
 
+if (Boost_FOUND)
+  include_directories(${Boost_INCLUDE_DIRS})
+
+  add_executable(testSmooth testSmooth.cpp)
+  target_link_libraries(testSmooth ${tolink})
+
+  add_executable(simple3DFilter simple3DFilter.cpp)
+  target_link_libraries(simple3DFilter ${tolink})
+endif (Boost_FOUND)
+
+add_executable(gadgetToDensity gadgetToDensity.cpp)
+target_link_libraries(gadgetToDensity ${tolink})

sample/gadgetToDensity.cpp

@@ -0,0 +1,62 @@
+#include <cmath>
+#include <iostream>
+#include <cstdlib>
+#include "cic.hpp"
+#include "loadGadget.hpp"
+#include "miniargs.hpp"
+#include "yorick.hpp"
+
+using namespace std;
+using namespace CosmoTool;
+
+int main(int argc, char **argv)
+{
+  uint32_t res;
+  char *fname;
+  int id;
+
+  MiniArgDesc desc[] = {
+    { "SNAPSHOT", &fname, MINIARG_STRING },
+    { "ID", &id, MINIARG_INT },
+    { "RESOLUTION", &res, MINIARG_INT },
+    { 0, 0, MINIARG_NULL }
+  };
+
+  if (!parseMiniArgs(argc, argv, desc))
+    return 1;
+
+  SimuData *p = loadGadgetMulti(fname, 0, 0);
+  double L0 = p->BoxSize/1000;
+  CICFilter filter(res, L0);
+
+  delete p;
+
+  try {
+  for (int cpuid=0;;cpuid++) {
+    p = loadGadgetMulti(fname, cpuid, NEED_POSITION);
+    for (uint32_t i = 0; i < p->NumPart; i++)
+    {
+      CICParticles a;
+      
+      a.mass = 1.0;
+      a.coords[0] = p->Pos[0][i]/1000;
+      a.coords[1] = p->Pos[1][i]/1000;
+      a.coords[2] = p->Pos[2][i]/1000;
+      filter.putParticles(&a, 1);
+    }
+  
+    delete p;
+  }
+  } catch (const NoSuchFileException& e) {}
+
+  CICType *denField;
+  uint32_t Ntot;
+  filter.getDensityField(denField, Ntot);
+  
+  cout << "L0=" << L0 << endl;
+  cout << "Saving density field" << endl;
+  uint32_t dimList[] = { res, res, res};
+  saveArray("densityField.nc", denField, dimList, 3);
+
+  return 0;
+}

sample/simple3DFilter.cpp

@@ -0,0 +1,193 @@
+#include <cassert>
+#include "yorick.hpp"
+#include "sphSmooth.hpp"
+#include "mykdtree.hpp"
+#include "miniargs.hpp"
+
+using namespace std;
+using namespace CosmoTool;
+
+#define N_SPH 16
+
+struct VCoord{
+  float v[3];
+};
+
+template<int i>
+ComputePrecision getVelocity(const VCoord& v)
+{
+  return v.v[i];
+}
+
+ComputePrecision getUnity(const VCoord& v)
+{
+  return 1.0;
+}
+
+typedef SPHSmooth<VCoord> MySmooth;
+typedef MySmooth::SPHTree MyTree;
+typedef MyTree::Cell MyCell;
+
+int main(int argc, char **argv)
+{
+  char *fname1, *fname2;
+  double rLimit, boxsize, rLimit2, cx, cy, cz;
+  int Nres;
+
+  MiniArgDesc args[] = {
+    { "INPUT DATA1", &fname1, MINIARG_STRING },
+    { "RADIUS LIMIT", &rLimit, MINIARG_DOUBLE },
+    { "BOXSIZE", &boxsize, MINIARG_DOUBLE },
+    { "RESOLUTION", &Nres, MINIARG_INT },
+    { "CX", &cx, MINIARG_DOUBLE },
+    { "CY", &cy, MINIARG_DOUBLE },
+    { "CZ", &cz, MINIARG_DOUBLE },
+    { 0, 0, MINIARG_NULL }
+  };
+
+  if (!parseMiniArgs(argc, argv, args))
+    return 1;
+
+  float *v1_data;
+  uint32_t *dimList;
+  uint32_t rank;
+  uint32_t N1_points, N2_points;
+  int *bins = new int[Nres*Nres*Nres];
+
+  rLimit2 = rLimit*rLimit;
+
+  loadArray(fname1, v1_data, dimList, rank);
+  assert(rank == 2);
+  assert(dimList[1] == 6);
+
+  N1_points = dimList[0];
+  delete[] dimList;
+  
+  cout << "Got " << N1_points << " in the first file." << endl;
+
+  MyCell *allCells_1 = new MyCell[N1_points];
+  
+  for (long i = 0; i < Nres*Nres*Nres; i++)
+    bins[i] = 0;
+
+  cout << "Shuffling data in cells..." << endl;
+  for (int i = 0 ; i < N1_points; i++)
+    {
+      for (int j = 0; j < 3; j++)
+	allCells_1[i].coord[j] = v1_data[i*6 + j];
+      for (int k = 0; k < 3; k++)
+	allCells_1[i].val.pValue.v[k] = v1_data[i*6 + 3 + k];
+      allCells_1[i].active = true;
+      allCells_1[i].val.weight = 0.0;
+
+      long rx = floor((allCells_1[i].coord[0]+cx)*Nres/boxsize+0.5);
+      long ry = floor((allCells_1[i].coord[1]+cy)*Nres/boxsize+0.5);
+      long rz = floor((allCells_1[i].coord[2]+cz)*Nres/boxsize+0.5);
+      
+      if (rx < 0 || rx >= Nres || ry < 0 || ry >= Nres || rz < 0 || rz >= Nres)
+	  continue;
+	  
+      bins[rx + ry*Nres + rz*Nres*Nres]++;
+    }
+  delete[] v1_data;
+   uint32_t dims[3] = { Nres, Nres, Nres } ;
+   saveArray("num_in_cell.nc", bins, dims, 3);
+
+  cout << "Building trees..." << endl;
+  MyTree tree1(allCells_1, N1_points);
+
+  cout << "Creating smoothing filter..." << endl;
+  MySmooth smooth1(&tree1, N_SPH);
+
+  uint32_t outDimList[3] = { Nres, Nres, Nres };
+  uint32_t outDimList2[4] = { 3, Nres, Nres, Nres };
+  ProgressiveOutput<float> out_den_1 = 
+    ProgressiveOutput<float>::saveArrayProgressive("density.nc", outDimList, 3);
+  ProgressiveOutput<float> out_vel_1 = 
+    ProgressiveOutput<float>::saveArrayProgressive("v3d.nc", outDimList2, 4);
+  ProgressiveOutput<float> out_rad_1 = 
+    ProgressiveOutput<float>::saveArrayProgressive("rad.nc", outDimList, 3);
+ 
+  cout << "Weighing..." << endl;
+  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)
+		{
+		  continue;
+		}
+
+	      uint32_t numInCell = bins[rx + Nres*ry + Nres*Nres*rz];
+	      if (numInCell > N_SPH)
+		smooth1.fetchNeighbours(c, numInCell);
+	      else
+		smooth1.fetchNeighbours(c);
+	      smooth1.addGridSite(c);
+	    }
+	}
+    }
+
+  
+  
+
+  cout << "Interpolating..." << endl;
+  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_vel_1.put(0);
+		  out_vel_1.put(0);
+		  out_vel_1.put(0);
+		  out_den_1.put(0);
+		  out_rad_1.put(0);
+		  continue;
+		}
+
+	      uint32_t numInCell = bins[rx + ry*Nres + rz*Nres*Nres];
+	      if (numInCell > N_SPH)
+		smooth1.fetchNeighbours(c, numInCell);
+	      else
+		smooth1.fetchNeighbours(c);
+
+	      float val;
+
+	      out_rad_1.put(smooth1.getSmoothingLen());
+	      val = smooth1.computeSmoothedValue(c, getVelocity<0>);
+	      out_vel_1.put(val);
+	      val = smooth1.computeSmoothedValue(c, getVelocity<1>);
+	      out_vel_1.put(val);
+	      val = smooth1.computeSmoothedValue(c, getVelocity<2>);
+	      out_vel_1.put(val);
+	      val = smooth1.computeSmoothedValue(c, getUnity);
+	      out_den_1.put(val);
+	    }
+	}
+    }
+
+  return 0;
+};

src/CMakeLists.txt

@@ -9,6 +9,7 @@ SET(CosmoTool_SRCS
 	miniargs.cpp
 	growthFactor.cpp
 	cosmopower.cpp
+  cic.cpp
 )
 
 IF(FOUND_NETCDF3)

src/sphSmooth.tcc

@@ -47,12 +47,12 @@ SPHSmooth<ValType,Ndims>::fetchNeighbours(const typename SPHTree::coords& c, uin
   if (requested > maxNgb)
     {
       maxNgb = requested;
-      internal.ngb = new SPHCell *[maxNgb];
-      internal.distances = new CoordType[maxNgb];
+      internal.ngb = boost::shared_ptr<SPHCell*[]>(new SPHCell *[maxNgb]);
+      internal.distances = boost::shared_ptr<CoordType[]>(new CoordType[maxNgb]);
     }
 
   memcpy(internal.currentCenter, c, sizeof(c));
-  tree->getNearestNeighbours(c, requested, internal.ngb, internal.distances);
+  tree->getNearestNeighbours(c, requested, internal.ngb.get(), internal.distances.get());
 
   internal.currentNgb = 0;
   for (uint32_t i = 0; i < requested && internal.ngb[i] != 0; i++,internal.currentNgb++)