Changeset 845 for XIOS

Timestamp:

04/27/16 17:25:57 (8 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:

• intersection_ym.cpp (modified) (4 diffs)
• polyg.cpp (modified) (5 diffs)

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

@@ -24 +24 @@
   vector<Coord> srcPolygon ;
   createGreatCirclePolygon(*b, srcGrid.pole, srcPolygon) ;
+//  b->area=polygonarea(&srcPolygon[0],srcPolygon.size()) ;
   vector<Coord> dstPolygon ;
   createGreatCirclePolygon(*a, tgtGrid.pole, dstPolygon) ;
+  a->area=polygonarea(&dstPolygon[0],dstPolygon.size()) ; // just for target
 
 // compute coordinates of the polygons into the gnomonique plane tangent to barycenter C of dst polygon
@@ -109 +111 @@
 
   double area=0 ;
-  if (intersection.size()==1)
-  {
-// go back into real coordinate on the sphere
-    Coord* intersectPolygon=new Coord[intersection[0].size()] ;
-//    Coord* intersect2D=new Coord[intersection[0].size()] ;
-    for(int n=0; n < intersection[0].size(); n++)
-    {
-      double x=intersection[0][n].X/xscale+xoffset ;
-      double y=intersection[0][n].Y/yscale+yoffset ;
-//      intersect2D[n].x=x  ;
-//      intersect2D[n].y=y  ;
-//      intersect2D[n].z=1. ;
+
+  for(int ni=0;ni<intersection.size(); ni++)
+  {
+    // go back into real coordinate on the sphere
+    Coord* intersectPolygon=new Coord[intersection[ni].size()] ;
+    for(int n=0; n < intersection[ni].size(); n++)
+    {
+      double x=intersection[ni][n].X/xscale+xoffset ;
+      double y=intersection[ni][n].Y/yscale+yoffset ;
 
       intersectPolygon[n]=Ox*x+Oy*y+Oz ;
@@ -128 +127 @@
 // remove redondants vertex
     int nv=0 ;
-    for(int n=0; n < intersection[0].size(); n++)
-    {
-      if (norm(intersectPolygon[n]-intersectPolygon[(n+1)%intersection[0].size()])>fusion_vertex)
+    for(int n=0; n < intersection[ni].size(); n++)
+    {
+      if (norm(intersectPolygon[n]-intersectPolygon[(n+1)%intersection[ni].size()])>fusion_vertex)
       {
         intersectPolygon[nv]=intersectPolygon[n] ;
@@ -151 +150 @@
        (b->is).push_back(is);
        area=is->area ;
-     }
-     delete[] intersectPolygon ;
-  }
-  else if (intersection.size()>1)
-  {
-
-    cout<<"Intersection Size > 1 : "<< intersection.size()<<endl ;
+    }
+    delete[] intersectPolygon ;
   }
 

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

@@ -51 +51 @@
         /* both distances can be equal down to roundoff when norm(bc) < mashineepsilon 
            which can occur when weighted with tiny area */
+
         assert(squaredist(bc, proj(bc)) <= squaredist(bc, proj(bc * (-1.0))));
         //if (squaredist(bc, proj(bc)) > squaredist(bc, proj(bc * (-1.0)))) return proj(bc * (-1.0));
@@ -266 +267 @@
 int packIntersectionSize(const Elt& elt) 
 {
-        return elt.is.size() * (2*sizeof(int)+ sizeof(GloId) + 4*sizeof(double));
+        return elt.is.size() * (2*sizeof(int)+ sizeof(GloId) + 5*sizeof(double));
 }
 
 void packIntersection(const Elt& e, char* buffer,int& pos) 
 {
-        for (list<Polyg *>::const_iterator it = e.is.begin(); it != e.is.end(); ++it)
+  for (list<Polyg *>::const_iterator it = e.is.begin(); it != e.is.end(); ++it)
         {
                 *((int *) &(buffer[0])) += 1;
@@ -278 +279 @@
                 pos += sizeof(int);
 
+    *((double *) &(buffer[pos])) = e.area;
+    pos += sizeof(double);
+
                 *((GloId *) &(buffer[pos])) = (*it)->id;
                 pos += sizeof(GloId);
-
+  
                 *((int *) &(buffer[pos])) = (*it)->n;
                 pos += sizeof(int);
@@ -295 +299 @@
         int ind;
         int pos = 0;
+  
         int n = *((int *) & (buffer[pos]));
         pos += sizeof(int);
@@ -303 +308 @@
 
                 Elt& elt= e[ind];
+
+    elt.area=*((double *) & (buffer[pos]));
+                pos += sizeof(double);
+
                 Polyg *polygon = new Polyg;