Changeset 849


Ignore:
Timestamp:
05/11/16 14:56:07 (6 years ago)
Author:
ymipsl
Message:

Remapper : manage some floating rounding error leading to some NAN

YM

Location:
XIOS/trunk/extern/remap/src
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • XIOS/trunk/extern/remap/src/intersection_ym.cpp

    r845 r849  
    144144       is->area = polygonarea(intersectPolygon,nv) ; 
    145145//        if (is->area < 1e-12) cout<<"Small intersection : "<<is->area<<endl ; 
    146        is->id = b->id; /* intersection holds id of corresponding source element (see Elt class definition for details about id) */ 
    147        is->src_id = b->src_id; 
    148        is->n = nv; 
    149        (a->is).push_back(is); 
    150        (b->is).push_back(is); 
    151        area=is->area ; 
     146       if (is->area==0.) delete is ; 
     147       else 
     148       {   
     149         is->id = b->id; /* intersection holds id of corresponding source element (see Elt class definition for details about id) */ 
     150         is->src_id = b->src_id; 
     151         is->n = nv; 
     152         (a->is).push_back(is); 
     153         (b->is).push_back(is); 
     154         area=is->area ; 
     155       } 
    152156    } 
    153157    delete[] intersectPolygon ; 
  • XIOS/trunk/extern/remap/src/polyg.cpp

    r845 r849  
    131131        double s = 0.5 * (a + b + c); 
    132132        double t = tan(0.5*s) * tan(0.5*(s - a)) * tan(0.5*(s - b)) * tan(0.5*(s - c)); 
    133         assert(t >= 0); 
    134 //        if (t<0) return 0. ; 
     133//      assert(t >= 0); 
     134  if (t<1e-20) return 0. ; 
    135135        return 4 * atan(sqrt(t)); 
    136136} 
  • XIOS/trunk/extern/remap/src/triple.cpp

    r688 r849  
    4040double arcdist(const Coord &x, const Coord &y) 
    4141{ // et angles aigus non orientes 
    42         Coord n = crossprod(x, y); 
    43         double a = asin(norm(n)); 
     42        double n = norm(crossprod(x, y)); 
     43  if (n>1.) n=1. ; 
     44        double a = asin(n); 
    4445        if (squaredist(x,y) > 2.0) a = M_PI - a; 
    4546        return a; 
Note: See TracChangeset for help on using the changeset viewer.