Changeset 845


Ignore:
Timestamp:
04/27/16 17:25:57 (6 years ago)
Author:
ymipsl
Message:

Correction of area incoherency between supermesh and target cells due to approximation of small circle by piece of great circle.

YM

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

Legend:

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

    r688 r845  
    2424  vector<Coord> srcPolygon ; 
    2525  createGreatCirclePolygon(*b, srcGrid.pole, srcPolygon) ; 
     26//  b->area=polygonarea(&srcPolygon[0],srcPolygon.size()) ; 
    2627  vector<Coord> dstPolygon ; 
    2728  createGreatCirclePolygon(*a, tgtGrid.pole, dstPolygon) ; 
     29  a->area=polygonarea(&dstPolygon[0],dstPolygon.size()) ; // just for target 
    2830 
    2931// compute coordinates of the polygons into the gnomonique plane tangent to barycenter C of dst polygon 
     
    109111 
    110112  double area=0 ; 
    111   if (intersection.size()==1) 
    112   { 
    113 // go back into real coordinate on the sphere 
    114     Coord* intersectPolygon=new Coord[intersection[0].size()] ; 
    115 //    Coord* intersect2D=new Coord[intersection[0].size()] ; 
    116     for(int n=0; n < intersection[0].size(); n++) 
    117     { 
    118       double x=intersection[0][n].X/xscale+xoffset ; 
    119       double y=intersection[0][n].Y/yscale+yoffset ; 
    120 //      intersect2D[n].x=x  ; 
    121 //      intersect2D[n].y=y  ; 
    122 //      intersect2D[n].z=1. ; 
     113 
     114  for(int ni=0;ni<intersection.size(); ni++) 
     115  { 
     116    // go back into real coordinate on the sphere 
     117    Coord* intersectPolygon=new Coord[intersection[ni].size()] ; 
     118    for(int n=0; n < intersection[ni].size(); n++) 
     119    { 
     120      double x=intersection[ni][n].X/xscale+xoffset ; 
     121      double y=intersection[ni][n].Y/yscale+yoffset ; 
    123122 
    124123      intersectPolygon[n]=Ox*x+Oy*y+Oz ; 
     
    128127// remove redondants vertex 
    129128    int nv=0 ; 
    130     for(int n=0; n < intersection[0].size(); n++) 
    131     { 
    132       if (norm(intersectPolygon[n]-intersectPolygon[(n+1)%intersection[0].size()])>fusion_vertex) 
     129    for(int n=0; n < intersection[ni].size(); n++) 
     130    { 
     131      if (norm(intersectPolygon[n]-intersectPolygon[(n+1)%intersection[ni].size()])>fusion_vertex) 
    133132      { 
    134133        intersectPolygon[nv]=intersectPolygon[n] ; 
     
    151150       (b->is).push_back(is); 
    152151       area=is->area ; 
    153      } 
    154      delete[] intersectPolygon ; 
    155   } 
    156   else if (intersection.size()>1) 
    157   { 
    158  
    159     cout<<"Intersection Size > 1 : "<< intersection.size()<<endl ; 
     152    } 
     153    delete[] intersectPolygon ; 
    160154  } 
    161155 
  • XIOS/trunk/extern/remap/src/polyg.cpp

    r688 r845  
    5151        /* both distances can be equal down to roundoff when norm(bc) < mashineepsilon  
    5252           which can occur when weighted with tiny area */ 
     53 
    5354        assert(squaredist(bc, proj(bc)) <= squaredist(bc, proj(bc * (-1.0)))); 
    5455        //if (squaredist(bc, proj(bc)) > squaredist(bc, proj(bc * (-1.0)))) return proj(bc * (-1.0)); 
     
    266267int packIntersectionSize(const Elt& elt)  
    267268{ 
    268         return elt.is.size() * (2*sizeof(int)+ sizeof(GloId) + 4*sizeof(double)); 
     269        return elt.is.size() * (2*sizeof(int)+ sizeof(GloId) + 5*sizeof(double)); 
    269270} 
    270271 
    271272void packIntersection(const Elt& e, char* buffer,int& pos)  
    272273{ 
    273         for (list<Polyg *>::const_iterator it = e.is.begin(); it != e.is.end(); ++it) 
     274  for (list<Polyg *>::const_iterator it = e.is.begin(); it != e.is.end(); ++it) 
    274275        { 
    275276                *((int *) &(buffer[0])) += 1; 
     
    278279                pos += sizeof(int); 
    279280 
     281    *((double *) &(buffer[pos])) = e.area; 
     282    pos += sizeof(double); 
     283 
    280284                *((GloId *) &(buffer[pos])) = (*it)->id; 
    281285                pos += sizeof(GloId); 
    282  
     286   
    283287                *((int *) &(buffer[pos])) = (*it)->n; 
    284288                pos += sizeof(int); 
     
    295299        int ind; 
    296300        int pos = 0; 
     301   
    297302        int n = *((int *) & (buffer[pos])); 
    298303        pos += sizeof(int); 
     
    303308 
    304309                Elt& elt= e[ind]; 
     310 
     311    elt.area=*((double *) & (buffer[pos])); 
     312                pos += sizeof(double); 
     313 
    305314                Polyg *polygon = new Polyg; 
    306315 
Note: See TracChangeset for help on using the changeset viewer.