fix bug (using hack...) due to ambiguous distance function (same abs value, different sign)

Numeric/cartesian.h

@@ -140,7 +140,8 @@ class cartesianBox {
       _childBox = new cartesianBox<scalar>(X, Y, Z, dxi, deta, dzeta,
                                            2 * Nxi, 2 * Neta, 2 * Nzeta,
                                            level - 1);
-  }  
+  }
+  double getLC(){ return sqrt(_dxi * _dxi + _deta * _deta + _dzeta * _dzeta); }
   int getNxi(){ return _Nxi; }
   int getNeta(){ return _Neta; }
   int getNzeta(){ return _Nzeta; }

utils/api_demos/mainCartesian.cpp

@@ -29,6 +29,8 @@ void insertActiveCells(double x, double y, double z, double rmax,
 
 void computeLevelset(GModel *gm, cartesianBox<double> &box)
 {
+  // tolerance for desambiguation
+  const double tol = box.getLC() * 1.e-12;
   std::vector<SPoint3> nodes;
   std::vector<int> indices;
   for (cartesianBox<double>::valIter it = box.nodalValuesBegin(); 
@@ -60,9 +62,21 @@ void computeLevelset(GModel *gm, cartesianBox<double> &box)
 	signedDistancesPointsTriangle(localdist, dummy, nodes, P2, P1, P3);
       if(dist.empty()) 
         dist = localdist;
-      else 
-        for (unsigned int j = 0; j < localdist.size(); j++)
-          dist[j] = (fabs(dist[j]) < fabs(localdist[j])) ? dist[j] : localdist[j];
+      else{ 
+        for (unsigned int j = 0; j < localdist.size(); j++){
+          // FIXME: if there is an ambiguity assume we are inside (to
+          // avoid holes in the structure). This is definitely just a
+          // hack, as it could create pockets of matter outside the
+          // structure...
+          if(dist[j] * localdist[j] < 0 && 
+             fabs(fabs(dist[j]) - fabs(localdist[j])) < tol){
+            dist[j] = std::max(dist[j], localdist[j]);
+          }
+          else{
+            dist[j] = (fabs(dist[j]) < fabs(localdist[j])) ? dist[j] : localdist[j];
+          }
+        }
+      }
     }
   }
   for (unsigned int j = 0; j < dist.size(); j++)