diff --git a/src/dinterpolate.hpp b/src/dinterpolate.hpp index 2777764..5109ce9 100644 --- a/src/dinterpolate.hpp +++ b/src/dinterpolate.hpp @@ -45,6 +45,7 @@ namespace CosmoTool { */ DelaunayInterpolate(CoordType *positions, IType *values, uint32_t *simplex_list, uint32_t numPoints, uint32_t numSimplex) + throw (InvalidArgumentException) { this->positions = positions; this->values = values; @@ -69,7 +70,8 @@ namespace CosmoTool { gsl_eigen_symmv_free(eigen_work); } - void buildPreweight(); + void buildPreweight() + throw (InvalidArgumentException); void buildQuickAccess(); void buildHyperplane(const PType *v, CoordType& hyper); diff --git a/src/dinterpolate.tcc b/src/dinterpolate.tcc index 810e1e7..d545788 100644 --- a/src/dinterpolate.tcc +++ b/src/dinterpolate.tcc @@ -70,6 +70,7 @@ namespace CosmoTool { template void DelaunayInterpolate::buildPreweight() + throw(InvalidArgumentException) { double preweight[N*N]; double preweight_inverse[N*N]; @@ -97,6 +98,8 @@ namespace CosmoTool { int signum; gsl_linalg_LU_decomp(&M.matrix, p, &signum); + if (fabs(gsl_linalg_LU_det(&M.matrix, signum)) < 1e-10) + throw InvalidArgumentException("Invalid tesselation. One simplex is coplanar."); gsl_linalg_LU_invert(&M.matrix, p, &iM.matrix); for (int j = 0; j < N*N; j++) @@ -239,7 +242,7 @@ namespace CosmoTool { // The point does not belong to any simplex. if (i != N_ngb) - throw InvalidArgumentException(); + throw InvalidArgumentException("the given point does not belong to any simplex"); N_ngb *= 2; cell_Ngb = new QuickCell *[N_ngb]; @@ -271,13 +274,13 @@ namespace CosmoTool { for (int j = 0; j < N; j++) weight[i] += preweight[(i-1)*N+j]*(c[j]-p0[j]); - assert(weight[i] >= 0); - assert(weight[i] <= 1); + assert(weight[i] > -1e-7); + assert(weight[i] < 1+1e-7); sum_weight += weight[i]; } weight[0] = 1-sum_weight; - assert(weight[0] >= 0); - assert(weight[0] <= 1); + assert(weight[0] > -1e-7); + assert(weight[0] < (1+1e-7)); // We compute the final value by weighing the value at the N+1 // points by the proper weight.