Changeset 849 for XIOS/trunk/extern

Timestamp:

05/11/16 14:56:07 (8 years ago)

Author:

ymipsl

Message:

Remapper : manage some floating rounding error leading to some NAN

YM

Location:

XIOS/trunk/extern/remap/src

Files:

• intersection_ym.cpp (modified) (1 diff)
• polyg.cpp (modified) (1 diff)
• triple.cpp (modified) (1 diff)

XIOS/trunk/extern/remap/src/intersection_ym.cpp

@@ -144 +144 @@
        is->area = polygonarea(intersectPolygon,nv) ;
 //        if (is->area < 1e-12) cout<<"Small intersection : "<<is->area<<endl ;
-       is->id = b->id; /* intersection holds id of corresponding source element (see Elt class definition for details about id) */
-       is->src_id = b->src_id;
-       is->n = nv;
-       (a->is).push_back(is);
-       (b->is).push_back(is);
-       area=is->area ;
+       if (is->area==0.) delete is ;
+       else
+       {  
+         is->id = b->id; /* intersection holds id of corresponding source element (see Elt class definition for details about id) */
+         is->src_id = b->src_id;
+         is->n = nv;
+         (a->is).push_back(is);
+         (b->is).push_back(is);
+         area=is->area ;
+       }
     }
     delete[] intersectPolygon ;

XIOS/trunk/extern/remap/src/polyg.cpp

@@ -131 +131 @@
         double s = 0.5 * (a + b + c);
         double t = tan(0.5*s) * tan(0.5*(s - a)) * tan(0.5*(s - b)) * tan(0.5*(s - c));
-        assert(t >= 0);
-//        if (t<0) return 0. ;
+//      assert(t >= 0);
+  if (t<1e-20) return 0. ;
         return 4 * atan(sqrt(t));
 }

XIOS/trunk/extern/remap/src/triple.cpp

@@ -40 +40 @@
 double arcdist(const Coord &x, const Coord &y)
 { // et angles aigus non orientes
-        Coord n = crossprod(x, y);
-        double a = asin(norm(n));
+        double n = norm(crossprod(x, y));
+  if (n>1.) n=1. ;
+        double a = asin(n);
         if (squaredist(x,y) > 2.0) a = M_PI - a;
         return a;