guilhem_lavaux/cosmotool
0
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Better error reporting

Browse Source
This commit is contained in:
Guilhem Lavaux 2009-03-08 22:05:15 -07:00
parent c7a5c9c193
commit 9df2340837
2 changed files with 11 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
src/dinterpolate.hpp
View File

@ -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);

13
src/dinterpolate.tcc
View File

@ -70,6 +70,7 @@ namespace CosmoTool {
template<typename PType, typename IType, int N>
void DelaunayInterpolate<PType,IType,N>::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.