Merge branch 'master' of bitbucket.org:glavaux/cosmotool

src/dinterpolate.hpp

@@ -29,6 +29,7 @@ namespace CosmoTool {
     uint32_t numSimplex;
     uint32_t *simplex_list;
     gsl_eigen_symmv_workspace *eigen_work;
+    bool *disable_simplex;
     
     /**
      * This construct the interpolator. The construction is time consuming so
@@ -53,6 +54,7 @@ namespace CosmoTool {
       this->simplex_list = simplex_list;
       this->numPoints = numPoints;
       this->numSimplex = numSimplex;
+      this->disable_simplex = new bool[numSimplex];
       
       buildPreweight();
       buildQuickAccess();
@@ -67,6 +69,7 @@ namespace CosmoTool {
       delete quickAccess;
       delete[] point_to_simplex_list_base;      
       delete[] all_preweight;
+      delete[] disable_simplex;
 
       gsl_eigen_symmv_free(eigen_work);
     }

src/dinterpolate.tcc

@@ -1,3 +1,6 @@
+#include <fstream>
+#include <iostream>
+#include <iomanip>
 #include <cstdlib>
 #include <cassert>
 #include <gsl/gsl_matrix.h>
@@ -21,6 +24,7 @@ namespace CosmoTool {
     for (uint32_t i = 0; i < (N+1)*numSimplex; i++)
       {
 	assert(simplex_list[i] < numPoints);
+	if (!disable_simplex[i/(N+1)])
 	  numSimplex_by_point[simplex_list[i]]++;
       }
 
@@ -37,7 +41,11 @@ namespace CosmoTool {
       {
 	for (int j = 0; j <= N; j++)
 	  {
-	    uint32_t p = simplex_list[(N+1)*i+j];
+	    uint32_t s = (N+1)*i+j;
+	    if (disable_simplex[i])
+	      continue;
+
+	    uint32_t p = simplex_list[s];
 	    assert(index_by_point[p] < point_to_simplex_size);
 	    point_to_simplex_list_base[index_by_point[p]] = i;
 	    ++index_by_point[p];
@@ -48,7 +56,9 @@ namespace CosmoTool {
     for (uint32_t i = 0; i < numPoints; i++)
       {
 	// check assertion
-	assert((i==0 && index_by_point[0]==numSimplex_by_point[0]) || ((index_by_point[i]-index_by_point[i-1]) == (numSimplex_by_point[i]+1)));
+	assert((i==0 && index_by_point[0]==numSimplex_by_point[0]) 
+	       || 
+	       ((index_by_point[i]-index_by_point[i-1]) == (numSimplex_by_point[i]+1)));
 	assert(index_by_point[i] < point_to_simplex_size);
 	point_to_simplex_list_base[index_by_point[i]] = -1;
       }
@@ -80,6 +90,7 @@ namespace CosmoTool {
     double preweight[N*N];
     double preweight_inverse[N*N];
     gsl_permutation *p = gsl_permutation_alloc(N);
+    uint32_t numDisabled = 0;
 
     all_preweight = new PType[N*N*numSimplex];
 
@@ -105,13 +116,38 @@ namespace CosmoTool {
 	gsl_linalg_LU_decomp(&M.matrix, p, &signum);
         double a = fabs(gsl_linalg_LU_det(&M.matrix, signum));
 	if (a < 1e-10)
+         {
+#ifdef DEBUG
+           for (int j = 0; j < N; j++)
+            {
+              PType xref = positions[pref][j];
+
+              for (int k = 0; k < N; k++)
+                {
+                  preweight[j*N + k] = positions[simplex_list[k+base+1]][j] - xref;
+                }
+            }
+           std::ofstream f("matrix.txt");
+           for (int j = 0; j < N*N; j++)
+             f << std::setprecision(12) << preweight[j] << std::endl;
 	   throw InvalidArgumentException("Invalid tesselation. One simplex is coplanar.");
+#else
+           gsl_matrix_set_zero(&iM.matrix);
+	   disable_simplex[i] = true;
+	   numDisabled++;
+#endif
+         }
+        else {
 	  gsl_linalg_LU_invert(&M.matrix, p, &iM.matrix);
+	  disable_simplex[i] = false;
+        }
  
 	for (int j = 0; j < N*N; j++)
 	  all_preweight[N*N*i + j] = preweight_inverse[j];
       }
 
+    std::cout << "Number of disabled simplices: " << numDisabled << std::endl;
+
     gsl_permutation_free(p);
   }
 
@@ -166,6 +202,9 @@ namespace CosmoTool {
   template<typename PType, typename IType, int N>
   bool DelaunayInterpolate<PType,IType,N>::checkPointInSimplex(const CoordType& pos, uint32_t simplex)
   {
+    if (disable_simplex[simplex])
+      return false;
+
     uint32_t *desc_simplex = &simplex_list[simplex*(N+1)];
     CoordType *p[N+1], v[N], hyper;
 

src/loadRamses.cpp

@@ -318,6 +318,9 @@ CosmoTool::SimuData *CosmoTool::loadRamsesSimu(const char *basename, int outputI
   data->Omega_M = info.omega_m;
   data->Omega_Lambda = info.omega_lambda;
 
+  if (flags == 0)
+    return data;
+
   if (cpuid < 0)
     cpuid = 1;
   else cpuid++;