diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index f67dc6119e49da312adb0299da6ddea443cef817..44dc67b23e83f085e7ac9408dd5d465edc0cedda 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -1321,12 +1321,17 @@ void GFaceCompound::buildOct() const
 //Verify topological conditions for computing the harmonic map
 bool GFaceCompound::checkTopology() const
 {
+  // FIXME!!! I think those things are wrong with cross-patch reparametrization
+  if ((*(_compound.begin()))->geomType() != GEntity::DiscreteSurface)return true;
+  
   bool correctTopo = true;
 
   int Nb = _interior_loops.size();
   int G  = genusGeom() ;
   if( G != 0 || Nb < 1) correctTopo = false;
 
+  //  printf("%d %d\n",Nb,G);
+
   if (!correctTopo){
     Msg::Info("--> Wrong topology: Genus=%d and N boundary=%d", G, Nb);
     nbSplit = std::max(G+1, 2);
diff --git a/Numeric/Numeric.cpp b/Numeric/Numeric.cpp
index ccc75801f4bcad70c1847c29a264c0bdcf9f903f..a9007871b0d317f64fcdd9db99393349957ce7d5 100644
--- a/Numeric/Numeric.cpp
+++ b/Numeric/Numeric.cpp
@@ -718,18 +718,25 @@ void signedDistancesPointsTriangle (std::vector<double>&distances,
   SVector3 t3 = p3-p2;
   SVector3 n = crossprod(t1,t2);
   n.normalize();
+
   double mat[3][3] = {{t1.x(),t2.x(),-n.x()},
 		      {t1.y(),t2.y(),-n.y()},
 		      {t1.z(),t2.z(),-n.z()}};
   double inv[3][3];
-  inv3x3(mat,inv);
+  double det = inv3x3(mat,inv);
   
+  distances.clear();
+  distances.resize(pts.size());
+
+  for (int i=0; i<pts.size();i++)
+    distances[i] = 1.e22;
+
+  if (det = 0.0)return;
+
   const double n2t1 = dot(t1,t1);
   const double n2t2 = dot(t2,t3);
   const double n2t3 = dot(t3,t3);
 
-  distances.clear();
-  distances.resize(pts.size());
 
   double u,v,d;
   for (int i=0; i<pts.size();i++){
@@ -739,7 +746,7 @@ void signedDistancesPointsTriangle (std::vector<double>&distances,
     v = (inv[1][0] * pp1.x() + inv[1][1] * pp1.y() + inv[1][2] * pp1.z());
     d = (inv[2][0] * pp1.x() + inv[2][1] * pp1.y() + inv[2][2] * pp1.z());
     double sign = (d>0) ? 1.0 : -1.0;
-    if (d == 0) sign = 1.e22;
+    if (d == 0) sign = 1.e10;
     
     if (u >= 0 && v >=0 && 1.-u-v >= 0.0){
       distances[i] = d;
@@ -749,7 +756,7 @@ void signedDistancesPointsTriangle (std::vector<double>&distances,
       const double t13 = dot(pp1,t2) / n2t2;
       SVector3 pp2 = p - p2;
       const double t23 = dot(pp2,t3) / n2t3;
-      d = 1.e22;
+      d = 1.e10;
       bool found = false;
       if (t12 >=0 && t12 <= 1.0){
 	d = sign * std::min(fabs(d),p.distance(p1+(p2-p1)*t12));            
diff --git a/Numeric/cartesian.h b/Numeric/cartesian.h
index f2877aff6e3fbb3d3dfee3aaa94f2a5640748ede..a78963ecd67bec0e61bb67bd189ba36143fd509a 100644
--- a/Numeric/cartesian.h
+++ b/Numeric/cartesian.h
@@ -63,7 +63,7 @@ private:
   }  
 
   typename std::map<INDEXTYPE,scalar>::const_iterator begin() const {return _nodal_values.begin();}
-  typename std::map<INDEXTYPE,scalar>::const_iterator ens() const {return _nodal_values.endn();}
+  typename std::map<INDEXTYPE,scalar>::const_iterator end() const {return _nodal_values.end();}
 
   // add that in the ann search tool
   void insert_point (double x, double y, double z){
@@ -72,19 +72,27 @@ private:
   // compute distance
   double distance (double x, double y, double z) const;
   
+  // set the value
+  void setValue(INDEXTYPE i, scalar s){
+    _nodal_values[i] = s;
+  }
   inline INDEXTYPE index_of_element (double x, double y, double z) const{
     
     // P = P_0 + xi * _vdx + eta * _vdy + zeta *vdz
     // DP = P-P_0 * _vdx = xi
     SVector3 DP (x-_X,y-_Y,z-_Z);
 
-    const double xi   = dot(DP,_xiAxis);
-    const double eta  = dot(DP,_etaAxis);
-    const double zeta = dot(DP,_zetaAxis);
-    
-    const int i = xi/_dxi * _Nxi;
-    const int j = eta/_deta * _Neta;
-    const int k = zeta/_dzeta * _Nzeta;
+    double xi   = dot(DP,_xiAxis);
+    double eta  = dot(DP,_etaAxis);
+    double zeta = dot(DP,_zetaAxis);    
+
+    int i = xi/_dxi * _Nxi;
+    int j = eta/_deta * _Neta;
+    int k = zeta/_dzeta * _Nzeta;
+
+    if (i < 0) i = 0;if (i >= _Nxi) i = _Nxi-1;
+    if (j < 0) j = 0;if (j >= _Neta) j = _Neta-1;
+    if (k < 0) k = 0;if (k >= _Nzeta) k = _Nzeta-1;
 
     return element_index(i,j,k);
   }
@@ -179,7 +187,8 @@ private:
       typename std::map<INDEXTYPE,scalar>::const_iterator it = _nodal_values.begin();
       for ( ; it!=_nodal_values.end();++it){
 	SPoint3 p = coordinates_of_node(it->first);
-	fprintf(f,"%d %g\n",it->first,distance(p.x(),p.y(),p.z()));
+	//	fprintf(f,"%d %g\n",it->first,distance(p.x(),p.y(),p.z()));
+	fprintf(f,"%d %g\n",it->first,it->second);
       }    
       fprintf(f,"$EndNodeData\n");
     }