Changeset 1642 for XIOS/dev/branch_openmp/extern

Timestamp:

01/23/19 10:31:44 (5 years ago)

Author:

yushan

Message:

dev on ADA. add flag switch _usingEP/_usingMPI

Location:

XIOS/dev/branch_openmp/extern

Files:

• remap/src/cputime.cpp (modified) (1 diff)
• remap/src/elt.hpp (modified) (2 diffs)
• remap/src/gridRemap.cpp (modified) (1 diff)
• remap/src/gridRemap.hpp (modified) (1 diff)
• remap/src/intersect.cpp (modified) (1 diff)
• remap/src/intersection_ym.cpp (modified) (7 diffs)
• remap/src/libmapper.cpp (modified) (3 diffs)
• remap/src/mapper.cpp (modified) (28 diffs)
• remap/src/mapper.hpp (modified) (1 diff)
• remap/src/meshutil.cpp (modified) (3 diffs)
• remap/src/mpi_cascade.cpp (modified) (2 diffs)
• remap/src/mpi_cascade.hpp (modified) (2 diffs)
• remap/src/mpi_routing.cpp (modified) (25 diffs)
• remap/src/parallel_tree.cpp (modified) (11 diffs)
• remap/src/parallel_tree.hpp (modified) (1 diff)
• remap/src/polyg.cpp (modified) (8 diffs)
• remap/src/polyg.hpp (modified) (1 diff)
• remap/src/timerRemap.cpp (modified) (3 diffs)
• remap/src/timerRemap.hpp (modified) (1 diff)
• src_ep_dev/ep_allgather.cpp (modified) (1 diff)
• src_ep_dev/ep_allgatherv.cpp (modified) (2 diffs)
• src_ep_dev/ep_allocate.cpp (modified) (1 diff)
• src_ep_dev/ep_create.cpp (modified) (2 diffs)
• src_ep_dev/ep_declaration.cpp (modified) (7 diffs)
• src_ep_dev/ep_declaration.hpp (modified) (3 diffs)
• src_ep_dev/ep_exscan.cpp (modified) (21 diffs)
• src_ep_dev/ep_finalize.cpp (modified) (1 diff)
• src_ep_dev/ep_fortran.cpp (modified) (3 diffs)
• src_ep_dev/ep_gatherv.cpp (modified) (1 diff)
• src_ep_dev/ep_intercomm.cpp (modified) (12 diffs)
• src_ep_dev/ep_lib.cpp (modified) (2 diffs)
• src_ep_dev/ep_lib.hpp (modified) (1 diff)
• src_ep_dev/ep_lib_fortran.hpp (modified) (1 diff)
• src_ep_dev/ep_merge.cpp (modified) (2 diffs)
• src_ep_dev/ep_message.cpp (modified) (2 diffs)
• src_ep_dev/ep_mpi.hpp (modified) (1 diff)
• src_ep_dev/ep_probe.cpp (modified) (4 diffs)
• src_ep_dev/ep_recv.cpp (modified) (1 diff)
• src_ep_dev/ep_reduce.cpp (modified) (19 diffs)
• src_ep_dev/ep_reduce_scatter.cpp (modified) (1 diff)
• src_ep_dev/ep_scan.cpp (modified) (42 diffs)
• src_ep_dev/ep_scatter.cpp (modified) (3 diffs)
• src_ep_dev/ep_scatterv.cpp (modified) (5 diffs)
• src_ep_dev/ep_split.cpp (modified) (5 diffs)
• src_ep_dev/ep_type.hpp (modified) (1 diff)
• src_ep_dev/ep_win.cpp (modified) (3 diffs)

XIOS/dev/branch_openmp/extern/remap/src/cputime.cpp

@@ -1 +1 @@
 #include "mpi.hpp"
-using namespace ep_lib;
 
 namespace sphereRemap {

XIOS/dev/branch_openmp/extern/remap/src/elt.hpp

@@ -48 +48 @@
         int n; /* number of vertices */
         double area;
+  double given_area ;
         Coord x; /* barycentre */
 };
@@ -80 +81 @@
                 n    = rhs.n;
                 area = rhs.area;
+                given_area = rhs.given_area;
                 x    = rhs.x;
                 val  = rhs.val;

XIOS/dev/branch_openmp/extern/remap/src/gridRemap.cpp

@@ -11 +11 @@
 
 CRemapGrid srcGrid;
-#pragma omp threadprivate(srcGrid)
-
 CRemapGrid tgtGrid;
-#pragma omp threadprivate(tgtGrid)
 
 }

XIOS/dev/branch_openmp/extern/remap/src/gridRemap.hpp

@@ -14 +14 @@
 Coord readPole(std::istream&);
 
-//extern CRemapGrid srcGrid;
-//extern CRemapGrid tgtGrid;
+extern CRemapGrid srcGrid;
+extern CRemapGrid tgtGrid;
 
 }

XIOS/dev/branch_openmp/extern/remap/src/intersect.cpp

@@ -14 +14 @@
 
 namespace sphereRemap {
-
-extern CRemapGrid srcGrid;
-#pragma omp threadprivate(srcGrid)
-
-extern CRemapGrid tgtGrid;
-#pragma omp threadprivate(tgtGrid)
 
 using namespace std;

XIOS/dev/branch_openmp/extern/remap/src/intersection_ym.cpp

@@ -8 +8 @@
 #include <stdlib.h>
 #include <limits>
+#include <array>
+#include <cstdint> 
+#include "earcut.hpp"
+#include <fstream> 
+
 
 #define epsilon 1e-3  // epsilon distance ratio over side lenght for approximate small circle by great circle
@@ -14 +19 @@
 namespace sphereRemap {
 
-extern CRemapGrid srcGrid;
-#pragma omp threadprivate(srcGrid)
-
-extern CRemapGrid tgtGrid;
-#pragma omp threadprivate(tgtGrid)
-
-
 using namespace std;
 using namespace ClipperLib ;
+        
 
 double intersect_ym(Elt *a, Elt *b)
 {
+
+  using N = uint32_t;
+  using Point = array<double, 2>;
+  vector<Point> vect_points;
+  vector< vector<Point> > polyline;
 
 // transform small circle into piece of great circle if necessary
@@ -52 +56 @@
   double cos_alpha;
 
+  /// vector<p2t::Point*> polyline;
   for(int n=0; n<na;n++)
   {
@@ -58 +63 @@
     a_gno[n].y=scalarprod(dstPolygon[n],Oy)/cos_alpha ;
     a_gno[n].z=scalarprod(dstPolygon[n],Oz)/cos_alpha ; // must be equal to 1
-  }
+
+    vect_points.push_back( array<double, 2>() );
+    vect_points[n][0] = a_gno[n].x;
+    vect_points[n][1] = a_gno[n].y;
+
+  }
+
+  polyline.push_back(vect_points);
+  vector<N> indices_a_gno = mapbox::earcut<N>(polyline);
+  
+  double area_a_gno=0 ;
+  for(int i=0;i<indices_a_gno.size()/3;++i)
+    {
+      Coord x0 = Ox * polyline[0][indices_a_gno[3*i]][0] + Oy* polyline[0][indices_a_gno[3*i]][1] + Oz ;
+      Coord x1 = Ox * polyline[0][indices_a_gno[3*i+1]][0] + Oy* polyline[0][indices_a_gno[3*i+1]][1] + Oz ;
+      Coord x2 = Ox * polyline[0][indices_a_gno[3*i+2]][0] + Oy* polyline[0][indices_a_gno[3*i+2]][1] + Oz ;
+      area_a_gno+=triarea(x0 * (1./norm(x0)),x1* (1./norm(x1)), x2* (1./norm(x2))) ;
+    }
+
+  vect_points.clear();
+  polyline.clear();
+  indices_a_gno.clear();
+
+ 
 
   for(int n=0; n<nb;n++)
@@ -66 +94 @@
     b_gno[n].y=scalarprod(srcPolygon[n],Oy)/cos_alpha ;
     b_gno[n].z=scalarprod(srcPolygon[n],Oz)/cos_alpha ; // must be equal to 1
-  }
-
+
+    vect_points.push_back( array<double, 2>() );
+    vect_points[n][0] = b_gno[n].x;
+    vect_points[n][1] = b_gno[n].y;
+  }
+
+
+  polyline.push_back(vect_points);
+  vector<N> indices_b_gno = mapbox::earcut<N>(polyline);
+
+  double area_b_gno=0 ;
+  for(int i=0;i<indices_b_gno.size()/3;++i)
+    {
+      Coord x0 = Ox * polyline[0][indices_b_gno[3*i]][0] + Oy* polyline[0][indices_b_gno[3*i]][1] + Oz ;
+      Coord x1 = Ox * polyline[0][indices_b_gno[3*i+1]][0] + Oy* polyline[0][indices_b_gno[3*i+1]][1] + Oz ;
+      Coord x2 = Ox * polyline[0][indices_b_gno[3*i+2]][0] + Oy* polyline[0][indices_b_gno[3*i+2]][1] + Oz ;
+      area_b_gno+=triarea(x0 * (1./norm(x0)),x1* (1./norm(x1)), x2* (1./norm(x2))) ;
+    }
+
+  vect_points.clear();
+  polyline.clear();
+  indices_b_gno.clear();
 
 
@@ -116 +164 @@
   clip.AddPaths(dst, ptClip, true);
   clip.Execute(ctIntersection, intersection);
-
+ 
   double area=0 ;
 
   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 ;
-      intersectPolygon[n]=intersectPolygon[n]*(1./norm(intersectPolygon[n])) ;
-    }
-
-// remove redondants vertex
-    int nv=0 ;
+      intersectPolygon[n].x=intersection[ni][n].X/xscale+xoffset ;
+      intersectPolygon[n].y=intersection[ni][n].Y/yscale+yoffset ;
+    }
+    
+
+    int nv=0;
+
     for(int n=0; n < intersection[ni].size(); n++)
     {
-      if (norm(intersectPolygon[n]-intersectPolygon[(n+1)%intersection[ni].size()])>fusion_vertex)
+       double dx=intersectPolygon[n].x-intersectPolygon[(n+1)%intersection[ni].size()].x ;
+       double dy=intersectPolygon[n].y-intersectPolygon[(n+1)%intersection[ni].size()].y ;
+     
+       if (dx*dx+dy*dy>fusion_vertex*fusion_vertex)
+       {
+          intersectPolygon[nv]=intersectPolygon[n] ;
+
+          vect_points.push_back( array<double, 2>() );
+          vect_points[nv][0] = intersectPolygon[n].x;
+          vect_points[nv][1] = intersectPolygon[n].y;
+
+          nv++ ;
+       }
+      
+
+    }
+
+    polyline.push_back(vect_points);
+    vect_points.clear();
+
+    if (nv>2)
+    {
+ 
+      vector<N> indices = mapbox::earcut<N>(polyline);
+
+      double area2=0 ;
+      for(int i=0;i<indices.size()/3;++i)
       {
-        intersectPolygon[nv]=intersectPolygon[n] ;
-        nv++ ;
+          Coord x0 = Ox * polyline[0][indices[3*i]][0] + Oy* polyline[0][indices[3*i]][1] + Oz ;
+          Coord x1 = Ox * polyline[0][indices[3*i+1]][0] + Oy* polyline[0][indices[3*i+1]][1] + Oz ;
+          Coord x2 = Ox * polyline[0][indices[3*i+2]][0] + Oy* polyline[0][indices[3*i+2]][1] + Oz ;
+          area2+=triarea(x0 * (1./norm(x0)),x1* (1./norm(x1)), x2* (1./norm(x2))) ;
       }
-    }
-
-
-    if (nv>2)
-    {
+
+      polyline.clear();
+
+      for(int n=0; n < nv; n++)
+      {
+        intersectPolygon[n] = Ox*intersectPolygon[n].x+Oy*intersectPolygon[n].y+Oz;
+        intersectPolygon[n] = intersectPolygon[n]*(1./norm(intersectPolygon[n])) ;
+      }
+
+
 //     assign intersection to source and destination polygons
        Polyg *is = new Polyg;
        is->x = exact_barycentre(intersectPolygon,nv);
-       is->area = polygonarea(intersectPolygon,nv) ;
+//       is->area = polygonarea(intersectPolygon,nv) ;
+       is->area = area2 ;
+
 //        if (is->area < 1e-12) cout<<"Small intersection : "<<is->area<<endl ;
        if (is->area==0.) delete is ;
@@ -168 +248 @@
   delete[] b_gno ;
   return area ;
+
 }
+
+
 
 void createGreatCirclePolygon(const Elt& element, const Coord& pole, vector<Coord>& coordinates)

XIOS/dev/branch_openmp/extern/remap/src/libmapper.cpp

@@ -15 +15 @@
 #include "cputime.hpp" // cputime
 
-using namespace ep_lib;
-
 using namespace sphereRemap ;
 
@@ -22 +20 @@
    and deallocated during the second step (computing the weights) */
 Mapper *mapper;
-#pragma omp threadprivate(mapper)
+
 
 /** xxx_bounds_yyy is of length n_vert_per_cell_xxx*n_cell_xxx
@@ -43 +41 @@
         assert(n_cell_dst >= 4);
         assert(1 <= order && order <= 2);
-
-  mapper = new Mapper(MPI_COMM_WORLD);
+  double* src_area=NULL ;
+  double* dst_area=NULL ;
+  mapper = new Mapper(EP_COMM_WORLD);
   mapper->setVerbosity(PROGRESS) ;
-  mapper->setSourceMesh(src_bounds_lon, src_bounds_lat, n_vert_per_cell_src, n_cell_src, src_pole ) ;
-  mapper->setTargetMesh(dst_bounds_lon, dst_bounds_lat, n_vert_per_cell_dst, n_cell_dst, dst_pole ) ;
+  mapper->setSourceMesh(src_bounds_lon, src_bounds_lat, src_area, n_vert_per_cell_src, n_cell_src, src_pole ) ;
+  mapper->setTargetMesh(dst_bounds_lon, dst_bounds_lat, dst_area, n_vert_per_cell_dst, n_cell_dst, dst_pole ) ;
 
 /*

XIOS/dev/branch_openmp/extern/remap/src/mapper.cpp

@@ -12 +12 @@
 
 #include "mapper.hpp"
-using namespace ep_lib;
 
 namespace sphereRemap {
-
-extern CRemapGrid srcGrid;
-#pragma omp threadprivate(srcGrid)
-
-extern CRemapGrid tgtGrid;
-#pragma omp threadprivate(tgtGrid)
-
 
 /* A subdivition of an array into N sub-arays
@@ -35 +27 @@
 
 
-void Mapper::setSourceMesh(const double* boundsLon, const double* boundsLat, int nVertex, int nbCells, const double* pole, const long int* globalId)
+void Mapper::setSourceMesh(const double* boundsLon, const double* boundsLat, const double* area, int nVertex, int nbCells, const double* pole, const long int* globalId)
 {
   srcGrid.pole = Coord(pole[0], pole[1], pole[2]);
 
   int mpiRank, mpiSize;
-  MPI_Comm_rank(communicator, &mpiRank);
-  MPI_Comm_size(communicator, &mpiSize);
+  ep_lib::MPI_Comm_rank(communicator, &mpiRank);
+  ep_lib::MPI_Comm_size(communicator, &mpiSize);
 
   sourceElements.reserve(nbCells);
@@ -51 +43 @@
     long int offset ;
     long int nb=nbCells ;
-    MPI_Scan(&nb,&offset,1,MPI_LONG,MPI_SUM,communicator) ;
+    ep_lib::MPI_Scan(&nb,&offset,1,EP_LONG,EP_SUM,communicator) ;
     offset=offset-nb ;
     for(int i=0;i<nbCells;i++) sourceGlobalId[i]=offset+i ;
@@ -67 +59 @@
     sourceMesh.push_back(Node(elt.x, cptRadius(elt), &sourceElements.back()));
     cptEltGeom(sourceElements[i], Coord(pole[0], pole[1], pole[2]));
-  }
-
-}
-
-void Mapper::setTargetMesh(const double* boundsLon, const double* boundsLat, int nVertex, int nbCells, const double* pole, const long int* globalId)
+    if (area!=NULL) sourceElements[i].given_area=area[i] ;
+    else sourceElements[i].given_area=sourceElements[i].area ;
+  }
+
+}
+
+void Mapper::setTargetMesh(const double* boundsLon, const double* boundsLat, const double* area, int nVertex, int nbCells, const double* pole, const long int* globalId)
 {
   tgtGrid.pole = Coord(pole[0], pole[1], pole[2]);
 
   int mpiRank, mpiSize;
-  MPI_Comm_rank(communicator, &mpiRank);
-  MPI_Comm_size(communicator, &mpiSize);
+  ep_lib::MPI_Comm_rank(communicator, &mpiRank);
+  ep_lib::MPI_Comm_size(communicator, &mpiSize);
 
   targetElements.reserve(nbCells);
@@ -87 +81 @@
     long int offset ;
     long int nb=nbCells ;
-    MPI_Scan(&nb,&offset,1,MPI_LONG,MPI_SUM,communicator) ;
+    ep_lib::MPI_Scan(&nb,&offset,1,EP_LONG,EP_SUM,communicator) ;
     offset=offset-nb ;
     for(int i=0;i<nbCells;i++) targetGlobalId[i]=offset+i ;
@@ -100 +94 @@
     targetMesh.push_back(Node(elt.x, cptRadius(elt), &sourceElements.back()));
     cptEltGeom(targetElements[i], Coord(pole[0], pole[1], pole[2]));
+    if (area!=NULL) targetElements[i].given_area=area[i] ;
+    else targetElements[i].given_area=targetElements[i].area ;
   }
 
@@ -121 +117 @@
   vector<double> timings;
   int mpiSize, mpiRank;
-  MPI_Comm_size(communicator, &mpiSize);
-  MPI_Comm_rank(communicator, &mpiRank);
+  ep_lib::MPI_Comm_size(communicator, &mpiSize);
+  ep_lib::MPI_Comm_rank(communicator, &mpiRank);
 
   this->buildSSTree(sourceMesh, targetMesh);
@@ -132 +128 @@
 
   tic = cputime();
-  if (interpOrder == 2) 
-  {
+  if (interpOrder == 2) {
     if (mpiRank == 0 && verbose) cout << "Computing grads ..." << endl;
     buildMeshTopology();
@@ -142 +137 @@
   /* Prepare computation of weights */
   /* compute number of intersections which for the first order case
-     corresponds to the number of edges in the remap matrix */
+           corresponds to the number of edges in the remap matrix */
   int nIntersections = 0;
   for (int j = 0; j < targetElements.size(); j++)
@@ -178 +173 @@
 {
   int mpiSize, mpiRank;
-  MPI_Comm_size(communicator, &mpiSize);
-  MPI_Comm_rank(communicator, &mpiRank);
+  ep_lib::MPI_Comm_size(communicator, &mpiSize);
+  ep_lib::MPI_Comm_rank(communicator, &mpiRank);
 
   /* create list of intersections (super mesh elements) for each rank */
@@ -187 +182 @@
     Elt& e = elements[j];
     for (list<Polyg *>::iterator it = e.is.begin(); it != e.is.end(); it++)
-    {
       elementList[(*it)->id.rank].insert(pair<int, Polyg *>((*it)->id.ind, *it));
-    }
   }
 
@@ -196 +189 @@
   double **recvValue = new double*[mpiSize];
   double **recvArea = new double*[mpiSize];
+  double **recvGivenArea = new double*[mpiSize];
   Coord **recvGrad = new Coord*[mpiSize];
   GloId **recvNeighIds = new GloId*[mpiSize]; /* ids of the of the source neighbours which also contribute through gradient */
@@ -218 +212 @@
       recvValue[rank]   = new double[nbSendElement[rank]];
       recvArea[rank]    = new double[nbSendElement[rank]];
+      recvGivenArea[rank] = new double[nbSendElement[rank]];
       if (order == 2)
       {
@@ -240 +235 @@
   /* communicate sizes of source elements to be sent (index lists and later values and gradients) */
   int *nbRecvElement = new int[mpiSize];
-  MPI_Alltoall(nbSendElement, 1, MPI_INT, nbRecvElement, 1, MPI_INT, communicator);
-  
+  ep_lib::MPI_Alltoall(nbSendElement, 1, EP_INT, nbRecvElement, 1, EP_INT, communicator);
 
   /* communicate indices of source elements on other ranks whoes value and gradient we need (since intersection) */
@@ -249 +243 @@
   double **sendValue = new double*[mpiSize];
   double **sendArea = new double*[mpiSize];
+  double **sendGivenArea = new double*[mpiSize];
   Coord **sendGrad = new Coord*[mpiSize];
   GloId **sendNeighIds = new GloId*[mpiSize];
-  MPI_Request *sendRequest = new MPI_Request[4*mpiSize];
-  MPI_Request *recvRequest = new MPI_Request[4*mpiSize];
+  ep_lib::MPI_Request *sendRequest = new ep_lib::MPI_Request[5*mpiSize];
+  ep_lib::MPI_Request *recvRequest = new ep_lib::MPI_Request[5*mpiSize];
   for (int rank = 0; rank < mpiSize; rank++)
   {
     if (nbSendElement[rank] > 0)
     {
-      MPI_Issend(sendElement[rank], nbSendElement[rank], MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]);
+      ep_lib::MPI_Issend(sendElement[rank], nbSendElement[rank], EP_INT, rank, 0, communicator, &sendRequest[nbSendRequest]);
       nbSendRequest++;
     }
@@ -266 +261 @@
       sendValue[rank]   = new double[nbRecvElement[rank]];
       sendArea[rank]   = new double[nbRecvElement[rank]];
+      sendGivenArea[rank] = new double[nbRecvElement[rank]];
       if (order == 2)
       {
@@ -275 +271 @@
         sendNeighIds[rank] = new GloId[nbRecvElement[rank]];
       }
-      MPI_Irecv(recvElement[rank], nbRecvElement[rank], MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]);
+      ep_lib::MPI_Irecv(recvElement[rank], nbRecvElement[rank], EP_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]);
       nbRecvRequest++;
     }
   }
-  MPI_Status *status = new MPI_Status[4*mpiSize];
-
-  MPI_Waitall(nbSendRequest, sendRequest, status);
-  MPI_Waitall(nbRecvRequest, recvRequest, status);
+  ep_lib::MPI_Status *status = new ep_lib::MPI_Status[5*mpiSize];
+  
+  ep_lib::MPI_Waitall(nbSendRequest, sendRequest, status);
+        ep_lib::MPI_Waitall(nbRecvRequest, recvRequest, status);
 
   /* for all indices that have been received from requesting ranks: pack values and gradients, then send */
@@ -296 +292 @@
         sendValue[rank][j] = sstree.localElements[recvElement[rank][j]].val;
         sendArea[rank][j] = sstree.localElements[recvElement[rank][j]].area;
+        sendGivenArea[rank][j] = sstree.localElements[recvElement[rank][j]].given_area;
         if (order == 2)
         {
           sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].grad;
+//          cout<<"grad  "<<jj<<"  "<<recvElement[rank][j]<<"  "<<sendGrad[rank][jj]<<" "<<sstree.localElements[recvElement[rank][j]].grad<<endl ;
           sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].src_id;
           jj++;
@@ -304 +302 @@
           {
             sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].gradNeigh[i];
+//            cout<<"grad  "<<jj<<"  "<<sendGrad[rank][jj]<<" "<<sstree.localElements[recvElement[rank][j]].grad<<endl ;
             sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].neighId[i];
-                                                jj++;
-                                        }
-                                }
-                                else
-                                        sendNeighIds[rank][j] = sstree.localElements[recvElement[rank][j]].src_id;
-                        }
-                        MPI_Issend(sendValue[rank],  nbRecvElement[rank], MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]);
-                        nbSendRequest++;
-                        MPI_Issend(sendArea[rank],  nbRecvElement[rank], MPI_DOUBLE, rank, 1, communicator, &sendRequest[nbSendRequest]);
-                        nbSendRequest++;
-                        if (order == 2)
-                        {
-                                MPI_Issend(sendGrad[rank], 3*nbRecvElement[rank]*(NMAX+1),
-                                                                MPI_DOUBLE, rank, 2, communicator, &sendRequest[nbSendRequest]);
-                                nbSendRequest++;
-                                MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank]*(NMAX+1), MPI_INT, rank, 3, communicator, &sendRequest[nbSendRequest]);
+            jj++;
+          }
+        }
+        else
+          sendNeighIds[rank][j] = sstree.localElements[recvElement[rank][j]].src_id;
+      }
+      ep_lib::MPI_Issend(sendValue[rank],  nbRecvElement[rank], EP_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]);
+      nbSendRequest++;
+      ep_lib::MPI_Issend(sendArea[rank],  nbRecvElement[rank], EP_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]);
+      nbSendRequest++;
+      ep_lib::MPI_Issend(sendGivenArea[rank],  nbRecvElement[rank], EP_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]);
+      nbSendRequest++;
+      if (order == 2)
+      {
+        ep_lib::MPI_Issend(sendGrad[rank], 3*nbRecvElement[rank]*(NMAX+1), EP_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]);
+        nbSendRequest++;
+        ep_lib::MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank]*(NMAX+1), EP_INT, rank, 0, communicator, &sendRequest[nbSendRequest]);
 //ym  --> attention taille GloId
-                                nbSendRequest++;
-                        }
-                        else
-                        {
-                                MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank], MPI_INT, rank, 4, communicator, &sendRequest[nbSendRequest]);
+        nbSendRequest++;
+      }
+      else
+      {
+        ep_lib::MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank], EP_INT, rank, 0, communicator, &sendRequest[nbSendRequest]);
 //ym  --> attention taille GloId
-                                nbSendRequest++;
-                        }
-                }
-                if (nbSendElement[rank] > 0)
-                {
-                        MPI_Irecv(recvValue[rank],  nbSendElement[rank], MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]);
-                        nbRecvRequest++;
-                        MPI_Irecv(recvArea[rank],  nbSendElement[rank], MPI_DOUBLE, rank, 1, communicator, &recvRequest[nbRecvRequest]);
-                        nbRecvRequest++;
-                        if (order == 2)
-                        {
-                                MPI_Irecv(recvGrad[rank], 3*nbSendElement[rank]*(NMAX+1),
-                                                MPI_DOUBLE, rank, 2, communicator, &recvRequest[nbRecvRequest]);
-                                nbRecvRequest++;
-                                MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank]*(NMAX+1), MPI_INT, rank, 3, communicator, &recvRequest[nbRecvRequest]);
-//ym  --> attention taille GloId
-                                nbRecvRequest++;
-                        }
-                        else
-                        {
-                                MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank], MPI_INT, rank, 4, communicator, &recvRequest[nbRecvRequest]);
+        nbSendRequest++;
+      }
+    }
+    if (nbSendElement[rank] > 0)
+    {
+      ep_lib::MPI_Irecv(recvValue[rank],  nbSendElement[rank], EP_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]);
+      nbRecvRequest++;
+      ep_lib::MPI_Irecv(recvArea[rank],  nbSendElement[rank], EP_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]);
+      nbRecvRequest++;
+      ep_lib::MPI_Irecv(recvGivenArea[rank],  nbSendElement[rank], EP_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]);
+      nbRecvRequest++;
+      if (order == 2)
+      {
+        ep_lib::MPI_Irecv(recvGrad[rank], 3*nbSendElement[rank]*(NMAX+1),
+            EP_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]);
+        nbRecvRequest++;
+        ep_lib::MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank]*(NMAX+1), EP_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]);
 //ym  --> attention taille GloId
         nbRecvRequest++;
       }
+      else
+      {
+        ep_lib::MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank], EP_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]);
+//ym  --> attention taille GloId
+        nbRecvRequest++;
+      }
     }
   }
         
-  MPI_Waitall(nbSendRequest, sendRequest, status);
-  MPI_Waitall(nbRecvRequest, recvRequest, status);
-        
+        ep_lib::MPI_Waitall(nbSendRequest, sendRequest, status);
+  ep_lib::MPI_Waitall(nbRecvRequest, recvRequest, status);
+  
 
   /* now that all values and gradients are available use them to computed interpolated values on target
-    and also to compute weights */
+     and also to compute weights */
   int i = 0;
   for (int j = 0; j < nbElements; j++)
@@ -367 +369 @@
 
     /* since for the 2nd order case source grid elements can contribute to a destination grid element over several "paths"
-    (step1: gradient is computed using neighbours on same grid, step2: intersection uses several elements on other grid)
-    accumulate them so that there is only one final weight between two elements */
+       (step1: gradient is computed using neighbours on same grid, step2: intersection uses several elements on other grid)
+       accumulate them so that there is only one final weight between two elements */
     map<GloId,double> wgt_map;
 
@@ -380 +382 @@
       double fk = recvValue[rank][n1];
       double srcArea = recvArea[rank][n1];
+      double srcGivenArea = recvGivenArea[rank][n1];
       double w = (*it)->area;
       if (quantity) w/=srcArea ;
+      else w=w*srcGivenArea/srcArea*e.area/e.given_area ;
 
       /* first order: src value times weight (weight = supermesh area), later divide by target area */
@@ -402 +406 @@
     double renorm=0;
     if (renormalize) 
-      for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++) renorm+=it->second / e.area;
+    {
+      if (quantity) for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++) renorm+=it->second ;
+      else for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++) renorm+=it->second / e.area;
+    }
     else renorm=1. ;
 
@@ -417 +424 @@
     }
   }
-        
-        //MPI_Barrier(communicator);
 
   /* free all memory allocated in this function */
@@ -428 +433 @@
       delete[] recvValue[rank];
       delete[] recvArea[rank];
+      delete[] recvGivenArea[rank];
       if (order == 2)
       {
@@ -439 +445 @@
       delete[] sendValue[rank];
       delete[] sendArea[rank];
+      delete[] sendGivenArea[rank];
       if (order == 2)
         delete[] sendGrad[rank];
@@ -463 +470 @@
 void Mapper::computeGrads()
 {
-        /* array of pointers to collect local elements and elements received from other cpu */
-        vector<Elt*> globalElements(sstree.nbLocalElements + nbNeighbourElements);
-        int index = 0;
-        for (int i = 0; i < sstree.nbLocalElements; i++, index++)
-                globalElements[index] = &(sstree.localElements[i]);
-        for (int i = 0; i < nbNeighbourElements; i++, index++)
-                globalElements[index] = &neighbourElements[i];
-
-        update_baryc(sstree.localElements, sstree.nbLocalElements);
-        computeGradients(&globalElements[0], sstree.nbLocalElements);
+  /* array of pointers to collect local elements and elements received from other cpu */
+  vector<Elt*> globalElements(sstree.nbLocalElements + nbNeighbourElements);
+  int index = 0;
+  for (int i = 0; i < sstree.nbLocalElements; i++, index++)
+    globalElements[index] = &(sstree.localElements[i]);
+  for (int i = 0; i < nbNeighbourElements; i++, index++)
+    globalElements[index] = &neighbourElements[i];
+
+  update_baryc(sstree.localElements, sstree.nbLocalElements);
+  computeGradients(&globalElements[0], sstree.nbLocalElements);
 }
 
@@ -479 +486 @@
 void Mapper::buildMeshTopology()
 {
-        int mpiSize, mpiRank;
-        MPI_Comm_size(communicator, &mpiSize);
-        MPI_Comm_rank(communicator, &mpiRank);
-
-        vector<Node> *routingList = new vector<Node>[mpiSize];
-        vector<vector<int> > routes(sstree.localTree.leafs.size());
-
-        sstree.routeIntersections(routes, sstree.localTree.leafs);
-
-        for (int i = 0; i < routes.size(); ++i)
-                for (int k = 0; k < routes[i].size(); ++k)
-                        routingList[routes[i][k]].push_back(sstree.localTree.leafs[i]);
-        routingList[mpiRank].clear();
-
-
-        CMPIRouting mpiRoute(communicator);
-        mpiRoute.init(routes);
-        int nRecv = mpiRoute.getTotalSourceElement();
-
-        int *nbSendNode = new int[mpiSize];
-        int *nbRecvNode = new int[mpiSize];
-        int *sendMessageSize = new int[mpiSize];
-        int *recvMessageSize = new int[mpiSize];
-
-        for (int rank = 0; rank < mpiSize; rank++)
-        {
-                nbSendNode[rank] = routingList[rank].size();
-                sendMessageSize[rank] = 0;
-                for (size_t j = 0; j < routingList[rank].size(); j++)
-                {
-                        Elt *elt = (Elt *) (routingList[rank][j].data);
-                        sendMessageSize[rank] += packedPolygonSize(*elt);
-                }
-        }
-
-        MPI_Alltoall(nbSendNode, 1, MPI_INT, nbRecvNode, 1, MPI_INT, communicator);
-        MPI_Alltoall(sendMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator);
-
-        char **sendBuffer = new char*[mpiSize];
-        char **recvBuffer = new char*[mpiSize];
-        int *pos = new int[mpiSize];
-
-        for (int rank = 0; rank < mpiSize; rank++)
-        {
-                if (nbSendNode[rank] > 0) sendBuffer[rank] = new char[sendMessageSize[rank]];
-                if (nbRecvNode[rank] > 0) recvBuffer[rank] = new char[recvMessageSize[rank]];
-        }
-
-        for (int rank = 0; rank < mpiSize; rank++)
-        {
-                pos[rank] = 0;
-                for (size_t j = 0; j < routingList[rank].size(); j++)
-                {
-                        Elt *elt = (Elt *) (routingList[rank][j].data);
-                        packPolygon(*elt, sendBuffer[rank], pos[rank]);
-                }
-        }
-        delete [] routingList;
-
-
-        int nbSendRequest = 0;
-        int nbRecvRequest = 0;
-        MPI_Request *sendRequest = new MPI_Request[mpiSize];
-        MPI_Request *recvRequest = new MPI_Request[mpiSize];
-        MPI_Status  *status      = new MPI_Status[mpiSize];
-
-        for (int rank = 0; rank < mpiSize; rank++)
-        {
-                if (nbSendNode[rank] > 0)
-                {
-                        MPI_Issend(sendBuffer[rank], sendMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]);
-                        nbSendRequest++;
-                }
-                if (nbRecvNode[rank] > 0)
-                {
-                        MPI_Irecv(recvBuffer[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]);
-                        nbRecvRequest++;
-                }
-        }
-
-        MPI_Waitall(nbRecvRequest, recvRequest, status);
-        MPI_Waitall(nbSendRequest, sendRequest, status);
-
-        for (int rank = 0; rank < mpiSize; rank++)
-                if (nbSendNode[rank] > 0) delete [] sendBuffer[rank];
-        delete [] sendBuffer;
-
-        char **sendBuffer2 = new char*[mpiSize];
-        char **recvBuffer2 = new char*[mpiSize];
-
-        for (int rank = 0; rank < mpiSize; rank++)
-        {
-                nbSendNode[rank] = 0;
-                sendMessageSize[rank] = 0;
-
-                if (nbRecvNode[rank] > 0)
-                {
-                        set<NodePtr> neighbourList;
-                        pos[rank] = 0;
-                        for (int j = 0; j < nbRecvNode[rank]; j++)
-                        {
-                                Elt elt;
-                                unpackPolygon(elt, recvBuffer[rank], pos[rank]);
-                                Node node(elt.x, cptRadius(elt), &elt);
-                                findNeighbour(sstree.localTree.root, &node, neighbourList);
-                        }
-                        nbSendNode[rank] = neighbourList.size();
-                        for (set<NodePtr>::iterator it = neighbourList.begin(); it != neighbourList.end(); it++)
-                        {
-                                Elt *elt = (Elt *) ((*it)->data);
-                                sendMessageSize[rank] += packedPolygonSize(*elt);
-                        }
-
-                        sendBuffer2[rank] = new char[sendMessageSize[rank]];
-                        pos[rank] = 0;
-
-                        for (set<NodePtr>::iterator it = neighbourList.begin(); it != neighbourList.end(); it++)
-                        {
-                                Elt *elt = (Elt *) ((*it)->data);
-                                packPolygon(*elt, sendBuffer2[rank], pos[rank]);
-                        }
-                }
-        }
-        for (int rank = 0; rank < mpiSize; rank++)
-                if (nbRecvNode[rank] > 0) delete [] recvBuffer[rank];
-        delete [] recvBuffer;
-
-
-        MPI_Barrier(communicator);
-        MPI_Alltoall(nbSendNode, 1, MPI_INT, nbRecvNode, 1, MPI_INT, communicator);
-        MPI_Alltoall(sendMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator);
-
-        for (int rank = 0; rank < mpiSize; rank++)
-                if (nbRecvNode[rank] > 0) recvBuffer2[rank] = new char[recvMessageSize[rank]];
-
-        nbSendRequest = 0;
-        nbRecvRequest = 0;
-
-        for (int rank = 0; rank < mpiSize; rank++)
-        {
-                if (nbSendNode[rank] > 0)
-                {
-                        MPI_Issend(sendBuffer2[rank], sendMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]);
-                        nbSendRequest++;
-                }
-                if (nbRecvNode[rank] > 0)
-                {
-                        MPI_Irecv(recvBuffer2[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]);
-                        nbRecvRequest++;
-                }
-        }
-
-        MPI_Waitall(nbRecvRequest, recvRequest, status);
-        MPI_Waitall(nbSendRequest, sendRequest, status);
-
-        int nbNeighbourNodes = 0;
-        for (int rank = 0; rank < mpiSize; rank++)
-                nbNeighbourNodes += nbRecvNode[rank];
-
-        neighbourElements = new Elt[nbNeighbourNodes];
-        nbNeighbourElements = nbNeighbourNodes;
-
-        int index = 0;
-        for (int rank = 0; rank < mpiSize; rank++)
-        {
-                pos[rank] = 0;
-                for (int j = 0; j < nbRecvNode[rank]; j++)
-                {
-                        unpackPolygon(neighbourElements[index], recvBuffer2[rank], pos[rank]);
-                        neighbourElements[index].id.ind = sstree.localTree.leafs.size() + index;
-                        index++;
-                }
-        }
-        for (int rank = 0; rank < mpiSize; rank++)
-        {
-                if (nbRecvNode[rank] > 0) delete [] recvBuffer2[rank];
-                if (nbSendNode[rank] > 0) delete [] sendBuffer2[rank];
-        }
-        delete [] recvBuffer2;
-        delete [] sendBuffer2;
-        delete [] sendMessageSize;
-        delete [] recvMessageSize;
-        delete [] nbSendNode;
-        delete [] nbRecvNode;
-        delete [] sendRequest;
-        delete [] recvRequest;
-        delete [] status;
-        delete [] pos;
-
-        /* re-compute on received elements to avoid having to send this information */
-        neighbourNodes.resize(nbNeighbourNodes);
-        setCirclesAndLinks(neighbourElements, neighbourNodes);
-        cptAllEltsGeom(neighbourElements, nbNeighbourNodes, srcGrid.pole);
-
-        /* the local SS tree must include nodes from other cpus if they are potential
+  int mpiSize, mpiRank;
+  ep_lib::MPI_Comm_size(communicator, &mpiSize);
+  ep_lib::MPI_Comm_rank(communicator, &mpiRank);
+
+  vector<Node> *routingList = new vector<Node>[mpiSize];
+  vector<vector<int> > routes(sstree.localTree.leafs.size());
+
+  sstree.routeIntersections(routes, sstree.localTree.leafs);
+
+  for (int i = 0; i < routes.size(); ++i)
+    for (int k = 0; k < routes[i].size(); ++k)
+      routingList[routes[i][k]].push_back(sstree.localTree.leafs[i]);
+  routingList[mpiRank].clear();
+
+
+  CMPIRouting mpiRoute(communicator);
+  mpiRoute.init(routes);
+  int nRecv = mpiRoute.getTotalSourceElement();
+// cout << mpiRank << " NRECV " << nRecv << "(" << routes.size() << ")"<< endl;
+
+  int *nbSendNode = new int[mpiSize];
+  int *nbRecvNode = new int[mpiSize];
+  int *sendMessageSize = new int[mpiSize];
+  int *recvMessageSize = new int[mpiSize];
+
+  for (int rank = 0; rank < mpiSize; rank++)
+  {
+    nbSendNode[rank] = routingList[rank].size();
+    sendMessageSize[rank] = 0;
+    for (size_t j = 0; j < routingList[rank].size(); j++)
+    {
+      Elt *elt = (Elt *) (routingList[rank][j].data);
+      sendMessageSize[rank] += packedPolygonSize(*elt);
+    }
+  }
+
+  ep_lib::MPI_Alltoall(nbSendNode, 1, EP_INT, nbRecvNode, 1, EP_INT, communicator);
+  ep_lib::MPI_Alltoall(sendMessageSize, 1, EP_INT, recvMessageSize, 1, EP_INT, communicator);
+
+  char **sendBuffer = new char*[mpiSize];
+  char **recvBuffer = new char*[mpiSize];
+  int *pos = new int[mpiSize];
+
+  for (int rank = 0; rank < mpiSize; rank++)
+  {
+    if (nbSendNode[rank] > 0) sendBuffer[rank] = new char[sendMessageSize[rank]];
+    if (nbRecvNode[rank] > 0) recvBuffer[rank] = new char[recvMessageSize[rank]];
+  }
+
+  for (int rank = 0; rank < mpiSize; rank++)
+  {
+    pos[rank] = 0;
+    for (size_t j = 0; j < routingList[rank].size(); j++)
+    {
+      Elt *elt = (Elt *) (routingList[rank][j].data);
+      packPolygon(*elt, sendBuffer[rank], pos[rank]);
+    }
+  }
+  delete [] routingList;
+
+
+  int nbSendRequest = 0;
+  int nbRecvRequest = 0;
+  ep_lib::MPI_Request *sendRequest = new ep_lib::MPI_Request[mpiSize];
+  ep_lib::MPI_Request *recvRequest = new ep_lib::MPI_Request[mpiSize];
+  ep_lib::MPI_Status  *status      = new ep_lib::MPI_Status[mpiSize];
+
+  for (int rank = 0; rank < mpiSize; rank++)
+  {
+    if (nbSendNode[rank] > 0)
+    {
+      ep_lib::MPI_Issend(sendBuffer[rank], sendMessageSize[rank], EP_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]);
+      nbSendRequest++;
+    }
+    if (nbRecvNode[rank] > 0)
+    {
+      ep_lib::MPI_Irecv(recvBuffer[rank], recvMessageSize[rank], EP_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]);
+      nbRecvRequest++;
+    }
+  }
+
+  ep_lib::MPI_Waitall(nbRecvRequest, recvRequest, status);
+  ep_lib::MPI_Waitall(nbSendRequest, sendRequest, status);
+
+  for (int rank = 0; rank < mpiSize; rank++)
+    if (nbSendNode[rank] > 0) delete [] sendBuffer[rank];
+  delete [] sendBuffer;
+
+  char **sendBuffer2 = new char*[mpiSize];
+  char **recvBuffer2 = new char*[mpiSize];
+
+  for (int rank = 0; rank < mpiSize; rank++)
+  {
+    nbSendNode[rank] = 0;
+    sendMessageSize[rank] = 0;
+
+    if (nbRecvNode[rank] > 0)
+    {
+      set<NodePtr> neighbourList;
+      pos[rank] = 0;
+      for (int j = 0; j < nbRecvNode[rank]; j++)
+      {
+        Elt elt;
+        unpackPolygon(elt, recvBuffer[rank], pos[rank]);
+        Node node(elt.x, cptRadius(elt), &elt);
+        findNeighbour(sstree.localTree.root, &node, neighbourList);
+      }
+      nbSendNode[rank] = neighbourList.size();
+      for (set<NodePtr>::iterator it = neighbourList.begin(); it != neighbourList.end(); it++)
+      {
+        Elt *elt = (Elt *) ((*it)->data);
+        sendMessageSize[rank] += packedPolygonSize(*elt);
+      }
+
+      sendBuffer2[rank] = new char[sendMessageSize[rank]];
+      pos[rank] = 0;
+
+      for (set<NodePtr>::iterator it = neighbourList.begin(); it != neighbourList.end(); it++)
+      {
+        Elt *elt = (Elt *) ((*it)->data);
+        packPolygon(*elt, sendBuffer2[rank], pos[rank]);
+      }
+    }
+  }
+  for (int rank = 0; rank < mpiSize; rank++)
+    if (nbRecvNode[rank] > 0) delete [] recvBuffer[rank];
+  delete [] recvBuffer;
+
+
+  ep_lib::MPI_Barrier(communicator);
+  ep_lib::MPI_Alltoall(nbSendNode, 1, EP_INT, nbRecvNode, 1, EP_INT, communicator);
+  ep_lib::MPI_Alltoall(sendMessageSize, 1, EP_INT, recvMessageSize, 1, EP_INT, communicator);
+
+  for (int rank = 0; rank < mpiSize; rank++)
+    if (nbRecvNode[rank] > 0) recvBuffer2[rank] = new char[recvMessageSize[rank]];
+
+  nbSendRequest = 0;
+  nbRecvRequest = 0;
+
+  for (int rank = 0; rank < mpiSize; rank++)
+  {
+    if (nbSendNode[rank] > 0)
+    {
+      ep_lib::MPI_Issend(sendBuffer2[rank], sendMessageSize[rank], EP_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]);
+      nbSendRequest++;
+    }
+    if (nbRecvNode[rank] > 0)
+    {
+      ep_lib::MPI_Irecv(recvBuffer2[rank], recvMessageSize[rank], EP_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]);
+      nbRecvRequest++;
+    }
+  }
+
+  ep_lib::MPI_Waitall(nbRecvRequest, recvRequest, status);
+  ep_lib::MPI_Waitall(nbSendRequest, sendRequest, status);
+
+  int nbNeighbourNodes = 0;
+  for (int rank = 0; rank < mpiSize; rank++)
+    nbNeighbourNodes += nbRecvNode[rank];
+
+  neighbourElements = new Elt[nbNeighbourNodes];
+  nbNeighbourElements = nbNeighbourNodes;
+
+  int index = 0;
+  for (int rank = 0; rank < mpiSize; rank++)
+  {
+    pos[rank] = 0;
+    for (int j = 0; j < nbRecvNode[rank]; j++)
+    {
+      unpackPolygon(neighbourElements[index], recvBuffer2[rank], pos[rank]);
+      neighbourElements[index].id.ind = sstree.localTree.leafs.size() + index;
+      index++;
+    }
+  }
+  for (int rank = 0; rank < mpiSize; rank++)
+  {
+    if (nbRecvNode[rank] > 0) delete [] recvBuffer2[rank];
+    if (nbSendNode[rank] > 0) delete [] sendBuffer2[rank];
+  }
+  delete [] recvBuffer2;
+  delete [] sendBuffer2;
+  delete [] sendMessageSize;
+  delete [] recvMessageSize;
+  delete [] nbSendNode;
+  delete [] nbRecvNode;
+  delete [] sendRequest;
+  delete [] recvRequest;
+  delete [] status;
+  delete [] pos;
+
+  /* re-compute on received elements to avoid having to send this information */
+  neighbourNodes.resize(nbNeighbourNodes);
+  setCirclesAndLinks(neighbourElements, neighbourNodes);
+  cptAllEltsGeom(neighbourElements, nbNeighbourNodes, srcGrid.pole);
+
+  /* the local SS tree must include nodes from other cpus if they are potential
            intersector of a local node */
-        sstree.localTree.insertNodes(neighbourNodes);
-
-        /* for every local element,
+  sstree.localTree.insertNodes(neighbourNodes);
+
+  /* for every local element,
            use the SS-tree to find all elements (including neighbourElements)
            who are potential neighbours because their circles intersect,
-           then check all canditates for common edges to build up connectivity information
-        */
-        for (int j = 0; j < sstree.localTree.leafs.size(); j++)
-        {
-                Node& node = sstree.localTree.leafs[j];
-
-                /* find all leafs whoes circles that intersect node's circle and save into node->intersectors */
-                node.search(sstree.localTree.root);
-
-                Elt *elt = (Elt *)(node.data);
-
-                for (int i = 0; i < elt->n; i++) elt->neighbour[i] = NOT_FOUND;
-
-                /* for element `elt` loop through all nodes in the SS-tree
+     then check all canditates for common edges to build up connectivity information
+  */
+  for (int j = 0; j < sstree.localTree.leafs.size(); j++)
+  {
+    Node& node = sstree.localTree.leafs[j];
+
+    /* find all leafs whoes circles that intersect node's circle and save into node->intersectors */
+    node.search(sstree.localTree.root);
+
+    Elt *elt = (Elt *)(node.data);
+
+    for (int i = 0; i < elt->n; i++) elt->neighbour[i] = NOT_FOUND;
+
+    /* for element `elt` loop through all nodes in the SS-tree
                    whoes circles intersect with the circle around `elt` (the SS intersectors)
                    and check if they are neighbours in the sense that the two elements share an edge.
                    If they do, save this information for elt */
-                for (list<NodePtr>::iterator it = (node.intersectors).begin(); it != (node.intersectors).end(); ++it)
-                {
-                        Elt *elt2 = (Elt *)((*it)->data);
-                        set_neighbour(*elt, *elt2);
-                }
+    for (list<NodePtr>::iterator it = (node.intersectors).begin(); it != (node.intersectors).end(); ++it)
+    {
+      Elt *elt2 = (Elt *)((*it)->data);
+      set_neighbour(*elt, *elt2);
+    }
 
 /*
-                for (int i = 0; i < elt->n; i++)
-                {
-                        if (elt->neighbour[i] == NOT_FOUND)
-                                error_exit("neighbour not found");
-                }
+    for (int i = 0; i < elt->n; i++)
+    {
+      if (elt->neighbour[i] == NOT_FOUND)
+        error_exit("neighbour not found");
+    }
 */
-        }
+  }
 }
 
@@ -716 +724 @@
 void Mapper::computeIntersection(Elt *elements, int nbElements)
 {
-        int mpiSize, mpiRank;
-        MPI_Comm_size(communicator, &mpiSize);
-        MPI_Comm_rank(communicator, &mpiRank);
-
-        MPI_Barrier(communicator);
-
-        vector<Node> *routingList = new vector<Node>[mpiSize];
-
-        vector<Node> routeNodes;  routeNodes.reserve(nbElements);
-        for (int j = 0; j < nbElements; j++)
-        {
-                elements[j].id.ind = j;
-                elements[j].id.rank = mpiRank;
-                routeNodes.push_back(Node(elements[j].x, cptRadius(elements[j]), &elements[j]));
-        }
-
-        vector<vector<int> > routes(routeNodes.size());
-        sstree.routeIntersections(routes, routeNodes);
-        for (int i = 0; i < routes.size(); ++i)
-                for (int k = 0; k < routes[i].size(); ++k)
-                        routingList[routes[i][k]].push_back(routeNodes[i]);
-
-        if (verbose >= 2)
-        {
-                cout << " --> rank  " << mpiRank << " nbElements " << nbElements << " : ";
-                for (int rank = 0; rank < mpiSize; rank++)
-                        cout << routingList[rank].size() << "   ";
-                cout << endl;
-        }
-        MPI_Barrier(communicator);
-
-        int *nbSendNode = new int[mpiSize];
-        int *nbRecvNode = new int[mpiSize];
-        int *sentMessageSize = new int[mpiSize];
-        int *recvMessageSize = new int[mpiSize];
-
-        for (int rank = 0; rank < mpiSize; rank++)
-        {
-                nbSendNode[rank] = routingList[rank].size();
-                sentMessageSize[rank] = 0;
-                for (size_t j = 0; j < routingList[rank].size(); j++)
-                {
-                        Elt *elt = (Elt *) (routingList[rank][j].data);
-                        sentMessageSize[rank] += packedPolygonSize(*elt);
-                }
-        }
-
-        MPI_Alltoall(nbSendNode, 1, MPI_INT, nbRecvNode, 1, MPI_INT, communicator);
-        MPI_Alltoall(sentMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator);
-
-        int total = 0;
-
-        for (int rank = 0; rank < mpiSize; rank++)
-        {
-                total = total + nbRecvNode[rank];
-        }
-
-        if (verbose >= 2) cout << "---> rank " << mpiRank << " : compute intersection : total received nodes  " << total << endl;
-        char **sendBuffer = new char*[mpiSize];
-        char **recvBuffer = new char*[mpiSize];
-        int *pos = new int[mpiSize];
-
-        for (int rank = 0; rank < mpiSize; rank++)
-        {
-                if (nbSendNode[rank] > 0) sendBuffer[rank] = new char[sentMessageSize[rank]];
-                if (nbRecvNode[rank] > 0) recvBuffer[rank] = new char[recvMessageSize[rank]];
-        }
-
-        for (int rank = 0; rank < mpiSize; rank++)
-        {
-                pos[rank] = 0;
-                for (size_t j = 0; j < routingList[rank].size(); j++)
-                {
-                        Elt* elt = (Elt *) (routingList[rank][j].data);
-                        packPolygon(*elt, sendBuffer[rank], pos[rank]);
-                }
-        }
-        delete [] routingList;
-
-        int nbSendRequest = 0;
-        int nbRecvRequest = 0;
-        MPI_Request *sendRequest = new MPI_Request[mpiSize];
-        MPI_Request *recvRequest = new MPI_Request[mpiSize];
-        MPI_Status   *status = new MPI_Status[mpiSize];
-
-        for (int rank = 0; rank < mpiSize; rank++)
-        {
-                if (nbSendNode[rank] > 0)
-                {
-                        MPI_Issend(sendBuffer[rank], sentMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]);
-                        nbSendRequest++;
-                }
-                if (nbRecvNode[rank] > 0)
-                {
-                        MPI_Irecv(recvBuffer[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]);
-                        nbRecvRequest++;
-                }
-        }
-
-        MPI_Waitall(nbRecvRequest, recvRequest, status);
-        MPI_Waitall(nbSendRequest, sendRequest, status);
-        char **sendBuffer2 = new char*[mpiSize];
-        char **recvBuffer2 = new char*[mpiSize];
-
-        double tic = cputime();
-        for (int rank = 0; rank < mpiSize; rank++)
-        {
-                sentMessageSize[rank] = 0;
-
-                if (nbRecvNode[rank] > 0)
-                {
-                        Elt *recvElt = new Elt[nbRecvNode[rank]];
-                        pos[rank] = 0;
-                        for (int j = 0; j < nbRecvNode[rank]; j++)
-                        {
-                                unpackPolygon(recvElt[j], recvBuffer[rank], pos[rank]);
-                                cptEltGeom(recvElt[j], tgtGrid.pole);
-                                Node recvNode(recvElt[j].x, cptRadius(recvElt[j]), &recvElt[j]);
-                                recvNode.search(sstree.localTree.root);
-                                /* for a node holding an element of the target, loop throught candidates for intersecting source */
-                                for (list<NodePtr>::iterator it = (recvNode.intersectors).begin(); it != (recvNode.intersectors).end(); ++it)
-                                {
-                                        Elt *elt2 = (Elt *) ((*it)->data);
-                                        /* recvElt is target, elt2 is source */
-//                                      intersect(&recvElt[j], elt2);
-                                        intersect_ym(&recvElt[j], elt2);
-                                }
-
-                                if (recvElt[j].is.size() > 0) sentMessageSize[rank] += packIntersectionSize(recvElt[j]);
-
-                                // here recvNode goes out of scope
-                        }
-
-                        if (sentMessageSize[rank] > 0)
-                        {
-                                sentMessageSize[rank] += sizeof(int);
-                                sendBuffer2[rank] = new char[sentMessageSize[rank]];
-                                *((int *) sendBuffer2[rank]) = 0;
-                                pos[rank] = sizeof(int);
-                                for (int j = 0; j < nbRecvNode[rank]; j++)
-                                {
-                                        packIntersection(recvElt[j], sendBuffer2[rank], pos[rank]);
-                                        //FIXME should be deleted: recvElt[j].delete_intersections(); // intersection areas have been packed to buffer and won't be used any more
-                                }
-                        }
-                        delete [] recvElt;
-
-                }
-        }
-        delete [] pos;
-
-        for (int rank = 0; rank < mpiSize; rank++)
-        {
-                if (nbSendNode[rank] > 0) delete [] sendBuffer[rank];
-                if (nbRecvNode[rank] > 0) delete [] recvBuffer[rank];
-                nbSendNode[rank] = 0;
-        }
-
-        if (verbose >= 2) cout << "Rank " << mpiRank << "  Compute (internal) intersection " << cputime() - tic << " s" << endl;
-        MPI_Alltoall(sentMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator);
-
-        for (int rank = 0; rank < mpiSize; rank++)
-                if (recvMessageSize[rank] > 0)
-                        recvBuffer2[rank] = new char[recvMessageSize[rank]];
-
-        nbSendRequest = 0;
-        nbRecvRequest = 0;
-
-        for (int rank = 0; rank < mpiSize; rank++)
-        {
-                if (sentMessageSize[rank] > 0)
-                {
-                        MPI_Issend(sendBuffer2[rank], sentMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]);
-                        nbSendRequest++;
-                }
-                if (recvMessageSize[rank] > 0)
-                {
-                        MPI_Irecv(recvBuffer2[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]);
-                        nbRecvRequest++;
-                }
-        }
-
-        MPI_Waitall(nbRecvRequest, recvRequest, status);
-        MPI_Waitall(nbSendRequest, sendRequest, status);
-
-        delete [] sendRequest;
-        delete [] recvRequest;
-        delete [] status;
-        for (int rank = 0; rank < mpiSize; rank++)
-        {
-                if (nbRecvNode[rank] > 0)
-                {
-                        if (sentMessageSize[rank] > 0)
-                                delete [] sendBuffer2[rank];
-                }
-
-                if (recvMessageSize[rank] > 0)
-                {
-                        unpackIntersection(elements, recvBuffer2[rank]);
-                        delete [] recvBuffer2[rank];
-                }
-        }
-        delete [] sendBuffer2;
-        delete [] recvBuffer2;
-        delete [] sendBuffer;
-        delete [] recvBuffer;
-
-        delete [] nbSendNode;
-        delete [] nbRecvNode;
-        delete [] sentMessageSize;
-        delete [] recvMessageSize;
+  int mpiSize, mpiRank;
+  ep_lib::MPI_Comm_size(communicator, &mpiSize);
+  ep_lib::MPI_Comm_rank(communicator, &mpiRank);
+
+  ep_lib::MPI_Barrier(communicator);
+
+  vector<Node> *routingList = new vector<Node>[mpiSize];
+
+  vector<Node> routeNodes;  routeNodes.reserve(nbElements);
+  for (int j = 0; j < nbElements; j++)
+  {
+    elements[j].id.ind = j;
+    elements[j].id.rank = mpiRank;
+    routeNodes.push_back(Node(elements[j].x, cptRadius(elements[j]), &elements[j]));
+  }
+
+  vector<vector<int> > routes(routeNodes.size());
+  sstree.routeIntersections(routes, routeNodes);
+  for (int i = 0; i < routes.size(); ++i)
+    for (int k = 0; k < routes[i].size(); ++k)
+      routingList[routes[i][k]].push_back(routeNodes[i]);
+
+  if (verbose >= 2)
+  {
+    cout << " --> rank  " << mpiRank << " nbElements " << nbElements << " : ";
+    for (int rank = 0; rank < mpiSize; rank++)
+      cout << routingList[rank].size() << "   ";
+    cout << endl;
+  }
+  ep_lib::MPI_Barrier(communicator);
+
+  int *nbSendNode = new int[mpiSize];
+  int *nbRecvNode = new int[mpiSize];
+  int *sentMessageSize = new int[mpiSize];
+  int *recvMessageSize = new int[mpiSize];
+
+  for (int rank = 0; rank < mpiSize; rank++)
+  {
+    nbSendNode[rank] = routingList[rank].size();
+    sentMessageSize[rank] = 0;
+    for (size_t j = 0; j < routingList[rank].size(); j++)
+    {
+      Elt *elt = (Elt *) (routingList[rank][j].data);
+      sentMessageSize[rank] += packedPolygonSize(*elt);
+    }
+  }
+
+  ep_lib::MPI_Alltoall(nbSendNode, 1, EP_INT, nbRecvNode, 1, EP_INT, communicator);
+  ep_lib::MPI_Alltoall(sentMessageSize, 1, EP_INT, recvMessageSize, 1, EP_INT, communicator);
+
+  int total = 0;
+
+  for (int rank = 0; rank < mpiSize; rank++)
+  {
+    total = total + nbRecvNode[rank];
+  }
+
+  if (verbose >= 2) cout << "---> rank " << mpiRank << " : compute intersection : total received nodes  " << total << endl;
+  char **sendBuffer = new char*[mpiSize];
+  char **recvBuffer = new char*[mpiSize];
+  int *pos = new int[mpiSize];
+
+  for (int rank = 0; rank < mpiSize; rank++)
+  {
+    if (nbSendNode[rank] > 0) sendBuffer[rank] = new char[sentMessageSize[rank]];
+    if (nbRecvNode[rank] > 0) recvBuffer[rank] = new char[recvMessageSize[rank]];
+  }
+
+  for (int rank = 0; rank < mpiSize; rank++)
+  {
+    pos[rank] = 0;
+    for (size_t j = 0; j < routingList[rank].size(); j++)
+    {
+      Elt* elt = (Elt *) (routingList[rank][j].data);
+      packPolygon(*elt, sendBuffer[rank], pos[rank]);
+    }
+  }
+  delete [] routingList;
+
+  int nbSendRequest = 0;
+  int nbRecvRequest = 0;
+  ep_lib::MPI_Request *sendRequest = new ep_lib::MPI_Request[mpiSize];
+  ep_lib::MPI_Request *recvRequest = new ep_lib::MPI_Request[mpiSize];
+  ep_lib::MPI_Status   *status = new ep_lib::MPI_Status[mpiSize];
+
+  for (int rank = 0; rank < mpiSize; rank++)
+  {
+    if (nbSendNode[rank] > 0)
+    {
+      ep_lib::MPI_Issend(sendBuffer[rank], sentMessageSize[rank], EP_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]);
+      nbSendRequest++;
+    }
+    if (nbRecvNode[rank] > 0)
+    {
+      ep_lib::MPI_Irecv(recvBuffer[rank], recvMessageSize[rank], EP_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]);
+      nbRecvRequest++;
+    }
+  }
+
+  ep_lib::MPI_Waitall(nbRecvRequest, recvRequest, status);
+  ep_lib::MPI_Waitall(nbSendRequest, sendRequest, status);
+  char **sendBuffer2 = new char*[mpiSize];
+  char **recvBuffer2 = new char*[mpiSize];
+
+  double tic = cputime();
+  for (int rank = 0; rank < mpiSize; rank++)
+  {
+    sentMessageSize[rank] = 0;
+
+    if (nbRecvNode[rank] > 0)
+    {
+      Elt *recvElt = new Elt[nbRecvNode[rank]];
+      pos[rank] = 0;
+      for (int j = 0; j < nbRecvNode[rank]; j++)
+      {
+        unpackPolygon(recvElt[j], recvBuffer[rank], pos[rank]);
+        cptEltGeom(recvElt[j], tgtGrid.pole);
+        Node recvNode(recvElt[j].x, cptRadius(recvElt[j]), &recvElt[j]);
+        recvNode.search(sstree.localTree.root);
+        /* for a node holding an element of the target, loop throught candidates for intersecting source */
+        for (list<NodePtr>::iterator it = (recvNode.intersectors).begin(); it != (recvNode.intersectors).end(); ++it)
+        {
+          Elt *elt2 = (Elt *) ((*it)->data);
+          /* recvElt is target, elt2 is source */
+//          intersect(&recvElt[j], elt2);
+          intersect_ym(&recvElt[j], elt2);
+        }
+
+        if (recvElt[j].is.size() > 0) sentMessageSize[rank] += packIntersectionSize(recvElt[j]);
+
+        // here recvNode goes out of scope
+      }
+
+      if (sentMessageSize[rank] > 0)
+      {
+        sentMessageSize[rank] += sizeof(int);
+        sendBuffer2[rank] = new char[sentMessageSize[rank]];
+        *((int *) sendBuffer2[rank]) = 0;
+        pos[rank] = sizeof(int);
+        for (int j = 0; j < nbRecvNode[rank]; j++)
+        {
+          packIntersection(recvElt[j], sendBuffer2[rank], pos[rank]);
+          //FIXME should be deleted: recvElt[j].delete_intersections(); // intersection areas have been packed to buffer and won't be used any more
+        }
+      }
+      delete [] recvElt;
+
+    }
+  }
+  delete [] pos;
+
+  for (int rank = 0; rank < mpiSize; rank++)
+  {
+    if (nbSendNode[rank] > 0) delete [] sendBuffer[rank];
+    if (nbRecvNode[rank] > 0) delete [] recvBuffer[rank];
+    nbSendNode[rank] = 0;
+  }
+
+  if (verbose >= 2) cout << "Rank " << mpiRank << "  Compute (internal) intersection " << cputime() - tic << " s" << endl;
+  ep_lib::MPI_Alltoall(sentMessageSize, 1, EP_INT, recvMessageSize, 1, EP_INT, communicator);
+
+  for (int rank = 0; rank < mpiSize; rank++)
+    if (recvMessageSize[rank] > 0)
+      recvBuffer2[rank] = new char[recvMessageSize[rank]];
+
+  nbSendRequest = 0;
+  nbRecvRequest = 0;
+
+  for (int rank = 0; rank < mpiSize; rank++)
+  {
+    if (sentMessageSize[rank] > 0)
+    {
+      ep_lib::MPI_Issend(sendBuffer2[rank], sentMessageSize[rank], EP_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]);
+      nbSendRequest++;
+    }
+    if (recvMessageSize[rank] > 0)
+    {
+      ep_lib::MPI_Irecv(recvBuffer2[rank], recvMessageSize[rank], EP_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]);
+      nbRecvRequest++;
+    }
+  }
+
+  ep_lib::MPI_Waitall(nbRecvRequest, recvRequest, status);
+  ep_lib::MPI_Waitall(nbSendRequest, sendRequest, status);
+
+  delete [] sendRequest;
+  delete [] recvRequest;
+  delete [] status;
+  for (int rank = 0; rank < mpiSize; rank++)
+  {
+    if (nbRecvNode[rank] > 0)
+    {
+      if (sentMessageSize[rank] > 0)
+        delete [] sendBuffer2[rank];
+    }
+
+    if (recvMessageSize[rank] > 0)
+    {
+      unpackIntersection(elements, recvBuffer2[rank]);
+      delete [] recvBuffer2[rank];
+    }
+  }
+  delete [] sendBuffer2;
+  delete [] recvBuffer2;
+  delete [] sendBuffer;
+  delete [] recvBuffer;
+
+  delete [] nbSendNode;
+  delete [] nbRecvNode;
+  delete [] sentMessageSize;
+  delete [] recvMessageSize;
 }
 
 Mapper::~Mapper()
 {
-        delete [] remapMatrix;
-        delete [] srcAddress;
-        delete [] srcRank;
-        delete [] dstAddress;
-        if (neighbourElements) delete [] neighbourElements;
-}
-
-}
+  delete [] remapMatrix;
+  delete [] srcAddress;
+  delete [] srcRank;
+  delete [] dstAddress;
+  if (neighbourElements) delete [] neighbourElements;
+}
+
+}

XIOS/dev/branch_openmp/extern/remap/src/mapper.hpp

@@ -18 +18 @@
 {
 public:
-       Mapper(ep_lib::MPI_Comm comm) : communicator(comm), verbose(SILENT), neighbourElements(NULL), sstree(comm) {}
-
+       Mapper(ep_lib::MPI_Comm comm=EP_COMM_WORLD) : communicator(comm), verbose(SILENT), neighbourElements(NULL), sstree(comm) {}
        ~Mapper();
        void setVerbosity(verbosity v) {verbose=v ;}
 
-       void setSourceMesh(const double* boundsLon, const double* boundsLat, int nVertex, int nbCells, const double* pole, const long int* globalId=NULL) ;
-       void setTargetMesh(const double* boundsLon, const double* boundsLat, int nVertex, int nbCells, const double* pole, const long int* globalId=NULL) ;
+       void setSourceMesh(const double* boundsLon, const double* boundsLat, const double* area, int nVertex, int nbCells, const double* pole, const long int* globalId=NULL) ;
+       void setTargetMesh(const double* boundsLon, const double* boundsLat, const double* area, int nVertex, int nbCells, const double* pole, const long int* globalId=NULL) ;
        void setSourceValue(const double* val) ;
        void getTargetValue(double* val) ;

XIOS/dev/branch_openmp/extern/remap/src/meshutil.cpp

@@ -2 +2 @@
 #include "elt.hpp"
 #include "polyg.hpp"
+#include "intersection_ym.hpp"
+#include "earcut.hpp"
+#include <vector>
 
 namespace sphereRemap {
 
 using namespace std;
+
+double computePolygoneArea(Elt& a, const Coord &pole)
+{
+  using N = uint32_t;
+  using Point = array<double, 2>;
+  vector<Point> vect_points;
+  vector< vector<Point> > polyline;
+  
+  vector<Coord> dstPolygon ;
+  createGreatCirclePolygon(a, pole, dstPolygon) ;
+
+  int na=dstPolygon.size() ;
+  Coord *a_gno   = new Coord[na];
+  
+  Coord OC=barycentre(a.vertex,a.n) ;
+  Coord Oz=OC ;
+  Coord Ox=crossprod(Coord(0,0,1),Oz) ;
+// choose Ox not too small to avoid rounding error
+  if (norm(Ox)< 0.1) Ox=crossprod(Coord(0,1,0),Oz) ;
+  Ox=Ox*(1./norm(Ox)) ;
+  Coord Oy=crossprod(Oz,Ox) ;
+  double cos_alpha;
+
+  for(int n=0; n<na;n++)
+  {
+    cos_alpha=scalarprod(OC,dstPolygon[n]) ;
+    a_gno[n].x=scalarprod(dstPolygon[n],Ox)/cos_alpha ;
+    a_gno[n].y=scalarprod(dstPolygon[n],Oy)/cos_alpha ;
+    a_gno[n].z=scalarprod(dstPolygon[n],Oz)/cos_alpha ; // must be equal to 1
+
+    vect_points.push_back( array<double, 2>() );
+    vect_points[n][0] = a_gno[n].x;
+    vect_points[n][1] = a_gno[n].y;
+
+  }
+
+  polyline.push_back(vect_points);
+  vector<N> indices_a_gno = mapbox::earcut<N>(polyline);
+  
+  double area_a_gno=0 ;
+  for(int i=0;i<indices_a_gno.size()/3;++i)
+    {
+      Coord x0 = Ox * polyline[0][indices_a_gno[3*i]][0] + Oy* polyline[0][indices_a_gno[3*i]][1] + Oz ;
+      Coord x1 = Ox * polyline[0][indices_a_gno[3*i+1]][0] + Oy* polyline[0][indices_a_gno[3*i+1]][1] + Oz ;
+      Coord x2 = Ox * polyline[0][indices_a_gno[3*i+2]][0] + Oy* polyline[0][indices_a_gno[3*i+2]][1] + Oz ;
+      area_a_gno+=triarea(x0 * (1./norm(x0)),x1* (1./norm(x1)), x2* (1./norm(x2))) ;
+    }
+
+  vect_points.clear();
+  polyline.clear();
+  indices_a_gno.clear();
+  return area_a_gno ;
+}
+
 
 void cptEltGeom(Elt& elt, const Coord &pole)
@@ -14 +71 @@
   elt.area = airbar(elt.n, elt.vertex, elt.edge, elt.d, pole, gg);
   elt.x = gg;
-}
+// overwrite area computation 
+
+  elt.area =  computePolygoneArea(elt, pole) ;
+}
+
 
 void cptAllEltsGeom(Elt *elt, int N, const Coord &pole)
@@ -133 +194 @@
   {
     for (int i = 0; i < elts[j]->n; i++) 
-      if ( elts[j]->neighbour[i]== NOT_FOUND) neighbours[i]=NULL ;
+      if ( elts[j]->neighbour[i]== NOT_FOUND) neighbours[i]=NULL ; // no neighbour
+      else if (elts[elts[j]->neighbour[i]]->is.size() == 0) neighbours[i]=NULL ; // neighbour has none supermesh cell 
       else  neighbours[i] = elts[elts[j]->neighbour[i]];
 

XIOS/dev/branch_openmp/extern/remap/src/mpi_cascade.cpp

@@ -1 +1 @@
 #include "mpi_cascade.hpp"
 #include <iostream>
-using namespace ep_lib;
 
 namespace sphereRemap {
 
-CMPICascade::CMPICascade(int nodes_per_level, MPI_Comm comm)
+CMPICascade::CMPICascade(int nodes_per_level, ep_lib::MPI_Comm comm)
 {
         int remaining_levels;
-        MPI_Comm intraComm;
+        ep_lib::MPI_Comm intraComm;
         int l = 0; // current level
         do {
@@ -16 +15 @@
                 level[l].p_grp_size = level[l].size/level[l].group_size;
         
-                MPI_Comm_split(comm, level[l].colour(), level[l].key(), &intraComm);
-                MPI_Comm_split(comm, level[l].p_colour(), level[l].p_key(), &(level[l].pg_comm));
+                ep_lib::MPI_Comm_split(comm, level[l].colour(), level[l].key(), &intraComm);
+                ep_lib::MPI_Comm_split(comm, level[l].p_colour(), level[l].p_key(), &(level[l].pg_comm));
                 comm = intraComm;
                 l++;

XIOS/dev/branch_openmp/extern/remap/src/mpi_cascade.hpp

@@ -12 +12 @@
 {
 public:
-  CCascadeLevel(ep_lib::MPI_Comm comm) : comm(comm)
-  {
-    ep_lib::MPI_Comm_size(comm, &size);
-    ep_lib::MPI_Comm_rank(comm, &rank);
-  }
-  int colour() const { return rank % group_size; };
-  int key() const { return p_colour() + rank/(p_grp_size*group_size)*p_grp_size; }
+        CCascadeLevel(ep_lib::MPI_Comm comm) : comm(comm)
+        {
+                ep_lib::MPI_Comm_size(comm, &size);
+                ep_lib::MPI_Comm_rank(comm, &rank);
+        }
+        int colour() const { return rank % group_size; };
+        int key() const { return p_colour() + rank/(p_grp_size*group_size)*p_grp_size; }
 
-  // perpendicular group
-  int p_colour() const { return (rank%group_size + rank/group_size) % p_grp_size; }
-  int p_key() const { return colour() + rank/(p_grp_size*group_size)*group_size; }
+        // perpendicular group
+        int p_colour() const { return (rank%group_size + rank/group_size) % p_grp_size; }
+        int p_key() const { return colour() + rank/(p_grp_size*group_size)*group_size; }
 
-  ep_lib::MPI_Comm comm, pg_comm;
-  int rank;
-  int size;
-  int group_size; // group_size and p_grp_size are interchanged?? FIXME
-  int p_grp_size;
+        ep_lib::MPI_Comm comm, pg_comm;
+        int rank;
+        int size;
+        int group_size; // group_size and p_grp_size are interchanged?? FIXME
+        int p_grp_size;
 };
 
@@ -34 +34 @@
 {
 public:
-  CMPICascade(int nodes_per_level, ep_lib::MPI_Comm comm);
+        // 
+        CMPICascade(int nodes_per_level, ep_lib::MPI_Comm comm);
 
-  int num_levels;
-  std::vector<CCascadeLevel> level;
+        int num_levels;
+        std::vector<CCascadeLevel> level;
 };
 

XIOS/dev/branch_openmp/extern/remap/src/mpi_routing.cpp

@@ -5 +5 @@
 #include "timerRemap.hpp"
 #include <iostream>
-using namespace ep_lib;
 
 namespace sphereRemap {
@@ -11 +10 @@
 const int verbose = 0;
 
-CMPIRouting::CMPIRouting(MPI_Comm comm) : communicator(comm)
-{
-        MPI_Comm_rank(comm, &mpiRank);
-        MPI_Comm_size(comm, &mpiSize);
+CMPIRouting::CMPIRouting(ep_lib::MPI_Comm comm) : communicator(comm)
+{
+        ep_lib::MPI_Comm_rank(comm, &mpiRank);
+        ep_lib::MPI_Comm_size(comm, &mpiSize);
 }
 
@@ -20 +19 @@
     but message lengths are *known* to receiver */
 template <typename T>
-void alltoalls_known(const vector<vector<T> >& send, vector<vector<T> >& recv, const vector<int>& ranks, MPI_Comm communicator)
-{
-        vector<MPI_Request> request(ranks.size() * 2);
-        vector<MPI_Status>  status(ranks.size() * 2);
+void alltoalls_known(const vector<vector<T> >& send, vector<vector<T> >& recv, const vector<int>& ranks, ep_lib::MPI_Comm communicator)
+{
+        vector<ep_lib::MPI_Request> request(ranks.size() * 2);
+        vector<ep_lib::MPI_Status>  status(ranks.size() * 2);
 
         // communicate data
@@ -29 +28 @@
         for (int i = 0; i < ranks.size(); i++)
                 if (recv[i].size())
-                        MPI_Irecv(&recv[i][0], recv[i].size()*sizeof(T), MPI_CHAR, ranks[i], 0, communicator, &request[nbRequest++]);
+                        ep_lib::MPI_Irecv(&recv[i][0], recv[i].size()*sizeof(T), EP_CHAR, ranks[i], 0, communicator, &request[nbRequest++]);
         for (int i = 0; i < ranks.size(); i++)
                 if (send[i].size())
-                        MPI_Isend((void *) &send[i][0], send[i].size()*sizeof(T), MPI_CHAR, ranks[i], 0, communicator, &request[nbRequest++]);
-        MPI_Waitall(nbRequest, &request[0], &status[0]);
+                        ep_lib::MPI_Isend((void *) &send[i][0], send[i].size()*sizeof(T), EP_CHAR, ranks[i], 0, communicator, &request[nbRequest++]);
+        ep_lib::MPI_Waitall(nbRequest, &request[0], &status[0]);
 }
 
@@ -39 +38 @@
     but message lengths are *unknown* to receiver */
 template <typename T>
-void alltoalls_unknown(const vector<vector<T> >& send, vector<vector<T> >& recv, const vector<int>& ranks, MPI_Comm communicator)
-{
-        vector<MPI_Request> request(ranks.size() * 2);
-        vector<MPI_Status>  status(ranks.size() * 2);
+void alltoalls_unknown(const vector<vector<T> >& send, vector<vector<T> >& recv, const vector<int>& ranks, ep_lib::MPI_Comm communicator)
+{
+        vector<ep_lib::MPI_Request> request(ranks.size() * 2);
+        vector<ep_lib::MPI_Status>  status(ranks.size() * 2);
 
         // communicate sizes
@@ -51 +50 @@
                 sendSizes[i] = send[i].size();
         for (int i = 0; i < ranks.size(); i++)
-                MPI_Irecv(&recvSizes[i], 1, MPI_INT, ranks[i], 0, communicator, &request[nbRequest++]);
+                ep_lib::MPI_Irecv(&recvSizes[i], 1, EP_INT, ranks[i], 0, communicator, &request[nbRequest++]);
         for (int i = 0; i < ranks.size(); i++)
-                MPI_Isend(&sendSizes[i], 1, MPI_INT, ranks[i], 0, communicator, &request[nbRequest++]);
-        MPI_Waitall(nbRequest, &request[0], &status[0]);
+                ep_lib::MPI_Isend(&sendSizes[i], 1, EP_INT, ranks[i], 0, communicator, &request[nbRequest++]);
+        ep_lib::MPI_Waitall(nbRequest, &request[0], &status[0]);
 
         // allocate
@@ -119 +118 @@
         CTimer::get("CMPIRouting::init(reduce_scatter)").reset();
         CTimer::get("CMPIRouting::init(reduce_scatter)").resume();
-        MPI_Reduce_scatter(toSend, &nbSource, recvCount, MPI_INT, MPI_SUM, communicator);
+        ep_lib::MPI_Reduce_scatter(toSend, &nbSource, recvCount, EP_INT, EP_SUM, communicator);
         CTimer::get("CMPIRouting::init(reduce_scatter)").suspend();
         CTimer::get("CMPIRouting::init(reduce_scatter)").print();
 
-        MPI_Alloc_mem(nbTarget *sizeof(int), MPI_INFO_NULL, &targetRank);
-        MPI_Alloc_mem(nbSource *sizeof(int), MPI_INFO_NULL, &sourceRank);
+        ep_lib::MPI_Alloc_mem(nbTarget *sizeof(int), EP_INFO_NULL, &targetRank);
+        ep_lib::MPI_Alloc_mem(nbSource *sizeof(int), EP_INFO_NULL, &sourceRank);
 
         targetRankToIndex = new int[mpiSize];
@@ -138 +137 @@
         }
 
-        MPI_Barrier(communicator);
+        ep_lib::MPI_Barrier(communicator);
         CTimer::get("CMPIRouting::init(get_source)").reset();
         CTimer::get("CMPIRouting::init(get_source)").resume();
 
-        MPI_Request *request = new MPI_Request[nbSource + nbTarget];
-        MPI_Status  *status = new MPI_Status[nbSource + nbTarget];
+        ep_lib::MPI_Request *request = new ep_lib::MPI_Request[nbSource + nbTarget];
+        ep_lib::MPI_Status  *status = new ep_lib::MPI_Status[nbSource + nbTarget];
 
         int indexRequest = 0;
@@ -151 +150 @@
         for (int i = 0; i < nbSource; i++)
         {
-                #ifdef _usingEP
-                MPI_Irecv(&sourceRank[i], 1, MPI_INT, -2, 0, communicator, &request[indexRequest++]);
-                #else
-                MPI_Irecv(&sourceRank[i], 1, MPI_INT, MPI_ANY_SOURCE, 0, communicator, &request[indexRequest++]);
+                #ifdef _usingMPI
+                ep_lib::MPI_Irecv(&sourceRank[i], 1, EP_INT, MPI_ANY_SOURCE, 0, communicator, &request[indexRequest]);
+                #elif _usingEP
+                ep_lib::MPI_Irecv(&sourceRank[i], 1, EP_INT, -2, 0, communicator, &request[indexRequest]);
                 #endif
+                indexRequest++;
         }
         MPI_Barrier(communicator);
         for (int i = 0; i < nbTarget; i++)
         {
-                MPI_Isend(&mpiRank, 1, MPI_INT, targetRank[i], 0, communicator, &request[indexRequest++]);
+                ep_lib::MPI_Isend(&mpiRank, 1, EP_INT, targetRank[i], 0, communicator, &request[indexRequest]);
+                indexRequest++;
         }
         MPI_Waitall(indexRequest, request, status);
@@ -173 +174 @@
         for (int i = 0; i < nbSource; i++)
         {
-                #ifdef _usingEP
-                MPI_Irecv(&sourceRank[i], 1, MPI_INT, -2, 0, communicator, &request[indexRequest]);
-                #else
-                MPI_Irecv(&sourceRank[i], 1, MPI_INT, MPI_ANY_SOURCE, 0, communicator, &request[indexRequest]);
+                #ifdef _usingMPI
+                ep_lib::MPI_Irecv(&sourceRank[i], 1, EP_INT, MPI_ANY_SOURCE, 0, communicator, &request[indexRequest]);
+                #elif _usingEP
+                ep_lib::MPI_Irecv(&sourceRank[i], 1, EP_INT, -2, 0, communicator, &request[indexRequest]);
                 #endif
-                indexRequest++;
-        }
-
-        for (int i = 0; i < nbTarget; i++)
-        {
-                MPI_Isend(&mpiRank, 1, MPI_INT, targetRank[i], 0, communicator, &request[indexRequest]);
+                indexRequest++;
+        }
+
+        for (int i = 0; i < nbTarget; i++)
+        {
+                ep_lib::MPI_Isend(&mpiRank, 1, EP_INT, targetRank[i], 0, communicator, &request[indexRequest]);
                 indexRequest++;
         }
@@ -208 +209 @@
         for (int i = 0; i < nbSource; i++)
         {
-                MPI_Irecv(&nbSourceElement[i], 1, MPI_INT, sourceRank[i], 0, communicator, &request[indexRequest]);
+                ep_lib::MPI_Irecv(&nbSourceElement[i], 1, EP_INT, sourceRank[i], 0, communicator, &request[indexRequest]);
                 indexRequest++;
         }
@@ -215 +216 @@
         {
                 totalTargetElement += nbTargetElement[i];
-                MPI_Isend(&nbTargetElement[i], 1, MPI_INT, targetRank[i], 0, communicator, &request[indexRequest]);
+                ep_lib::MPI_Isend(&nbTargetElement[i], 1, EP_INT, targetRank[i], 0, communicator, &request[indexRequest]);
                 indexRequest++;
         }
@@ -283 +284 @@
 
 
-        MPI_Request* request=new MPI_Request[nbSource+nbTarget];
-        MPI_Status*  status=new MPI_Status[nbSource+nbTarget];
+        ep_lib::MPI_Request* request=new ep_lib::MPI_Request[nbSource+nbTarget];
+        ep_lib::MPI_Status*  status=new ep_lib::MPI_Status[nbSource+nbTarget];
         int indexRequest=0;
 
-        MPI_Barrier(communicator);
+        ep_lib::MPI_Barrier(communicator);
         CTimer::get("CMPIRouting::transferToTarget").reset();
         CTimer::get("CMPIRouting::transferToTarget").resume();
@@ -293 +294 @@
         for(int i=0; i<nbSource; i++)
         {
-                MPI_Irecv(sourceBuffer[i],nbSourceElement[i]*sizeof(T),MPI_CHAR, sourceRank[i], 0, communicator, &request[indexRequest]);
+                ep_lib::MPI_Irecv(sourceBuffer[i],nbSourceElement[i]*sizeof(T),EP_CHAR, sourceRank[i], 0, communicator, &request[indexRequest]);
                 indexRequest++;
         }
@@ -299 +300 @@
         for(int i=0;i<nbTarget; i++)
         {
-                MPI_Isend(targetBuffer[i],nbTargetElement[i]*sizeof(T), MPI_CHAR, targetRank[i], 0, communicator, &request[indexRequest]);
-                indexRequest++;
-        }
-
-        MPI_Waitall(indexRequest,request,status);
+                ep_lib::MPI_Isend(targetBuffer[i],nbTargetElement[i]*sizeof(T), EP_CHAR, targetRank[i], 0, communicator, &request[indexRequest]);
+                indexRequest++;
+        }
+
+        ep_lib::MPI_Waitall(indexRequest,request,status);
 
         CTimer::get("CMPIRouting::transferToTarget").suspend();
         CTimer::get("CMPIRouting::transferToTarget").print();
-        MPI_Barrier(communicator);
+        ep_lib::MPI_Barrier(communicator);
 
         // unpack the data
@@ -347 +348 @@
         }
 
-        MPI_Request *request = new MPI_Request[nbSource + nbTarget];
-        MPI_Status  *status = new MPI_Status[nbSource + nbTarget];
+        ep_lib::MPI_Request *request = new ep_lib::MPI_Request[nbSource + nbTarget];
+        ep_lib::MPI_Status  *status = new ep_lib::MPI_Status[nbSource + nbTarget];
         int indexRequest = 0;
 
-        MPI_Barrier(communicator);
+        ep_lib::MPI_Barrier(communicator);
         CTimer::get("CMPIRouting::transferToTarget(messageSize)").reset();
         CTimer::get("CMPIRouting::transferToTarget(messageSize)").resume();
@@ -357 +358 @@
         for(int i=0; i<nbSource; i++)
         {
-                MPI_Irecv(&sourceMessageSize[i],1,MPI_INT, sourceRank[i], 0, communicator, &request[indexRequest]);
+                ep_lib::MPI_Irecv(&sourceMessageSize[i],1,EP_INT, sourceRank[i], 0, communicator, &request[indexRequest]);
                 indexRequest++;
         }
@@ -363 +364 @@
         for(int i=0; i<nbTarget; i++)
         {
-                MPI_Isend(&targetMessageSize[i],1, MPI_INT, targetRank[i], 0, communicator, &request[indexRequest]);
-                indexRequest++;
-        }
-
-        MPI_Waitall(indexRequest,request,status);
-
-        MPI_Barrier(communicator);
+                ep_lib::MPI_Isend(&targetMessageSize[i],1, EP_INT, targetRank[i], 0, communicator, &request[indexRequest]);
+                indexRequest++;
+        }
+
+        ep_lib::MPI_Waitall(indexRequest,request,status);
+
+        ep_lib::MPI_Barrier(communicator);
         CTimer::get("CMPIRouting::transferToTarget(messageSize)").suspend();
         CTimer::get("CMPIRouting::transferToTarget(messageSize)").print();
@@ -402 +403 @@
         for(int i=0; i<nbSource; i++)
         {
-                MPI_Irecv(sourceBuffer[i],sourceMessageSize[i],MPI_CHAR, sourceRank[i], 0, communicator, &request[indexRequest]);
+                ep_lib::MPI_Irecv(sourceBuffer[i],sourceMessageSize[i],EP_CHAR, sourceRank[i], 0, communicator, &request[indexRequest]);
                 indexRequest++;
         }
@@ -408 +409 @@
         for(int i=0;i<nbTarget; i++)
         {
-                MPI_Isend(targetBuffer[i],targetMessageSize[i], MPI_CHAR, targetRank[i], 0, communicator, &request[indexRequest]);
+                ep_lib::MPI_Isend(targetBuffer[i],targetMessageSize[i], EP_CHAR, targetRank[i], 0, communicator, &request[indexRequest]);
                 indexRequest++;
         }
@@ -467 +468 @@
         }
 
-        MPI_Request* request=new MPI_Request[nbSource+nbTarget];
-        MPI_Status*  status=new MPI_Status[nbSource+nbTarget];
+        ep_lib::MPI_Request* request=new ep_lib::MPI_Request[nbSource+nbTarget];
+        ep_lib::MPI_Status*  status=new ep_lib::MPI_Status[nbSource+nbTarget];
         int indexRequest=0;
 
         for(int i=0; i<nbSource; i++)
         {
-                MPI_Isend(sourceBuffer[i],nbSourceElement[i]*sizeof(T),MPI_CHAR, sourceRank[i], 0, communicator, &request[indexRequest]);
+                ep_lib::MPI_Isend(sourceBuffer[i],nbSourceElement[i]*sizeof(T),EP_CHAR, sourceRank[i], 0, communicator, &request[indexRequest]);
                 indexRequest++;
         }
@@ -479 +480 @@
         for(int i=0;i<nbTarget; i++)
         {
-                MPI_Irecv(targetBuffer[i],nbTargetElement[i]*sizeof(T), MPI_CHAR, targetRank[i], 0, communicator, &request[indexRequest]);
-                indexRequest++;
-        }
-
-        MPI_Waitall(indexRequest,request,status);
+                ep_lib::MPI_Irecv(targetBuffer[i],nbTargetElement[i]*sizeof(T), EP_CHAR, targetRank[i], 0, communicator, &request[indexRequest]);
+                indexRequest++;
+        }
+
+        ep_lib::MPI_Waitall(indexRequest,request,status);
 
         // unpack the data
@@ -523 +524 @@
         }
 
-        MPI_Request *request = new MPI_Request[nbSource + nbTarget];
-        MPI_Status  *status = new MPI_Status[nbSource + nbTarget];
+        ep_lib::MPI_Request *request = new ep_lib::MPI_Request[nbSource + nbTarget];
+        ep_lib::MPI_Status  *status = new ep_lib::MPI_Status[nbSource + nbTarget];
         int indexRequest = 0;
         for (int i = 0; i < nbSource; i++)
         {
-                MPI_Isend(&sourceMessageSize[i], 1, MPI_INT, sourceRank[i], 0, communicator, &request[indexRequest]);
-                indexRequest++;
-        }
-        for (int i = 0; i < nbTarget; i++)
-        {
-                MPI_Irecv(&targetMessageSize[i], 1, MPI_INT, targetRank[i], 0, communicator, &request[indexRequest]);
-                indexRequest++;
-        }
-        MPI_Waitall(indexRequest, request, status);
+                ep_lib::MPI_Isend(&sourceMessageSize[i], 1, EP_INT, sourceRank[i], 0, communicator, &request[indexRequest]);
+                indexRequest++;
+        }
+        for (int i = 0; i < nbTarget; i++)
+        {
+                ep_lib::MPI_Irecv(&targetMessageSize[i], 1, EP_INT, targetRank[i], 0, communicator, &request[indexRequest]);
+                indexRequest++;
+        }
+        ep_lib::MPI_Waitall(indexRequest, request, status);
 
         for (int i = 0; i < nbTarget; i++)
@@ -564 +565 @@
         for (int i = 0; i < nbSource; i++)
         {
-                MPI_Isend(sourceBuffer[i], sourceMessageSize[i], MPI_CHAR, sourceRank[i], 0, communicator, &request[indexRequest]);
-                indexRequest++;
-        }
-        for (int i = 0; i < nbTarget; i++)
-        {
-                MPI_Irecv(targetBuffer[i], targetMessageSize[i], MPI_CHAR, targetRank[i], 0, communicator, &request[indexRequest]);
-                indexRequest++;
-        }
-        MPI_Waitall(indexRequest, request, status);
+                ep_lib::MPI_Isend(sourceBuffer[i], sourceMessageSize[i], EP_CHAR, sourceRank[i], 0, communicator, &request[indexRequest]);
+                indexRequest++;
+        }
+        for (int i = 0; i < nbTarget; i++)
+        {
+                ep_lib::MPI_Irecv(targetBuffer[i], targetMessageSize[i], EP_CHAR, targetRank[i], 0, communicator, &request[indexRequest]);
+                indexRequest++;
+        }
+        ep_lib::MPI_Waitall(indexRequest, request, status);
 
         // unpack the data
@@ -612 +613 @@
 
 template void alltoalls_unknown(const std::vector<std::vector<NES> >& send, std::vector<std::vector<NES> >& recv,
-                                const std::vector<int>& ranks, MPI_Comm communicator);
+                                const std::vector<int>& ranks, ep_lib::MPI_Comm communicator);
 
 template void alltoalls_known(const std::vector<std::vector<int> >& send, std::vector<std::vector<int> >& recv,
-                              const std::vector<int>& ranks, MPI_Comm communicator);
-
-}
+                              const std::vector<int>& ranks, ep_lib::MPI_Comm communicator);
+
+}

XIOS/dev/branch_openmp/extern/remap/src/parallel_tree.cpp

@@ -12 +12 @@
 
 #include "parallel_tree.hpp"
-using namespace ep_lib;
 
 namespace sphereRemap {
-
-extern CRemapGrid srcGrid;
-#pragma omp threadprivate(srcGrid)
-
-extern CRemapGrid tgtGrid;
-#pragma omp threadprivate(tgtGrid)
 
 static const int assignLevel = 2;
@@ -121 +114 @@
 }
 
-CParallelTree::CParallelTree(MPI_Comm comm) : communicator(comm), cascade(MIN_NODE_SZ*MIN_NODE_SZ*2, comm)
-{
-  treeCascade.reserve(cascade.num_levels);
-  for (int lev = 0; lev < cascade.num_levels; lev++)
-    treeCascade.push_back(CSampleTree(cascade.level[lev].group_size, assignLevel));
+//CParallelTree::CParallelTree(MPI_Comm comm) : communicator(comm), cascade(MIN_NODE_SZ*MIN_NODE_SZ, comm)
+CParallelTree::CParallelTree(ep_lib::MPI_Comm comm) : communicator(comm), cascade(MAX_NODE_SZ*MAX_NODE_SZ*2, comm)
+{
+        treeCascade.reserve(cascade.num_levels);
+        for (int lev = 0; lev < cascade.num_levels; lev++)
+                treeCascade.push_back(CSampleTree(cascade.level[lev].group_size, assignLevel));
 }
 
@@ -157 +151 @@
 
         int nrecv; // global number of samples  THIS WILL BE THE NUMBER OF LEAFS IN THE SAMPLE TREE
-        MPI_Allreduce(&n, &nrecv, 1, MPI_INT, MPI_SUM, comm.comm); // => size of sample tree does not depend on keepNodes!
-
+        ep_lib::MPI_Allreduce(&n, &nrecv, 1, EP_INT, EP_SUM, comm.comm); // => size of sample tree does not depend on keepNodes!
         double ratio = blocSize / (1.0 * nrecv);
         int nsend = ratio * n + 1; // nsend = n_local_samples / n_global_samples * blocksize + 1 = blocksize/comm.size
@@ -164 +157 @@
 
         int *counts = new int[comm.size];
-        MPI_Allgather(&nsend, 1, MPI_INT, counts, 1, MPI_INT, comm.comm);
+        ep_lib::MPI_Allgather(&nsend, 1, EP_INT, counts, 1, EP_INT, comm.comm);
 
         nrecv = 0;
@@ -190 +183 @@
         /* each process needs the sample elements from all processes */
         double *recvBuffer = new double[nrecv*4];
-        MPI_Allgatherv(sendBuffer, 4 * nsend, MPI_DOUBLE, recvBuffer, counts, displs, MPI_DOUBLE, comm.comm);
+        ep_lib::MPI_Allgatherv(sendBuffer, 4 * nsend, EP_DOUBLE, recvBuffer, counts, displs, EP_DOUBLE, comm.comm);
         delete[] sendBuffer;
         delete[] counts;
@@ -248 +241 @@
          << "   node size : "<<node.size()<<"   bloc size : "<<blocSize<<"  total number of leaf : "<<tree.leafs.size()<<endl ;
 /*
-        MPI_Allreduce(&ok, &allok, 1, MPI_INT, MPI_PROD, communicator);
+        MPI_Allreduce(&ok, &allok, 1, EP_INT, MPI_PROD, communicator);
         if (!allok) {
                 MPI_Finalize();
@@ -254 +247 @@
         }
 */
-    MPI_Abort(MPI_COMM_WORLD,-1) ;
+    ep_lib::MPI_Abort(EP_COMM_WORLD,-1) ;
   }
 /*
@@ -272 +265 @@
 {
         CMPIRouting MPIRoute(communicator);
-        MPI_Barrier(communicator);
+        ep_lib::MPI_Barrier(communicator);
         CTimer::get("buildLocalTree(initRoute)").resume();
         MPIRoute.init(route);
@@ -297 +290 @@
 
         int mpiRank;
-        MPI_Comm_rank(communicator, &mpiRank);
+        ep_lib::MPI_Comm_rank(communicator, &mpiRank);
         localTree.leafs.reserve(nbLocalElements);
         for (int i = 0; i < nbLocalElements; i++)
@@ -323 +316 @@
   nb1=node.size() ; nb2=node2.size() ;
   nb=nb1+nb2 ;
-  MPI_Allreduce(&nb, &nbTot, 1, MPI_LONG, MPI_SUM, communicator) ;
+  ep_lib::MPI_Allreduce(&nb, &nbTot, 1, EP_LONG, EP_SUM, communicator) ;
   int commSize ;
-  MPI_Comm_size(communicator,&commSize) ;
+  ep_lib::MPI_Comm_size(communicator,&commSize) ;
   
         // make multiple of two
@@ -508 +501 @@
         // gather circles on this level of the cascade
         int pg_size;
-        MPI_Comm_size(cascade.level[level].pg_comm, &pg_size);
+        ep_lib::MPI_Comm_size(cascade.level[level].pg_comm, &pg_size);
         vector<Coord> allRootCentres(pg_size);
         vector<double> allRootRadia(pg_size);
-        MPI_Allgather(&rootCentre, 3, MPI_DOUBLE, &allRootCentres[0], 3, MPI_DOUBLE, cascade.level[level].pg_comm);
-        MPI_Allgather(&rootRadius, 1, MPI_DOUBLE, &allRootRadia[0],   1, MPI_DOUBLE, cascade.level[level].pg_comm);
+        ep_lib::MPI_Allgather(&rootCentre, 3, EP_DOUBLE, &allRootCentres[0], 3, EP_DOUBLE, cascade.level[level].pg_comm);
+        ep_lib::MPI_Allgather(&rootRadius, 1, EP_DOUBLE, &allRootRadia[0],   1, EP_DOUBLE, cascade.level[level].pg_comm);
 
         // now allRootsRadia and allRootCentres must be inserted into second levels of us and propagated to root

XIOS/dev/branch_openmp/extern/remap/src/parallel_tree.hpp

@@ -12 +12 @@
 {
 public:
-  CParallelTree(ep_lib::MPI_Comm comm);
-  ~CParallelTree();
+        CParallelTree(ep_lib::MPI_Comm comm);
+        ~CParallelTree();
 
-  void build(vector<Node>& node, vector<Node>& node2);
+        void build(vector<Node>& node, vector<Node>& node2);
 
-  void routeNodes(vector<int>& route, vector<Node>& nodes, int level = 0);
-  void routeIntersections(vector<vector<int> >& route, vector<Node>& nodes, int level = 0);
+        void routeNodes(vector<int>& route, vector<Node>& nodes, int level = 0);
+        void routeIntersections(vector<vector<int> >& route, vector<Node>& nodes, int level = 0);
 
-  int nbLocalElements;
-  Elt* localElements;
+        int nbLocalElements;
+        Elt* localElements;
 
-  CTree localTree;
+        CTree localTree;
 
 private:
-  void updateCirclesForRouting(Coord rootCentre, double rootRadius, int level = 0);
-  void buildSampleTreeCascade(vector<Node>& sampleNodes, int level = 0);
-  void buildLocalTree(const vector<Node>& node, const vector<int>& route);
-  void buildRouteTree();
+        void updateCirclesForRouting(Coord rootCentre, double rootRadius, int level = 0);
+        void buildSampleTreeCascade(vector<Node>& sampleNodes, int level = 0);
+        void buildLocalTree(const vector<Node>& node, const vector<int>& route);
+        void buildRouteTree();
 
-  //CSampleTree sampleTree;
-  vector<CSampleTree> treeCascade; // first for sample tree, then for routing tree
-  CMPICascade cascade;
+        //CSampleTree sampleTree;
+        vector<CSampleTree> treeCascade; // first for sample tree, then for routing tree
+        CMPICascade cascade;
   ep_lib::MPI_Comm communicator ;
 

XIOS/dev/branch_openmp/extern/remap/src/polyg.cpp

@@ -20 +20 @@
 void orient(int N, Coord *vertex, Coord *edge, double *d, const Coord &g)
 {
-        Coord ga = vertex[0] - g;
-        Coord gb = vertex[1] - g;
-        Coord vertical = crossprod(ga, gb);
-        if (N > 2 && scalarprod(g, vertical) < 0)  // (GAxGB).G
-        {
-                for (int i = 0; i < N/2; i++)
-                        swap(vertex[i], vertex[N-1-i]);
-
-                for (int i = 0; i < (N-1)/2; i++)
-                {
-                        swap(edge[N-2-i], edge[i]);
-                        swap(d[i], d[N-2-i]);
-                }
-        }
+  Coord ga = vertex[0] - g;
+  Coord gb = vertex[1] - g;
+  Coord vertical = crossprod(ga, gb);
+  if (N > 2 && scalarprod(g, vertical) < 0)  // (GAxGB).G
+  {
+    for (int i = 0; i < N/2; i++)
+      swap(vertex[i], vertex[N-1-i]);
+
+    for (int i = 0; i < (N-1)/2; i++)
+    {
+      swap(edge[N-2-i], edge[i]);
+      swap(d[i], d[N-2-i]);
+    }
+  }
 }
 
@@ -39 +39 @@
 void normals(Coord *x, int n, Coord *a)
 {
-        for (int i = 0; i<n; i++)
-                a[i] = crossprod(x[(i+1)%n], x[i]);
+  for (int i = 0; i<n; i++)
+    a[i] = crossprod(x[(i+1)%n], x[i]);
 }
 
 Coord barycentre(const Coord *x, int n)
 {
-        if (n == 0) return ORIGIN;
-        Coord bc = ORIGIN;
-        for (int i = 0; i < n; i++)
-                bc = bc + x[i];
-        /* 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));
-
-        return proj(bc);
+  if (n == 0) return ORIGIN;
+  Coord bc = ORIGIN;
+  for (int i = 0; i < n; i++)
+    bc = bc + x[i];
+  /* 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));
+
+  return proj(bc);
 }
 
@@ -62 +62 @@
 static Coord tetrah_side_diff_centre(Coord a, Coord b)
 {
-        Coord n = crossprod(a,b);
+  Coord n = crossprod(a,b);
         double sinc2 = n.x*n.x + n.y*n.y + n.z*n.z;
-        assert(sinc2 < 1.0 + EPS);
+  assert(sinc2 < 1.0 + EPS);
 
         // exact: u = asin(sinc)/sinc - 1; asin(sinc) = geodesic length of arc ab
-        // approx:
+  // approx:
         // double u = sinc2/6. + (3./40.)*sinc2*sinc2;
 
-        // exact  
-        if (sinc2 > 1.0 - EPS) /* if round-off leads to sinc > 1 asin produces NaN */
-                return n * (M_PI_2 - 1);
-        double sinc = sqrt(sinc2);
+  // exact  
+  if (sinc2 > 1.0 - EPS) /* if round-off leads to sinc > 1 asin produces NaN */
+    return n * (M_PI_2 - 1);
+  double sinc = sqrt(sinc2);
         double u = asin(sinc)/sinc - 1;
 
@@ -82 +82 @@
 Coord gc_normalintegral(const Coord *x, int n)
 {
-        Coord m = barycentre(x, n);
-        Coord bc = crossprod(x[n-1]-m, x[0]-m) + tetrah_side_diff_centre(x[n-1], x[0]);
-        for (int i = 1; i < n; i++)
-                bc = bc + crossprod(x[i-1]-m, x[i]-m) + tetrah_side_diff_centre(x[i-1], x[i]);
-        return bc*0.5;
+  Coord m = barycentre(x, n);
+  Coord bc = crossprod(x[n-1]-m, x[0]-m) + tetrah_side_diff_centre(x[n-1], x[0]);
+  for (int i = 1; i < n; i++)
+    bc = bc + crossprod(x[i-1]-m, x[i]-m) + tetrah_side_diff_centre(x[i-1], x[i]);
+  return bc*0.5;
 }
 
 Coord exact_barycentre(const Coord *x, int n)
 {
-        if (n >= 3)
-        {
-                return  proj(gc_normalintegral(x, n));
-        }
-        else if (n == 0) return ORIGIN;
-        else if (n == 2) return midpoint(x[0], x[1]);
-        else if (n == 1) return x[0];
+  if (n >= 3)
+  {
+    return  proj(gc_normalintegral(x, n));
+  }
+  else if (n == 0) return ORIGIN;
+  else if (n == 2) return midpoint(x[0], x[1]);
+  else if (n == 1) return x[0];
 }
 
 Coord sc_gc_moon_normalintegral(Coord a, Coord b, Coord pole)
 {
-        double hemisphere = (a.z > 0) ? 1: -1;
-
-        double lat = hemisphere * (M_PI_2 - acos(a.z));
-        double lon1 = atan2(a.y, a.x);
-        double lon2 = atan2(b.y, b.x);
-        double lon_diff = lon2 - lon1;
-
-        // wraparound at lon=-pi=pi
-        if (lon_diff < -M_PI) lon_diff += 2.0*M_PI;
-        else if (lon_diff > M_PI) lon_diff -= 2.0*M_PI;
-
-        // integral of the normal over the surface bound by great arcs a-pole and b-pole and small arc a-b 
-        Coord sc_normalintegral = Coord(0.5*(sin(lon2)-sin(lon1))*(M_PI_2 - lat - 0.5*sin(2.0*lat)),
-                                        0.5*(cos(lon1)-cos(lon2))*(M_PI_2 - lat - 0.5*sin(2.0*lat)),
-                                        hemisphere * lon_diff * 0.25 * (cos(2.0*lat) + 1.0));
-        Coord p = Coord(0,0,hemisphere); // TODO assumes north pole is (0,0,1)
-        Coord t[] = {a, b, p};
-        if (hemisphere < 0) swap(t[0], t[1]);
-        return (sc_normalintegral - gc_normalintegral(t, 3)) * hemisphere;
-}
-
-
-double triarea(Coord& A, Coord& B, Coord& C)
-{
-        double a = ds(B, C);
-        double b = ds(C, A);
-        double c = ds(A, B);
-        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<1e-20) return 0. ;
-        return 4 * atan(sqrt(t));
+  double hemisphere = (a.z > 0) ? 1: -1;
+
+  double lat = hemisphere * (M_PI_2 - acos(a.z));
+  double lon1 = atan2(a.y, a.x);
+  double lon2 = atan2(b.y, b.x);
+  double lon_diff = lon2 - lon1;
+
+  // wraparound at lon=-pi=pi
+  if (lon_diff < -M_PI) lon_diff += 2.0*M_PI;
+  else if (lon_diff > M_PI) lon_diff -= 2.0*M_PI;
+
+  // integral of the normal over the surface bound by great arcs a-pole and b-pole and small arc a-b 
+  Coord sc_normalintegral = Coord(0.5*(sin(lon2)-sin(lon1))*(M_PI_2 - lat - 0.5*sin(2.0*lat)),
+                                  0.5*(cos(lon1)-cos(lon2))*(M_PI_2 - lat - 0.5*sin(2.0*lat)),
+                                  hemisphere * lon_diff * 0.25 * (cos(2.0*lat) + 1.0));
+  Coord p = Coord(0,0,hemisphere); // TODO assumes north pole is (0,0,1)
+  Coord t[] = {a, b, p};
+  if (hemisphere < 0) swap(t[0], t[1]);
+  return (sc_normalintegral - gc_normalintegral(t, 3)) * hemisphere;
+}
+
+
+double triarea(const Coord& A, const Coord& B, const Coord& C)
+{
+  double a = ds(B, C);
+  double b = ds(C, A);
+  double c = ds(A, B);
+  double tmp ;
+
+  if (a<b) { tmp=a ; a=b ; b=tmp ; }
+  if (c > a) { tmp=a ; a=c ; c=b, b=tmp;  }
+  else if (c > b) { tmp=c ; c=b ; b=tmp ; }
+  
+  double s = 0.5 * (a + b + c);
+  double t = tan(0.25*(a+(b+c))) * tan(0.25*(c-(a-b))) * tan(0.25*( c + (a-b) )) * tan(0.25*( a + (b - c)));
+  if (t>0) return 4 * atan(sqrt(t));
+  else
+  {
+    std::cout<<"double triarea(const Coord& A, const Coord& B, const Coord& C) : t < 0 !!! t="<<t<<endl ;
+    return 0 ;
+  }
 }
 
@@ -142 +151 @@
 double alun(double b, double d)
 {
-        double a  = acos(d);
-        assert(b <= 2 * a);
-        double s  = a + 0.5 * b;
-        double t = tan(0.5 * s) * tan(0.5 * (s - a)) * tan(0.5 * (s - a)) * tan(0.5 * (s - b));
-        double r  = sqrt(1 - d*d);
-        double p  = 2 * asin(sin(0.5*b) / r);
-        return p*(1 - d) - 4*atan(sqrt(t));
+  double a  = acos(d);
+  assert(b <= 2 * a);
+  double s  = a + 0.5 * b;
+  double t = tan(0.5 * s) * tan(0.5 * (s - a)) * tan(0.5 * (s - a)) * tan(0.5 * (s - b));
+  double r  = sqrt(1 - d*d);
+  double p  = 2 * asin(sin(0.5*b) / r);
+  return p*(1 - d) - 4*atan(sqrt(t));
 }
 
@@ -160 +169 @@
 double airbar(int N, const Coord *x, const Coord *c, double *d, const Coord& pole, Coord& gg)
 {
-        if (N < 3)
-          return 0; /* polygons with less than three vertices have zero area */
-        Coord t[3];
-        t[0] = barycentre(x, N);
-        Coord *g = new Coord[N];
-        double area = 0;
-        Coord gg_exact = gc_normalintegral(x, N);
-        for (int i = 0; i < N; i++)
-        {
-                /* for "spherical circle segment" sum triangular part and "small moon" and => account for small circle */
-                int ii = (i + 1) % N;
-                t[1] = x[i];
-                t[2] = x[ii];
+  if (N < 3)
+    return 0; /* polygons with less then three vertices have zero area */
+  Coord t[3];
+  t[0] = barycentre(x, N);
+  Coord *g = new Coord[N];
+  double area = 0;
+  Coord gg_exact = gc_normalintegral(x, N);
+  for (int i = 0; i < N; i++)
+  {
+    /* for "spherical circle segment" sum triangular part and "small moon" and => account for small circle */
+    int ii = (i + 1) % N;
+    t[1] = x[i];
+    t[2] = x[ii];
                 double sc=scalarprod(crossprod(t[1] - t[0], t[2] - t[0]), t[0]) ;
-                assert(sc >= -1e-10); // Error: tri a l'env (wrong orientation)
-                double area_gc = triarea(t[0], t[1], t[2]);
-                double area_sc_gc_moon = 0;
-                if (d[i]) /* handle small circle case */
-                {
-                        Coord m = midpoint(t[1], t[2]);
-                        double mext = scalarprod(m, c[i]) - d[i];
-                        char sgl = (mext > 0) ? -1 : 1;
-                        area_sc_gc_moon = sgl * alun(arcdist(t[1], t[2]), fabs(scalarprod(t[1], pole)));
-                        gg_exact = gg_exact + sc_gc_moon_normalintegral(t[1], t[2], pole);
-                }
-                area += area_gc + area_sc_gc_moon; /* for "spherical circle segment" sum triangular part (at) and "small moon" and => account for small circle */
-                g[i] = barycentre(t, 3) * (area_gc + area_sc_gc_moon);
-        }
-        gg = barycentre(g, N);
-        gg_exact = proj(gg_exact);
-        delete[] g;
-        gg = gg_exact;
-        return area;
+    assert(sc >= -1e-10); // Error: tri a l'env (wrong orientation)
+    double area_gc = triarea(t[0], t[1], t[2]);
+    double area_sc_gc_moon = 0;
+    if (d[i]) /* handle small circle case */
+    {
+      Coord m = midpoint(t[1], t[2]);
+      double mext = scalarprod(m, c[i]) - d[i];
+      char sgl = (mext > 0) ? -1 : 1;
+      area_sc_gc_moon = sgl * alun(arcdist(t[1], t[2]), fabs(scalarprod(t[1], pole)));
+      gg_exact = gg_exact + sc_gc_moon_normalintegral(t[1], t[2], pole);
+    }
+    area += area_gc + area_sc_gc_moon; /* for "spherical circle segment" sum triangular part (at) and "small moon" and => account for small circle */
+    g[i] = barycentre(t, 3) * (area_gc + area_sc_gc_moon);
+  }
+  gg = barycentre(g, N);
+  gg_exact = proj(gg_exact);
+  delete[] g;
+  gg = gg_exact;
+  return area;
 }
 
 double polygonarea(Coord *vertices, int N)
 {
-        assert(N >= 3);
-
-        /* compute polygon area as sum of triangles */
-        Coord centre = barycentre(vertices, N);
-        double area = 0;
-        for (int i = 0; i < N; i++)
-                area += triarea(centre, vertices[i], vertices[(i+1)%N]);
-        return area;
+  assert(N >= 3);
+
+  /* compute polygon area as sum of triangles */
+  Coord centre = barycentre(vertices, N);
+  double area = 0;
+  for (int i = 0; i < N; i++)
+    area += triarea(centre, vertices[i], vertices[(i+1)%N]);
+  return area;
 }
 
 int packedPolygonSize(const Elt& e)
 {
-        return sizeof(e.id) + sizeof(e.src_id) + sizeof(e.x) + sizeof(e.val) +
-               sizeof(e.n) + e.n*(sizeof(double)+sizeof(Coord));
+  return sizeof(e.id) + sizeof(e.src_id) + sizeof(e.x) + sizeof(e.val)  + sizeof(e.given_area)+
+         sizeof(e.n) + e.n*(sizeof(double)+sizeof(Coord));
 }
 
 void packPolygon(const Elt& e, char *buffer, int& pos) 
 {
-        *((GloId *) &(buffer[pos])) = e.id;
-        pos += sizeof(e.id);
-        *((GloId *) &(buffer[pos])) = e.src_id;
-        pos += sizeof(e.src_id);
-
-        *((Coord *) &(buffer[pos])) = e.x;
-        pos += sizeof(e.x);
-
-        *((double*) &(buffer[pos])) = e.val;
-        pos += sizeof(e.val);
-
-        *((int *) & (buffer[pos])) = e.n;
-        pos += sizeof(e.n);
-
-        for (int i = 0; i < e.n; i++)
-        {
-                *((double *) & (buffer[pos])) = e.d[i];
-                pos += sizeof(e.d[i]);
-
-                *((Coord *) &(buffer[pos])) = e.vertex[i];
-                pos += sizeof(e.vertex[i]);
-        } 
+  *((GloId *) &(buffer[pos])) = e.id;
+  pos += sizeof(e.id);
+  *((GloId *) &(buffer[pos])) = e.src_id;
+  pos += sizeof(e.src_id);
+
+  *((Coord *) &(buffer[pos])) = e.x;
+  pos += sizeof(e.x);
+
+  *((double*) &(buffer[pos])) = e.val;
+  pos += sizeof(e.val);
+
+  *((double*) &(buffer[pos])) = e.given_area;
+  pos += sizeof(e.val);
+
+  *((int *) & (buffer[pos])) = e.n;
+  pos += sizeof(e.n);
+
+  for (int i = 0; i < e.n; i++)
+  {
+    *((double *) & (buffer[pos])) = e.d[i];
+    pos += sizeof(e.d[i]);
+
+    *((Coord *) &(buffer[pos])) = e.vertex[i];
+    pos += sizeof(e.vertex[i]);
+  } 
 
 }
@@ -242 +254 @@
 void unpackPolygon(Elt& e, const char *buffer, int& pos) 
 {
-        e.id = *((GloId *) &(buffer[pos]));
-        pos += sizeof(e.id);
-        e.src_id = *((GloId *) &(buffer[pos]));
-        pos += sizeof(e.src_id);
-
-        e.x = *((Coord *) & (buffer[pos]) );
-        pos += sizeof(e.x);
-
-        e.val = *((double *) & (buffer[pos]));
-        pos += sizeof(double);
-
-        e.n = *((int *) & (buffer[pos]));
-        pos += sizeof(int);
-
-        for (int i = 0; i < e.n; i++)
-        {
-                e.d[i] = *((double *) & (buffer[pos]));
-                pos += sizeof(double);
-
-                e.vertex[i] = *((Coord *) & (buffer[pos]));
-                pos += sizeof(Coord);
-        }
+  e.id = *((GloId *) &(buffer[pos]));
+  pos += sizeof(e.id);
+  e.src_id = *((GloId *) &(buffer[pos]));
+  pos += sizeof(e.src_id);
+
+  e.x = *((Coord *) & (buffer[pos]) );
+  pos += sizeof(e.x);
+
+  e.val = *((double *) & (buffer[pos]));
+  pos += sizeof(double);
+
+  e.given_area = *((double *) & (buffer[pos]));
+  pos += sizeof(double);
+
+  e.n = *((int *) & (buffer[pos]));
+  pos += sizeof(int);
+
+  for (int i = 0; i < e.n; i++)
+  {
+    e.d[i] = *((double *) & (buffer[pos]));
+    pos += sizeof(double);
+
+    e.vertex[i] = *((Coord *) & (buffer[pos]));
+    pos += sizeof(Coord);
+  }
 }
 
 int packIntersectionSize(const Elt& elt) 
 {
-        return elt.is.size() * (2*sizeof(int)+ sizeof(GloId) + 5*sizeof(double));
+  return elt.is.size() * (2*sizeof(int)+ sizeof(GloId) + 5*sizeof(double));
 }
 
@@ -274 +289 @@
 {
   for (list<Polyg *>::const_iterator it = e.is.begin(); it != e.is.end(); ++it)
-        {
-                *((int *) &(buffer[0])) += 1;
-
-                *((int *) &(buffer[pos])) = e.id.ind;
-                pos += sizeof(int);
+  {
+    *((int *) &(buffer[0])) += 1;
+
+    *((int *) &(buffer[pos])) = e.id.ind;
+    pos += sizeof(int);
 
     *((double *) &(buffer[pos])) = e.area;
     pos += sizeof(double);
 
-                *((GloId *) &(buffer[pos])) = (*it)->id;
-                pos += sizeof(GloId);
+
+    *((GloId *) &(buffer[pos])) = (*it)->id;
+    pos += sizeof(GloId);
   
-                *((int *) &(buffer[pos])) = (*it)->n;
-                pos += sizeof(int);
-                *((double *) &(buffer[pos])) = (*it)->area;
-                pos += sizeof(double);
-
-                *((Coord *) &(buffer[pos])) = (*it)->x;
-                pos += sizeof(Coord);
-        }
+    *((int *) &(buffer[pos])) = (*it)->n;
+    pos += sizeof(int);
+    *((double *) &(buffer[pos])) = (*it)->area;
+    pos += sizeof(double);
+
+    *((Coord *) &(buffer[pos])) = (*it)->x;
+    pos += sizeof(Coord);
+  }
 }
 
 void unpackIntersection(Elt* e, const char* buffer) 
 {
-        int ind;
-        int pos = 0;
+  int ind;
+  int pos = 0;
   
-        int n = *((int *) & (buffer[pos]));
-        pos += sizeof(int);
-        for (int i = 0; i < n; i++)
-        {
-                ind = *((int*) &(buffer[pos]));
-                pos+=sizeof(int);
-
-                Elt& elt= e[ind];
+  int n = *((int *) & (buffer[pos]));
+  pos += sizeof(int);
+  for (int i = 0; i < n; i++)
+  {
+    ind = *((int*) &(buffer[pos]));
+    pos+=sizeof(int);
+
+    Elt& elt= e[ind];
 
     elt.area=*((double *) & (buffer[pos]));
-                pos += sizeof(double);
-
-                Polyg *polygon = new Polyg;
-
-                polygon->id =  *((GloId *) & (buffer[pos]));
-                pos += sizeof(GloId);
-
-                polygon->n =  *((int *) & (buffer[pos]));
-                pos += sizeof(int);
-
-                polygon->area =  *((double *) & (buffer[pos]));
-                pos += sizeof(double);
-
-                polygon->x = *((Coord *) & (buffer[pos]));
-                pos += sizeof(Coord);
-
-                elt.is.push_back(polygon);
-        }
-}
-
-}
+    pos += sizeof(double);
+
+
+    Polyg *polygon = new Polyg;
+
+    polygon->id =  *((GloId *) & (buffer[pos]));
+    pos += sizeof(GloId);
+
+    polygon->n =  *((int *) & (buffer[pos]));
+    pos += sizeof(int);
+
+    polygon->area =  *((double *) & (buffer[pos]));
+    pos += sizeof(double);
+
+    polygon->x = *((Coord *) & (buffer[pos]));
+    pos += sizeof(Coord);
+
+    elt.is.push_back(polygon);
+  }
+}
+
+}

XIOS/dev/branch_openmp/extern/remap/src/polyg.hpp

@@ -23 +23 @@
 void unpackIntersection(Elt *e, const char *buffer);
 int packIntersectionSize(const Elt& e);
+double triarea( const Coord& A, const Coord& B, const Coord& C) ;
 
 }

XIOS/dev/branch_openmp/extern/remap/src/timerRemap.cpp

@@ -4 +4 @@
 #include <map>
 #include <iostream>
-using namespace ep_lib;
 
 namespace sphereRemap {
@@ -10 +9 @@
 using namespace std;
 
-//map<string,CTimer*> CTimer::allTimer;
-map<string,CTimer*> *CTimer::allTimer_ptr = 0;
+map<string,CTimer*> CTimer::allTimer;
 
 CTimer::CTimer(const string& name_) : name(name_)
@@ -58 +56 @@
 {
         map<string,CTimer*>::iterator it;
-        if(allTimer_ptr == 0) allTimer_ptr = new map<string,CTimer*>;
-        //it=allTimer.find(name);
-        it=allTimer_ptr->find(name);
-        //if (it==allTimer.end()) it=allTimer.insert(pair<string,CTimer*>(name,new CTimer(name))).first;
-        if (it==allTimer_ptr->end()) it=allTimer_ptr->insert(pair<string,CTimer*>(name,new CTimer(name))).first;
+        it=allTimer.find(name);
+        if (it==allTimer.end()) it=allTimer.insert(pair<string,CTimer*>(name,new CTimer(name))).first;
         return *(it->second);
 }

XIOS/dev/branch_openmp/extern/remap/src/timerRemap.hpp

@@ -26 +26 @@
     double getCumulatedTime(void);
     void print(void);
-    //static map<string,CTimer*> allTimer;
-    static map<string,CTimer*> *allTimer_ptr;
-    #pragma omp threadprivate(allTimer_ptr)
+    static map<string,CTimer*> allTimer;
     static double getTime(void);
     static CTimer& get(string name);

XIOS/dev/branch_openmp/extern/src_ep_dev/ep_allgather.cpp

@@ -92 +92 @@
       int local_sendcount = num_ep * count;
 
-      ::MPI_Allgather(&local_sendcount, 1, to_mpi_type(MPI_INT), mpi_recvcounts, 1, to_mpi_type(MPI_INT), to_mpi_comm(comm->mpi_comm));
+      ::MPI_Allgather(&local_sendcount, 1, to_mpi_type(EP_INT), mpi_recvcounts, 1, to_mpi_type(EP_INT), to_mpi_comm(comm->mpi_comm));
 
       mpi_displs[0] = 0;

XIOS/dev/branch_openmp/extern/src_ep_dev/ep_allgatherv.cpp

@@ -56 +56 @@
     vector<int>local_displs(num_ep, 0);
 
-    MPI_Gather_local(&sendcount, 1, MPI_INT, local_recvcounts.data(), 0, comm);
+    MPI_Gather_local(&sendcount, 1, EP_INT, local_recvcounts.data(), 0, comm);
     for(int i=1; i<num_ep; i++) local_displs[i] = local_displs[i-1] + local_recvcounts[i-1]; 
 
@@ -75 +75 @@
 
       int local_sendcount = std::accumulate(local_recvcounts.begin(), local_recvcounts.end(), 0);
-      ::MPI_Allgather(&local_sendcount, 1, to_mpi_type(MPI_INT), mpi_recvcounts.data(), 1, to_mpi_type(MPI_INT), to_mpi_comm(comm->mpi_comm));
+      ::MPI_Allgather(&local_sendcount, 1, to_mpi_type(EP_INT), mpi_recvcounts.data(), 1, to_mpi_type(EP_INT), to_mpi_comm(comm->mpi_comm));
 
       for(int i=1; i<mpi_size; i++)

XIOS/dev/branch_openmp/extern/src_ep_dev/ep_allocate.cpp

@@ -33 +33 @@
 
     int num_ep_max;
-    MPI_Allreduce(&num_ep, &num_ep_max, 1, MPI_INT, MPI_MAX, comm);
+    MPI_Allreduce(&num_ep, &num_ep_max, 1, EP_INT, MPI_MAX, comm);
 
     assert(num_ep_max > 1);

XIOS/dev/branch_openmp/extern/src_ep_dev/ep_create.cpp

@@ -73 +73 @@
     }
 
-    ::MPI_Allgather(&num_ep, 1, to_mpi_type(MPI_INT), &recv_num_ep[0], 1, to_mpi_type(MPI_INT), mpi_base_comm);
+    ::MPI_Allgather(&num_ep, 1, to_mpi_type(EP_INT), &recv_num_ep[0], 1, to_mpi_type(EP_INT), mpi_base_comm);
 
 
@@ -124 +124 @@
       for(int j=0; j<recv_num_ep[i]; j++)
       {
-        out_comm_hdls[0]->ep_rank_map->insert(std::pair< int, std::pair<int,int> >(ind, j, i));
+        //out_comm_hdls[0]->ep_rank_map->insert(std::pair< int, std::pair<int,int> >(ind, j, i));
+        (*(out_comm_hdls[0]->ep_rank_map))[ind] = std::make_pair(j,i);//->insert(std::pair< int, std::pair<int,int> >(ind, j, i));
         ind++;
       }

XIOS/dev/branch_openmp/extern/src_ep_dev/ep_declaration.cpp

@@ -4 +4 @@
 
 ::MPI_Comm MPI_COMM_WORLD_STD = MPI_COMM_WORLD;
-#undef MPI_COMM_WORLD
+//#undef MPI_COMM_WORLD
 
 
 ::MPI_Comm MPI_COMM_NULL_STD = MPI_COMM_NULL;
-#undef MPI_COMM_NULL
+//#undef MPI_COMM_NULL
 
 
 ::MPI_Request MPI_REQUEST_NULL_STD = MPI_REQUEST_NULL;
-#undef MPI_REQUEST_NULL
+//#undef MPI_REQUEST_NULL
 
 ::MPI_Info MPI_INFO_NULL_STD = MPI_INFO_NULL;
-#undef MPI_INFO_NULL
+//#undef MPI_INFO_NULL
 
 ::MPI_Datatype MPI_INT_STD = MPI_INT;
@@ -26 +26 @@
 ::MPI_Datatype MPI_UINT64_T_STD = MPI_UINT64_T;
 ::MPI_Datatype MPI_LONG_LONG_INT_STD = MPI_LONG_LONG_INT;
-
+::MPI_Datatype MPI_LONG_LONG_STD = MPI_LONG_LONG;
+/*
 #undef MPI_INT
 #undef MPI_FLOAT
@@ -36 +37 @@
 #undef MPI_UINT64_T
 #undef MPI_LONG_LONG_INT
-
+*/
 
 ::MPI_Op MPI_SUM_STD = MPI_SUM;
@@ -43 +44 @@
 ::MPI_Op MPI_LOR_STD = MPI_LOR;
 ::MPI_Op MPI_REPLACE_STD = MPI_REPLACE;
-
+/*
 #undef MPI_SUM
 #undef MPI_MAX
@@ -49 +50 @@
 #undef MPI_LOR
 #undef MPI_REPLACE
-
+*/
 
 // _STD defined in ep_type.cpp
@@ -62 +63 @@
 extern ::MPI_Datatype MPI_UINT64_T_STD;
 extern ::MPI_Datatype MPI_LONG_LONG_INT_STD;
+extern ::MPI_Datatype MPI_LONG_LONG_STD;
+
 
 
@@ -77 +80 @@
 extern ::MPI_Info MPI_INFO_NULL_STD;
 
-ep_lib::MPI_Datatype MPI_INT = &MPI_INT_STD;
-ep_lib::MPI_Datatype MPI_FLOAT = &MPI_FLOAT_STD;
-ep_lib::MPI_Datatype MPI_DOUBLE = &MPI_DOUBLE_STD;
-ep_lib::MPI_Datatype MPI_CHAR = &MPI_CHAR_STD;
-ep_lib::MPI_Datatype MPI_LONG = &MPI_LONG_STD;
-ep_lib::MPI_Datatype MPI_UNSIGNED_LONG = &MPI_UNSIGNED_LONG_STD;
-ep_lib::MPI_Datatype MPI_UNSIGNED_CHAR = &MPI_UNSIGNED_CHAR_STD;
-ep_lib::MPI_Datatype MPI_UINT64_T = &MPI_UINT64_T_STD;
-ep_lib::MPI_Datatype MPI_LONG_LONG_INT = &MPI_LONG_LONG_INT_STD;
-
-
-ep_lib::MPI_Op MPI_SUM = &MPI_SUM_STD;
-ep_lib::MPI_Op MPI_MAX = &MPI_MAX_STD;
-ep_lib::MPI_Op MPI_MIN = &MPI_MIN_STD;
-ep_lib::MPI_Op MPI_LOR = &MPI_LOR_STD;
-ep_lib::MPI_Op MPI_REPLACE = &MPI_REPLACE_STD;
-
-ep_lib::ep_comm EP_COMM_WORLD(&MPI_COMM_WORLD_STD);
-ep_lib::ep_comm EP_COMM_NULL(&MPI_COMM_NULL_STD);
-
-ep_lib::MPI_Comm MPI_COMM_WORLD = &EP_COMM_WORLD;
-ep_lib::MPI_Comm MPI_COMM_NULL = &EP_COMM_NULL;
-
-//ep_lib::ep_status EP_STATUS_IGNORE(&MPI_STATUS_IGNORE_STD);
-ep_lib::ep_request EP_REQUEST_NULL(&MPI_REQUEST_NULL_STD);
-ep_lib::ep_info EP_INFO_NULL(&MPI_INFO_NULL_STD);
-
-//ep_lib::MPI_Status MPI_STATUS_IGNORE = &EP_STATUS_IGNORE;
-ep_lib::MPI_Request MPI_REQUEST_NULL = &EP_REQUEST_NULL;
-ep_lib::MPI_Info MPI_INFO_NULL = &EP_INFO_NULL;
+ep_lib::MPI_Datatype EP_INT = &MPI_INT_STD;
+ep_lib::MPI_Datatype EP_FLOAT = &MPI_FLOAT_STD;
+ep_lib::MPI_Datatype EP_DOUBLE = &MPI_DOUBLE_STD;
+ep_lib::MPI_Datatype EP_CHAR = &MPI_CHAR_STD;
+ep_lib::MPI_Datatype EP_LONG = &MPI_LONG_STD;
+ep_lib::MPI_Datatype EP_UNSIGNED_LONG = &MPI_UNSIGNED_LONG_STD;
+ep_lib::MPI_Datatype EP_UNSIGNED_CHAR = &MPI_UNSIGNED_CHAR_STD;
+ep_lib::MPI_Datatype EP_UINT64_T = &MPI_UINT64_T_STD;
+ep_lib::MPI_Datatype EP_LONG_LONG_INT = &MPI_LONG_LONG_INT_STD;
+ep_lib::MPI_Datatype EP_LONG_LONG = &MPI_LONG_LONG_STD;
 
 
 
+ep_lib::MPI_Op EP_SUM = &MPI_SUM_STD;
+ep_lib::MPI_Op EP_MAX = &MPI_MAX_STD;
+ep_lib::MPI_Op EP_MIN = &MPI_MIN_STD;
+ep_lib::MPI_Op EP_LOR = &MPI_LOR_STD;
+ep_lib::MPI_Op EP_REPLACE = &MPI_REPLACE_STD;
+
+ep_lib::ep_comm EP_COMM_WORLD_t(&MPI_COMM_WORLD_STD);
+ep_lib::ep_comm EP_COMM_NULL_t(&MPI_COMM_NULL_STD);
+
+ep_lib::MPI_Comm EP_COMM_WORLD = &EP_COMM_WORLD_t;
+ep_lib::MPI_Comm EP_COMM_NULL = &EP_COMM_NULL_t;
+
+//ep_lib::ep_status EP_STATUS_IGNORE(&MPI_STATUS_IGNORE_STD);
+ep_lib::ep_request EP_REQUEST_NULL_t(&MPI_REQUEST_NULL_STD);
+ep_lib::ep_info EP_INFO_NULL_t(&MPI_INFO_NULL_STD);
+
+//ep_lib::MPI_Status MPI_STATUS_IGNORE = &EP_STATUS_IGNORE;
+ep_lib::MPI_Request EP_REQUEST_NULL = &EP_REQUEST_NULL_t;
+ep_lib::MPI_Info EP_INFO_NULL = &EP_INFO_NULL_t;
+
+
+

XIOS/dev/branch_openmp/extern/src_ep_dev/ep_declaration.hpp

@@ -1 +1 @@
 #ifndef EP_DECLARATION_HPP_INCLUDED
 #define EP_DECLARATION_HPP_INCLUDED
-
+/*
 #undef MPI_INT
 #undef MPI_FLOAT
@@ -18 +18 @@
 #undef MPI_REPLACE
 
-#undef MPI_COMM_WORLD
+//#undef MPI_COMM_WORLD
 #undef MPI_COMM_NULL
 
@@ -24 +24 @@
 #undef MPI_STATUS_IGNORE
 #undef MPI_INFO_NULL
+*/
+extern ep_lib::MPI_Datatype EP_INT;
+extern ep_lib::MPI_Datatype EP_FLOAT;
+extern ep_lib::MPI_Datatype EP_DOUBLE;
+extern ep_lib::MPI_Datatype EP_CHAR;
+extern ep_lib::MPI_Datatype EP_LONG;
+extern ep_lib::MPI_Datatype EP_UNSIGNED_LONG;
+extern ep_lib::MPI_Datatype EP_UNSIGNED_CHAR;
+extern ep_lib::MPI_Datatype EP_UINT64_T;
+extern ep_lib::MPI_Datatype EP_LONG_LONG_INT;
+extern ep_lib::MPI_Datatype EP_LONG_LONG;
 
-extern ep_lib::MPI_Datatype MPI_INT;
-extern ep_lib::MPI_Datatype MPI_FLOAT;
-extern ep_lib::MPI_Datatype MPI_DOUBLE;
-extern ep_lib::MPI_Datatype MPI_CHAR;
-extern ep_lib::MPI_Datatype MPI_LONG;
-extern ep_lib::MPI_Datatype MPI_UNSIGNED_LONG;
-extern ep_lib::MPI_Datatype MPI_UNSIGNED_CHAR;
-extern ep_lib::MPI_Datatype MPI_UINT64_T;
-extern ep_lib::MPI_Datatype MPI_LONG_LONG_INT;
 
-extern ep_lib::MPI_Op MPI_SUM;
-extern ep_lib::MPI_Op MPI_MAX;
-extern ep_lib::MPI_Op MPI_MIN;
-extern ep_lib::MPI_Op MPI_LOR;
-extern ep_lib::MPI_Op MPI_REPLACE;
+extern ep_lib::MPI_Op EP_SUM;
+extern ep_lib::MPI_Op EP_MAX;
+extern ep_lib::MPI_Op EP_MIN;
+extern ep_lib::MPI_Op EP_LOR;
+extern ep_lib::MPI_Op EP_REPLACE;
 
-extern ep_lib::MPI_Comm MPI_COMM_WORLD;
-extern ep_lib::MPI_Comm MPI_COMM_NULL;
+extern ep_lib::MPI_Comm EP_COMM_WORLD;
+extern ep_lib::MPI_Comm EP_COMM_NULL;
 
-extern ep_lib::MPI_Status MPI_STATUS_IGNORE;
-extern ep_lib::MPI_Request MPI_REQUEST_NULL;
-extern ep_lib::MPI_Info MPI_INFO_NULL;
-
+extern ep_lib::MPI_Status EP_STATUS_IGNORE;
+extern ep_lib::MPI_Request EP_REQUEST_NULL;
+extern ep_lib::MPI_Info EP_INFO_NULL;
 
 #endif // EP_DECLARATION_HPP_INCLUDED

XIOS/dev/branch_openmp/extern/src_ep_dev/ep_exscan.cpp

@@ -78 +78 @@
     MPI_Barrier_local(comm);
 
-    if(op == MPI_SUM)
-    {
-      if(datatype == MPI_INT )
+    if(op == EP_SUM)
+    {
+      if(datatype == EP_INT )
       {
         assert(datasize == sizeof(int));
@@ -87 +87 @@
       }
      
-      else if(datatype == MPI_FLOAT )
+      else if(datatype == EP_FLOAT )
       {
         assert(datasize == sizeof(float));
@@ -95 +95 @@
       
 
-      else if(datatype == MPI_DOUBLE )
+      else if(datatype == EP_DOUBLE )
       {
         assert(datasize == sizeof(double));
@@ -102 +102 @@
       }
 
-      else if(datatype == MPI_CHAR )
+      else if(datatype == EP_CHAR )
       {
         assert(datasize == sizeof(char));
@@ -109 +109 @@
       }
 
-      else if(datatype == MPI_LONG )
+      else if(datatype == EP_LONG )
       {
         assert(datasize == sizeof(long));
@@ -116 +116 @@
       }
 
-      else if(datatype == MPI_UNSIGNED_LONG )
+      else if(datatype == EP_UNSIGNED_LONG )
       {
         assert(datasize == sizeof(unsigned long));
@@ -123 +123 @@
       }
       
-      else if(datatype == MPI_LONG_LONG_INT )
+      else if(datatype == EP_LONG_LONG_INT )
       {
         assert(datasize == sizeof(long long int));
@@ -139 +139 @@
     }
 
-    else if(op == MPI_MAX)
-    {
-      if(datatype == MPI_INT )
+    else if(op == EP_MAX)
+    {
+      if(datatype == EP_INT )
       {
         assert(datasize == sizeof(int));
@@ -148 +148 @@
       }
 
-      else if(datatype == MPI_FLOAT )
+      else if(datatype == EP_FLOAT )
       {
         assert(datasize == sizeof(float));
@@ -155 +155 @@
       }
 
-      else if(datatype == MPI_DOUBLE )
+      else if(datatype == EP_DOUBLE )
       {
         assert(datasize == sizeof(double));
@@ -162 +162 @@
       }
 
-      else if(datatype == MPI_CHAR )
+      else if(datatype == EP_CHAR )
       {
         assert(datasize == sizeof(char));
@@ -169 +169 @@
       }
 
-      else if(datatype == MPI_LONG )
+      else if(datatype == EP_LONG )
       {
         assert(datasize == sizeof(long));
@@ -176 +176 @@
       }
 
-      else if(datatype == MPI_UNSIGNED_LONG )
+      else if(datatype == EP_UNSIGNED_LONG )
       {
         assert(datasize == sizeof(unsigned long));
@@ -183 +183 @@
       }
      
-      else if(datatype == MPI_LONG_LONG_INT )
+      else if(datatype == EP_LONG_LONG_INT )
       {
         assert(datasize == sizeof(long long int));
@@ -197 +197 @@
     }
 
-    else if(op == MPI_MIN)
-    {
-      if(datatype == MPI_INT )
+    else if(op == EP_MIN)
+    {
+      if(datatype == EP_INT )
       {
         assert(datasize == sizeof(int));
@@ -206 +206 @@
       }
 
-      else if(datatype == MPI_FLOAT )
+      else if(datatype == EP_FLOAT )
       {
         assert(datasize == sizeof(float));
@@ -213 +213 @@
       }
 
-      else if(datatype == MPI_DOUBLE )
+      else if(datatype == EP_DOUBLE )
       {
         assert(datasize == sizeof(double));
@@ -220 +220 @@
       }
 
-      else if(datatype == MPI_CHAR )
+      else if(datatype == EP_CHAR )
       {
         assert(datasize == sizeof(char));
@@ -227 +227 @@
       }
 
-      else if(datatype == MPI_LONG )
+      else if(datatype == EP_LONG )
       {
         assert(datasize == sizeof(long));
@@ -234 +234 @@
       }
 
-      else if(datatype == MPI_UNSIGNED_LONG )
+      else if(datatype == EP_UNSIGNED_LONG )
       {
         assert(datasize == sizeof(unsigned long));
@@ -241 +241 @@
       }
 
-      else if(datatype == MPI_LONG_LONG_INT )
+      else if(datatype == EP_LONG_LONG_INT )
       {
         assert(datasize == sizeof(long long int));

XIOS/dev/branch_openmp/extern/src_ep_dev/ep_finalize.cpp

@@ -9 +9 @@
   int MPI_Finalize()
   {
-    #pragma omp master
+    //pragma omp master
     {
       printf("calling EP Finalize\n");

XIOS/dev/branch_openmp/extern/src_ep_dev/ep_fortran.cpp

@@ -8 +8 @@
 
 
-namespace ep_lib
-{
+//namespace ep_lib
+//{
 
-  void* EP_Comm_c2f(MPI_Comm comm)
+  void* EP_Comm_c2f(ep_lib::MPI_Comm comm)
   {
     Debug("MPI_Comm_c2f");
-    void* fint = new ::MPI_Fint;
+    void* fint = new MPI_Fint;
     #ifdef _intelmpi
-    *static_cast< ::MPI_Fint*>(fint) = (::MPI_Fint)(to_mpi_comm(comm->mpi_comm));
+    //*static_cast< ::MPI_Fint*>(fint) = (::MPI_Fint)(to_mpi_comm(comm->mpi_comm));
+    *static_cast< MPI_Fint*>(fint) = MPI_Comm_c2f(to_mpi_comm(comm->mpi_comm));
     #elif _openmpi
-    *static_cast< ::MPI_Fint*>(fint) = MPI_Comm_c2f(to_mpi_comm(comm->mpi_comm));
+    *static_cast< MPI_Fint*>(fint) = MPI_Comm_c2f(to_mpi_comm(comm->mpi_comm));
     #endif
     
-    std::map<std::pair< ::MPI_Fint, int>, MPI_Comm > ::iterator it;
+    std::map<std::pair< MPI_Fint, int>, ep_lib::MPI_Comm > ::iterator it;
     
     #pragma omp critical (fc_comm_map)
     {
-      it = fc_comm_map.find(std::make_pair(*static_cast< ::MPI_Fint*>(fint), omp_get_thread_num()));
-      if(it == fc_comm_map.end())
+      it = ep_lib::fc_comm_map.find(std::make_pair(*static_cast< MPI_Fint*>(fint), omp_get_thread_num()));
+      if(it == ep_lib::fc_comm_map.end())
       {
-        fc_comm_map.insert(std::make_pair( std::make_pair( *static_cast< ::MPI_Fint*>(fint), omp_get_thread_num()) , comm));
-        printf("EP_Comm_c2f : MAP %p insert: %d, %d, %p\n", &fc_comm_map, *static_cast< ::MPI_Fint*>(fint), omp_get_thread_num(), comm->ep_comm_ptr);
+        ep_lib::fc_comm_map.insert(std::make_pair( std::make_pair( *static_cast< MPI_Fint*>(fint), omp_get_thread_num()) , comm));
+        printf("EP_Comm_c2f : MAP %p insert: %d, %d, %p\n", &(ep_lib::fc_comm_map), *static_cast< MPI_Fint*>(fint), omp_get_thread_num(), comm->ep_comm_ptr);
       }
     }
@@ -38 +39 @@
   }
 
-  MPI_Comm EP_Comm_f2c(void* comm)
+  ep_lib::MPI_Comm EP_Comm_f2c(void* comm)
   {
     Debug("MPI_Comm_f2c");
     
     
-    std::map<std::pair< ::MPI_Fint, int>, MPI_Comm > ::iterator it;
+    std::map<std::pair< MPI_Fint, int>, ep_lib::MPI_Comm > ::iterator it;
     
     #pragma omp critical (fc_comm_map)
-    it = fc_comm_map.find(std::make_pair(*static_cast< ::MPI_Fint*>(comm), omp_get_thread_num()));
+    it = ep_lib::fc_comm_map.find(std::make_pair(*static_cast< MPI_Fint*>(comm), omp_get_thread_num()));
     
-    if(it != fc_comm_map.end())
+    if(it != ep_lib::fc_comm_map.end())
     {
-      MPI_Comm comm_ptr;
+      ep_lib::MPI_Comm comm_ptr;
       comm_ptr = it->second;
-      printf("EP_Comm_f2c : MAP %p find: %d, %d, %p\n", &fc_comm_map, it->first.first, it->first.second, comm_ptr->ep_comm_ptr);
+      printf("EP_Comm_f2c : MAP %p find: %d, %d, %p\n", &(ep_lib::fc_comm_map), it->first.first, it->first.second, comm_ptr->ep_comm_ptr);
       return  comm_ptr;
     }
@@ -58 +59 @@
       
     
+    MPI_Comm *base_comm = new MPI_Comm;
     #ifdef _intelmpi
-    ::MPI_Comm *base_comm = new ::MPI_Comm;
-    *base_comm = (::MPI_Comm)(*static_cast< ::MPI_Fint*>(comm));
+    *base_comm = (MPI_Comm)(*static_cast< MPI_Fint*>(comm));
     #elif _openmpi
-    ::MPI_Comm *base_comm = ::MPI_Comm_f2c(*static_cast< ::MPI_Fint*>(comm));
+    *base_comm = MPI_Comm_f2c(*static_cast< MPI_Fint*>(comm));
     #endif
 
-    if(*base_comm != to_mpi_comm(MPI_COMM_NULL->mpi_comm))
+    if(*base_comm != to_mpi_comm(EP_COMM_NULL->mpi_comm))
     {
       if(omp_get_thread_num() == 0)
       {
         int num_ep = omp_get_num_threads();
-        MPI_Comm *new_comm;
-        MPI_Info info;
-        MPI_Comm_create_endpoints(base_comm, num_ep, info, new_comm);
-        passage = new_comm;
+        ep_lib::MPI_Comm *new_comm;
+        ep_lib::MPI_Info info;
+        ep_lib::MPI_Comm_create_endpoints(base_comm, num_ep, info, new_comm);
+        ep_lib::passage = new_comm;
       }
       #pragma omp barrier
 
-      MPI_Comm return_comm = passage[omp_get_thread_num()];
+      ep_lib::MPI_Comm return_comm = ep_lib::passage[omp_get_thread_num()];
       return return_comm;
         
     }
-    return MPI_COMM_NULL;
+    return EP_COMM_NULL;
 
   }
 
-}
+//}
 
 

XIOS/dev/branch_openmp/extern/src_ep_dev/ep_gatherv.cpp

@@ -87 +87 @@
     }
 
-    MPI_Bcast(recvcounts, ep_size, MPI_INT, root, comm);
-    MPI_Bcast(displs, ep_size, MPI_INT, root, comm);
+    MPI_Bcast(recvcounts, ep_size, EP_INT, root, comm);
+    MPI_Bcast(displs, ep_size, EP_INT, root, comm);
 
-    if(mpi_rank == root_mpi_rank) MPI_Gather_local(&sendcount, 1, MPI_INT, local_recvcounts.data(), root_ep_loc, comm);
-    else                          MPI_Gather_local(&sendcount, 1, MPI_INT, local_recvcounts.data(), 0, comm);
+    if(mpi_rank == root_mpi_rank) MPI_Gather_local(&sendcount, 1, EP_INT, local_recvcounts.data(), root_ep_loc, comm);
+    else                          MPI_Gather_local(&sendcount, 1, EP_INT, local_recvcounts.data(), 0, comm);
 
 

XIOS/dev/branch_openmp/extern/src_ep_dev/ep_intercomm.cpp

@@ -33 +33 @@
     }
     
-    MPI_Bcast(mpi_rank_of_leader, 2, MPI_INT, local_leader, local_comm);
+    MPI_Bcast(mpi_rank_of_leader, 2, EP_INT, local_leader, local_comm);
     
     if(mpi_rank_of_leader[0] != mpi_rank_of_leader[1])
@@ -86 +86 @@
       MPI_Comm_rank(peer_comm, &local_leader_rank_in_peer);
       ::MPI_Comm_rank(to_mpi_comm(peer_comm->mpi_comm), &local_leader_rank_in_peer_mpi);
-      ::MPI_Comm_rank(to_mpi_comm(MPI_COMM_WORLD->mpi_comm), &local_leader_rank_in_world);
+      ::MPI_Comm_rank(to_mpi_comm(EP_COMM_WORLD->mpi_comm), &local_leader_rank_in_world);
 
       send_quadruple[0] = ep_size;
@@ -99 +99 @@
       if(remote_leader > local_leader_rank_in_peer)
       {
-        MPI_Isend(send_quadruple, 4, MPI_INT, remote_leader, tag, peer_comm, &request);
+        MPI_Isend(send_quadruple, 4, EP_INT, remote_leader, tag, peer_comm, &request);
         MPI_Wait(&request, &status);
        
-        MPI_Irecv(recv_quadruple, 4, MPI_INT, remote_leader, tag, peer_comm, &request);
+        MPI_Irecv(recv_quadruple, 4, EP_INT, remote_leader, tag, peer_comm, &request);
         MPI_Wait(&request, &status);
       }
       else
       {
-        MPI_Irecv(recv_quadruple, 4, MPI_INT, remote_leader, tag, peer_comm, &request);
+        MPI_Irecv(recv_quadruple, 4, EP_INT, remote_leader, tag, peer_comm, &request);
         MPI_Wait(&request, &status);
           
-        MPI_Isend(send_quadruple, 4, MPI_INT, remote_leader, tag, peer_comm, &request);
+        MPI_Isend(send_quadruple, 4, EP_INT, remote_leader, tag, peer_comm, &request);
         MPI_Wait(&request, &status);
       }
@@ -123 +123 @@
     }
 
-    MPI_Bcast(send_quadruple, 4, MPI_INT, local_leader, local_comm);
-    MPI_Bcast(recv_quadruple, 4, MPI_INT, local_leader, local_comm);
+    MPI_Bcast(send_quadruple, 4, EP_INT, local_leader, local_comm);
+    MPI_Bcast(recv_quadruple, 4, EP_INT, local_leader, local_comm);
 
     if(!is_local_leader)
@@ -152 +152 @@
 
     int rank_in_world;
-    ::MPI_Comm_rank(to_mpi_comm(MPI_COMM_WORLD->mpi_comm), &rank_in_world);
+    ::MPI_Comm_rank(to_mpi_comm(EP_COMM_WORLD->mpi_comm), &rank_in_world);
 
     int *ranks_in_world_local  = new int[ep_size];
     int *ranks_in_world_remote = new int[remote_ep_size];
 
-    MPI_Allgather(&rank_in_world, 1, MPI_INT, ranks_in_world_local, 1, MPI_INT, local_comm);
+    MPI_Allgather(&rank_in_world, 1, EP_INT, ranks_in_world_local, 1, EP_INT, local_comm);
 
     if(is_local_leader)
@@ -166 +166 @@
       if(remote_leader > local_leader_rank_in_peer)
       {
-        MPI_Isend(ranks_in_world_local,  ep_size, MPI_INT, remote_leader, tag, peer_comm, &request);
+        MPI_Isend(ranks_in_world_local,  ep_size, EP_INT, remote_leader, tag, peer_comm, &request);
         MPI_Wait(&request, &status);
        
-        MPI_Irecv(ranks_in_world_remote, remote_ep_size, MPI_INT, remote_leader, tag, peer_comm, &request);
+        MPI_Irecv(ranks_in_world_remote, remote_ep_size, EP_INT, remote_leader, tag, peer_comm, &request);
         MPI_Wait(&request, &status);
       }
       else
       {
-        MPI_Irecv(ranks_in_world_remote, remote_ep_size, MPI_INT, remote_leader, tag, peer_comm, &request);
+        MPI_Irecv(ranks_in_world_remote, remote_ep_size, EP_INT, remote_leader, tag, peer_comm, &request);
         MPI_Wait(&request, &status);
           
-        MPI_Isend(ranks_in_world_local,  ep_size, MPI_INT, remote_leader, tag, peer_comm, &request);
+        MPI_Isend(ranks_in_world_local,  ep_size, EP_INT, remote_leader, tag, peer_comm, &request);
         MPI_Wait(&request, &status);
       }
@@ -185 +185 @@
     }
 
-    MPI_Bcast(ranks_in_world_remote, remote_ep_size, MPI_INT, local_leader, local_comm);
+    MPI_Bcast(ranks_in_world_remote, remote_ep_size, EP_INT, local_leader, local_comm);
 
 #ifdef _showinfo
@@ -313 +313 @@
       int *mpi_rank_list = new int[mpi_size];
 
-      ::MPI_Allgather(&ownership, 1, to_mpi_type(MPI_INT), ownership_list, 1, to_mpi_type(MPI_INT), to_mpi_comm(local_comm->mpi_comm));
-      ::MPI_Allgather(&mpi_rank,  1, to_mpi_type(MPI_INT), mpi_rank_list,  1, to_mpi_type(MPI_INT), to_mpi_comm(local_comm->mpi_comm));
+      ::MPI_Allgather(&ownership, 1, to_mpi_type(EP_INT), ownership_list, 1, to_mpi_type(EP_INT), to_mpi_comm(local_comm->mpi_comm));
+      ::MPI_Allgather(&mpi_rank,  1, to_mpi_type(EP_INT), mpi_rank_list,  1, to_mpi_type(EP_INT), to_mpi_comm(local_comm->mpi_comm));
 
       
@@ -347 +347 @@
       }
 
-      ::MPI_Bcast(&local_leader_rank_in_extracted_comm, 1, to_mpi_type(MPI_INT), local_comm->ep_rank_map->at(local_leader).second, to_mpi_comm(local_comm->mpi_comm));
+      ::MPI_Bcast(&local_leader_rank_in_extracted_comm, 1, to_mpi_type(EP_INT), local_comm->ep_rank_map->at(local_leader).second, to_mpi_comm(local_comm->mpi_comm));
 
       ::MPI_Comm *intracomm = new ::MPI_Comm;
-      bool is_real_involved = ownership && extracted_comm != to_mpi_comm(MPI_COMM_NULL->mpi_comm);
+      bool is_real_involved = ownership && extracted_comm != to_mpi_comm(EP_COMM_NULL->mpi_comm);
 
       if(is_real_involved)
@@ -489 +489 @@
     int *remote_rank_info = new int[2*remote_ep_size];
 
-    MPI_Allgather(rank_info, 2, MPI_INT, local_rank_info, 2, MPI_INT, local_comm);
+    MPI_Allgather(rank_info, 2, EP_INT, local_rank_info, 2, EP_INT, local_comm);
 
     if(is_local_leader)
@@ -498 +498 @@
       if(priority)
       {
-        MPI_Isend(local_rank_info, 2*ep_size, MPI_INT, remote_leader, tag, peer_comm, &request);
+        MPI_Isend(local_rank_info, 2*ep_size, EP_INT, remote_leader, tag, peer_comm, &request);
         MPI_Wait(&request, &status);
        
-        MPI_Irecv(remote_rank_info, 2*remote_ep_size, MPI_INT, remote_leader, tag, peer_comm, &request);
+        MPI_Irecv(remote_rank_info, 2*remote_ep_size, EP_INT, remote_leader, tag, peer_comm, &request);
         MPI_Wait(&request, &status);
       }
       else
       {
-        MPI_Irecv(remote_rank_info, 2*remote_ep_size, MPI_INT, remote_leader, tag, peer_comm, &request);
-        MPI_Wait(&request, &status);
-          
-        MPI_Isend(local_rank_info, 2*ep_size, MPI_INT, remote_leader, tag, peer_comm, &request);
-        MPI_Wait(&request, &status);
-      }
-    }
-
-    MPI_Bcast(remote_rank_info, 2*remote_ep_size, MPI_INT, local_leader, local_comm);
+        MPI_Irecv(remote_rank_info, 2*remote_ep_size, EP_INT, remote_leader, tag, peer_comm, &request);
+        MPI_Wait(&request, &status);
+        
+        MPI_Isend(local_rank_info, 2*ep_size, EP_INT, remote_leader, tag, peer_comm, &request);
+        MPI_Wait(&request, &status);
+      }
+    }
+
+    MPI_Bcast(remote_rank_info, 2*remote_ep_size, EP_INT, local_leader, local_comm);
 
     for(int i=0; i<remote_ep_size; i++)
@@ -537 +537 @@
     }
     */
-    
+   //MPI_Barrier(*newintercomm); 
+   //MPI_Barrier(*newintercomm); 
+   //MPI_Barrier(*newintercomm); 
  
   }

XIOS/dev/branch_openmp/extern/src_ep_dev/ep_lib.cpp

@@ -131 +131 @@
   bool valid_type(MPI_Datatype datatype)
   {
-    if(   datatype == MPI_INT 
-       || datatype == MPI_FLOAT
-       || datatype == MPI_DOUBLE 
-       || datatype == MPI_CHAR 
-       || datatype == MPI_LONG 
-       || datatype == MPI_UNSIGNED_LONG
-       || datatype == MPI_LONG_LONG)
+    if(   datatype == EP_INT 
+       || datatype == EP_FLOAT
+       || datatype == EP_DOUBLE 
+       || datatype == EP_CHAR 
+       || datatype == EP_LONG 
+       || datatype == EP_UNSIGNED_LONG
+       || datatype == EP_LONG_LONG)
     {
       return true;
@@ -148 +148 @@
   bool valid_op(MPI_Op op)
   {
-    if(   op == MPI_MAX 
-       || op == MPI_MIN
-       || op == MPI_SUM
-       || op == MPI_LOR)
+    if(   op == EP_MAX 
+       || op == EP_MIN
+       || op == EP_SUM
+       || op == EP_LOR)
     {
       return true;

XIOS/dev/branch_openmp/extern/src_ep_dev/ep_lib.hpp

@@ -51 +51 @@
   int MPI_Imrecv(void *buf, int count, MPI_Datatype datatype, MPI_Message *message, MPI_Request *request);
 
+  int my_recv  (void *buf, int count, MPI_Datatype datatype, int src, int tag, MPI_Comm comm, MPI_Status *status);
+  int my_irecv (void *buf, int count, MPI_Datatype datatype, int src, int tag, MPI_Comm comm, MPI_Request *request);
 
   int MPI_Test(MPI_Request *request, int *flag, MPI_Status *status);
   int MPI_Testall(int count, MPI_Request *array_of_requests, int *flag, MPI_Status *array_of_statuses);
 
+  int my_test(MPI_Request *request, int *flag, MPI_Status *status);
+  int my_testall(int count, MPI_Request *array_of_requests, int *flag, MPI_Status *array_of_statuses);
+  
   int MPI_Wait(MPI_Request *request, MPI_Status *status);
   int MPI_Waitall(int count, MPI_Request *array_of_requests, MPI_Status *array_of_statuses);
 
+  int my_wait(MPI_Request *request, MPI_Status *status);
+  int my_waitall(int count, MPI_Request *array_of_requests, MPI_Status *array_of_statuses);
 
   int MPI_Iprobe(int source, int tag, MPI_Comm comm, int *flag, MPI_Status *status);

XIOS/dev/branch_openmp/extern/src_ep_dev/ep_lib_fortran.hpp

@@ -4 +4 @@
 #include "ep_type.hpp"
 
-namespace ep_lib
-{
+//namespace ep_lib
+//{
 
-  void* EP_Comm_c2f(MPI_Comm comm);
-  MPI_Comm EP_Comm_f2c(void* comm);
-}
+  void* EP_Comm_c2f(ep_lib::MPI_Comm comm);
+  ep_lib::MPI_Comm EP_Comm_f2c(void* comm);
+//}
 
 

XIOS/dev/branch_openmp/extern/src_ep_dev/ep_merge.cpp

@@ -36 +36 @@
     int sum_high;
     
-    MPI_Allreduce(&int_high, &sum_high, 1, MPI_INT, MPI_SUM, *newintracomm);
+    MPI_Allreduce(&int_high, &sum_high, 1, EP_INT, EP_SUM, *newintracomm);
 
     if(sum_high==0 || sum_high==ep_size+remote_ep_size)
@@ -55 +55 @@
 
     
-    MPI_Allgather(my_triple, 3, MPI_INT, my_triple_list, 3, MPI_INT, *newintracomm); 
+    MPI_Allgather(my_triple, 3, EP_INT, my_triple_list, 3, EP_INT, *newintracomm); 
 
 

XIOS/dev/branch_openmp/extern/src_ep_dev/ep_message.cpp

@@ -147 +147 @@
     {
       Debug("Message probing for intracomm\n");
-     
+/*     
       #pragma omp critical (_mpi_call)
       {
@@ -158 +158 @@
       }
 
-      
+*/    
+      ::MPI_Improbe(MPI_ANY_SOURCE, MPI_ANY_TAG, to_mpi_comm(comm->mpi_comm), &flag, &message, &status);      
+
       if(flag)
       {

XIOS/dev/branch_openmp/extern/src_ep_dev/ep_mpi.hpp

@@ -1 +1 @@
 #ifndef EP_MPI_HPP_INCLUDED
 #define EP_MPI_HPP_INCLUDED
-
+#ifdef _usingEP
 #include "ep_type.hpp"
-
+#endif
 MPI_Datatype to_mpi_type(ep_lib::MPI_Datatype type);
 MPI_Op       to_mpi_op(ep_lib::MPI_Op op);

XIOS/dev/branch_openmp/extern/src_ep_dev/ep_probe.cpp

@@ -68 +68 @@
     *flag = false;
     
-    Message_Check(comm);
-
-    #pragma omp flush
-
     #pragma omp critical (_query)
     for(Message_list::iterator it = comm->ep_comm_ptr->message_queue->begin(); it!= comm->ep_comm_ptr->message_queue->end(); ++it)
@@ -98 +94 @@
       }
     }
+    if(*flag) return 0;
+    
+    Message_Check(comm);
+
+    #pragma omp flush
+
+    #pragma omp critical (_query)
+    for(Message_list::iterator it = comm->ep_comm_ptr->message_queue->begin(); it!= comm->ep_comm_ptr->message_queue->end(); ++it)
+    {
+      bool src_matched = src<0? true: (*it)->ep_src == src;
+      bool tag_matched = tag<0? true: (*it)->ep_tag == tag;
+      
+      if(src_matched && tag_matched)        
+      {
+        Debug("find message\n");
+         
+        status->mpi_status = new ::MPI_Status(*static_cast< ::MPI_Status*>((*it)->mpi_status));
+        status->ep_src = (*it)->ep_src;
+        status->ep_tag = (*it)->ep_tag;
+        
+        if(comm->is_intercomm)
+        {
+          for(INTER_RANK_MAP::iterator iter = comm->inter_rank_map->begin(); iter != comm->inter_rank_map->end(); iter++)
+          {
+            if(iter->second == (*it)->ep_src) status->ep_src=iter->first;
+          }
+        }
+
+        *flag = true;
+        break;
+      }
+    }
+    if(*flag) return 0;
   }
 
@@ -129 +158 @@
     *flag = false;
     
-    Message_Check(comm);
-    
-    #pragma omp flush
-
     #pragma omp critical (_query)
     if(! comm->ep_comm_ptr->message_queue->empty())
@@ -179 +204 @@
       }
     }
+
+    if(*flag) return 0;
+   
+    Message_Check(comm);
+    
+    #pragma omp flush
+
+    #pragma omp critical (_query)
+    if(! comm->ep_comm_ptr->message_queue->empty())
+    {
+      for(Message_list::iterator it = comm->ep_comm_ptr->message_queue->begin(); it!= comm->ep_comm_ptr->message_queue->end(); ++it)
+      {
+                                          
+        bool src_matched = src<0? true: (*it)->ep_src == src;
+        bool tag_matched = tag<0? true: (*it)->ep_tag == tag;
+        
+        if(src_matched && tag_matched)
+        {
+          *flag = true;
+
+          status->mpi_status = new ::MPI_Status(*static_cast< ::MPI_Status*>((*it)->mpi_status));
+          memcheck("new "<< status->mpi_status << " : in ep_lib::MPI_Improbe, status->mpi_status = new ::MPI_Status");
+          status->ep_src = (*it)->ep_src;
+          status->ep_tag = (*it)->ep_tag;
+
+          (*message)->mpi_message = new ::MPI_Message(*static_cast< ::MPI_Message*>((*it)->mpi_message));
+          memcheck("new "<< (*message)->mpi_message <<" : in ep_lib::MPI_Improbe, (*message)->mpi_message = new ::MPI_Message");
+          (*message)->ep_src = (*it)->ep_src;
+          (*message)->ep_tag = (*it)->ep_tag;
+                                      
+
+          #pragma omp critical (_query2)
+          {              
+            memcheck("delete "<< (*it)->mpi_message <<" : in ep_lib::Message_Check, delete (*it)->mpi_message");
+            memcheck("delete "<< (*it)->mpi_status <<" : in ep_lib::Message_Check, delete (*it)->mpi_status");
+            memcheck("delete "<< (*it) <<" : in ep_lib::Message_Check, delete (*it)");
+            
+            
+            delete (*it)->mpi_message;     
+            delete (*it)->mpi_status; 
+            delete *it;
+            
+                       
+            comm->ep_comm_ptr->message_queue->erase(it);
+            memcheck("message_queue["<<mpi_rank<<","<<ep_rank_loc<<"]->size = "<<comm->ep_comm_ptr->message_queue->size());
+            #pragma omp flush
+          }
+          
+          break;
+        }
+
+      }
+    }
+
+    if(*flag) return 0;
+
   }
 

XIOS/dev/branch_openmp/extern/src_ep_dev/ep_recv.cpp

@@ -22 +22 @@
   int MPI_Recv(void *buf, int count, MPI_Datatype datatype, int src, int tag, MPI_Comm comm, MPI_Status *status)
   {
-
     if(!comm->is_ep) return MPI_Recv_mpi(buf, count, datatype, src, tag, comm, status);
     

XIOS/dev/branch_openmp/extern/src_ep_dev/ep_reduce.cpp

@@ -80 +80 @@
       memcpy(recvbuf, comm->my_buffer->void_buffer[0], datasize * count);
 
-      if(op == MPI_MAX)
-      {
-        if(datatype == MPI_INT)
+      if(op == EP_MAX)
+      {
+        if(datatype == EP_INT)
         {
           assert(datasize == sizeof(int));
@@ -96 +96 @@
         }
 
-        else if(datatype == MPI_DOUBLE)
+        else if(datatype == EP_DOUBLE)
         {
           assert(datasize == sizeof(double));
@@ -103 +103 @@
         }
 
-        else if(datatype == MPI_CHAR)
+        else if(datatype == EP_CHAR)
         {
           assert(datasize == sizeof(char));
@@ -110 +110 @@
         }
 
-        else if(datatype == MPI_LONG)
+        else if(datatype == EP_LONG)
         {
           assert(datasize == sizeof(long));
@@ -117 +117 @@
         }
 
-        else if(datatype == MPI_UNSIGNED_LONG)
+        else if(datatype == EP_UNSIGNED_LONG)
         {
           assert(datasize == sizeof(unsigned long));
@@ -124 +124 @@
         }
 
-        else if(datatype == MPI_LONG_LONG_INT)
+        else if(datatype == EP_LONG_LONG_INT)
         {
           assert(datasize == sizeof(long long));
@@ -139 +139 @@
       }
 
-      else if(op == MPI_MIN)
-      {
-        if(datatype ==MPI_INT)
+      else if(op == EP_MIN)
+      {
+        if(datatype ==EP_INT)
         {
           assert(datasize == sizeof(int));
@@ -155 +155 @@
         }
 
-        else if(datatype == MPI_DOUBLE)
+        else if(datatype == EP_DOUBLE)
         {
           assert(datasize == sizeof(double));
@@ -162 +162 @@
         }
 
-        else if(datatype == MPI_CHAR)
+        else if(datatype == EP_CHAR)
         {
           assert(datasize == sizeof(char));
@@ -169 +169 @@
         }
 
-        else if(datatype == MPI_LONG)
+        else if(datatype == EP_LONG)
         {
           assert(datasize == sizeof(long));
@@ -176 +176 @@
         }
 
-        else if(datatype == MPI_UNSIGNED_LONG)
+        else if(datatype == EP_UNSIGNED_LONG)
         {
           assert(datasize == sizeof(unsigned long));
@@ -183 +183 @@
         }
 
-        else if(datatype == MPI_LONG_LONG_INT)
+        else if(datatype == EP_LONG_LONG_INT)
         {
           assert(datasize == sizeof(long long));
@@ -199 +199 @@
 
 
-      else if(op == MPI_SUM)
-      {
-        if(datatype==MPI_INT)
+      else if(op == EP_SUM)
+      {
+        if(datatype==EP_INT)
         {
           assert(datasize == sizeof(int));
@@ -215 +215 @@
         }
 
-        else if(datatype == MPI_DOUBLE)
+        else if(datatype == EP_DOUBLE)
         {
           assert(datasize == sizeof(double));
@@ -222 +222 @@
         }
 
-        else if(datatype == MPI_CHAR)
+        else if(datatype == EP_CHAR)
         {
           assert(datasize == sizeof(char));
@@ -229 +229 @@
         }
 
-        else if(datatype == MPI_LONG)
+        else if(datatype == EP_LONG)
         {
           assert(datasize == sizeof(long));
@@ -236 +236 @@
         }
 
-        else if(datatype ==MPI_UNSIGNED_LONG)
+        else if(datatype ==EP_UNSIGNED_LONG)
         {
           assert(datasize == sizeof(unsigned long));
@@ -243 +243 @@
         }
         
-        else if(datatype ==MPI_LONG_LONG_INT)
+        else if(datatype ==EP_LONG_LONG_INT)
         {
           assert(datasize == sizeof(long long));
@@ -258 +258 @@
       }
 
-      else if(op == MPI_LOR)
-      {
-        if(datatype != MPI_INT)
+      else if(op == EP_LOR)
+      {
+        if(datatype != EP_INT)
           printf("datatype Error, must be MPI_INT\n");
         else

XIOS/dev/branch_openmp/extern/src_ep_dev/ep_reduce_scatter.cpp

@@ -61 +61 @@
 
     int my_recvcount = recvcounts[ep_rank];
-    MPI_Gather_local(&my_recvcount, 1, MPI_INT, local_recvcounts.data(), 0, comm);
-    MPI_Bcast_local(local_recvcounts.data(), num_ep, MPI_INT, 0, comm);
+    MPI_Gather_local(&my_recvcount, 1, EP_INT, local_recvcounts.data(), 0, comm);
+    MPI_Bcast_local(local_recvcounts.data(), num_ep, EP_INT, 0, comm);
 
     int my_displs = std::accumulate(recvcounts, recvcounts+ep_rank, 0);
-    MPI_Gather_local(&my_displs, 1, MPI_INT, local_displs.data(), 0, comm);
-    MPI_Bcast_local(local_displs.data(), num_ep, MPI_INT, 0, comm);
+    MPI_Gather_local(&my_displs, 1, EP_INT, local_displs.data(), 0, comm);
+    MPI_Bcast_local(local_displs.data(), num_ep, EP_INT, 0, comm);
 
     

XIOS/dev/branch_openmp/extern/src_ep_dev/ep_scan.cpp

@@ -60 +60 @@
     if(ep_rank_loc == 0 && mpi_rank != 0)
     {
-      if(op == MPI_SUM)
-      {
-        if(datatype == MPI_INT)
+      if(op == EP_SUM)
+      {
+        if(datatype == EP_INT)
         {
           assert(datasize == sizeof(int));
@@ -68 +68 @@
         }
           
-        else if(datatype == MPI_FLOAT)
+        else if(datatype == EP_FLOAT)
         {
           assert( datasize == sizeof(float));
@@ -74 +74 @@
         } 
              
-        else if(datatype == MPI_DOUBLE )
+        else if(datatype == EP_DOUBLE )
         {
           assert( datasize == sizeof(double));
@@ -80 +80 @@
         }
       
-        else if(datatype == MPI_CHAR)
+        else if(datatype == EP_CHAR)
         {
           assert( datasize == sizeof(char));
@@ -86 +86 @@
         } 
           
-        else if(datatype == MPI_LONG)
+        else if(datatype == EP_LONG)
         {
           assert( datasize == sizeof(long));
@@ -93 +93 @@
           
             
-        else if(datatype == MPI_UNSIGNED_LONG)
+        else if(datatype == EP_UNSIGNED_LONG)
         {
           assert(datasize == sizeof(unsigned long));
@@ -99 +99 @@
         }
         
-        else if(datatype == MPI_LONG_LONG_INT)
+        else if(datatype == EP_LONG_LONG_INT)
         {
           assert(datasize == sizeof(long long int));
@@ -112 +112 @@
       }
 
-      else if(op == MPI_MAX)
-      {
-        if(datatype == MPI_INT)
+      else if(op == EP_MAX)
+      {
+        if(datatype == EP_INT)
         {
           assert( datasize == sizeof(int));
@@ -120 +120 @@
         } 
           
-        else if(datatype == MPI_FLOAT )
+        else if(datatype == EP_FLOAT )
         {
           assert( datasize == sizeof(float));
@@ -126 +126 @@
         }
 
-        else if(datatype == MPI_DOUBLE )
+        else if(datatype == EP_DOUBLE )
         {
           assert( datasize == sizeof(double));
@@ -132 +132 @@
         }
       
-        else if(datatype == MPI_CHAR )
+        else if(datatype == EP_CHAR )
         {
           assert(datasize == sizeof(char));
@@ -138 +138 @@
         }
       
-        else if(datatype == MPI_LONG)
+        else if(datatype == EP_LONG)
         {
           assert( datasize == sizeof(long));
@@ -144 +144 @@
         } 
             
-        else if(datatype == MPI_UNSIGNED_LONG)
+        else if(datatype == EP_UNSIGNED_LONG)
         {
           assert( datasize == sizeof(unsigned long));
@@ -150 +150 @@
         } 
             
-        else if(datatype == MPI_LONG_LONG_INT)
+        else if(datatype == EP_LONG_LONG_INT)
         {
           assert(datasize == sizeof(long long int));
@@ -163 +163 @@
       }
 
-      else if(op == MPI_MIN)
-      {
-        if(datatype == MPI_INT )
+      else if(op == EP_MIN)
+      {
+        if(datatype == EP_INT )
         {
           assert (datasize == sizeof(int));
@@ -171 +171 @@
         }
           
-        else if(datatype == MPI_FLOAT )
+        else if(datatype == EP_FLOAT )
         {
           assert( datasize == sizeof(float));
@@ -177 +177 @@
         }
              
-        else if(datatype == MPI_DOUBLE )
+        else if(datatype == EP_DOUBLE )
         {
           assert( datasize == sizeof(double));
@@ -183 +183 @@
         }
       
-        else if(datatype == MPI_CHAR )
+        else if(datatype == EP_CHAR )
         {
           assert( datasize == sizeof(char));
@@ -189 +189 @@
         }
       
-        else if(datatype == MPI_LONG )
+        else if(datatype == EP_LONG )
         { 
           assert( datasize == sizeof(long));
@@ -195 +195 @@
         }
             
-        else if(datatype == MPI_UNSIGNED_LONG )
+        else if(datatype == EP_UNSIGNED_LONG )
         {
           assert( datasize == sizeof(unsigned long));
@@ -201 +201 @@
         }
             
-        else if(datatype == MPI_LONG_LONG_INT)
+        else if(datatype == EP_LONG_LONG_INT)
         {
           assert(datasize == sizeof(long long int));
@@ -235 +235 @@
 
 
-    if(op == MPI_SUM)
-    {
-      if(datatype == MPI_INT )
+    if(op == EP_SUM)
+    {
+      if(datatype == EP_INT )
       {
         assert (datasize == sizeof(int));
@@ -244 +244 @@
       }
      
-      else if(datatype == MPI_FLOAT )
+      else if(datatype == EP_FLOAT )
       {
         assert(datasize == sizeof(float));
@@ -252 +252 @@
       
 
-      else if(datatype == MPI_DOUBLE )
+      else if(datatype == EP_DOUBLE )
       {
         assert(datasize == sizeof(double));
@@ -259 +259 @@
       }
 
-      else if(datatype == MPI_CHAR )
+      else if(datatype == EP_CHAR )
       {
         assert(datasize == sizeof(char));
@@ -266 +266 @@
       }
 
-      else if(datatype == MPI_LONG )
+      else if(datatype == EP_LONG )
       {
         assert(datasize == sizeof(long));
@@ -273 +273 @@
       }
 
-      else if(datatype == MPI_UNSIGNED_LONG )
+      else if(datatype == EP_UNSIGNED_LONG )
       {
         assert(datasize == sizeof(unsigned long));
@@ -280 +280 @@
       }
       
-      else if(datatype == MPI_LONG_LONG_INT )
+      else if(datatype == EP_LONG_LONG_INT )
       {
         assert(datasize == sizeof(long long int));
@@ -296 +296 @@
     }
 
-    else if(op == MPI_MAX)
-    {
-      if(datatype == MPI_INT)
+    else if(op == EP_MAX)
+    {
+      if(datatype == EP_INT)
       {
         assert(datasize == sizeof(int));
@@ -305 +305 @@
       }
 
-      else if(datatype == MPI_FLOAT )
+      else if(datatype == EP_FLOAT )
       {
         assert(datasize == sizeof(float));
@@ -312 +312 @@
       }
 
-      else if(datatype == MPI_DOUBLE )
+      else if(datatype == EP_DOUBLE )
       {
         assert(datasize == sizeof(double));
@@ -319 +319 @@
       }
 
-      else if(datatype == MPI_CHAR )
+      else if(datatype == EP_CHAR )
       {
         assert(datasize == sizeof(char));
@@ -326 +326 @@
       }
 
-      else if(datatype == MPI_LONG )
+      else if(datatype == EP_LONG )
       {
         assert(datasize == sizeof(long));
@@ -333 +333 @@
       }
 
-      else if(datatype == MPI_UNSIGNED_LONG )
+      else if(datatype == EP_UNSIGNED_LONG )
       {
         assert(datasize == sizeof(unsigned long));
@@ -340 +340 @@
       }
      
-      else if(datatype == MPI_LONG_LONG_INT )
+      else if(datatype == EP_LONG_LONG_INT )
       {
         assert(datasize == sizeof(long long int));
@@ -355 +355 @@
     }
 
-    else if(op == MPI_MIN)
-    {
-      if(datatype == MPI_INT )
+    else if(op == EP_MIN)
+    {
+      if(datatype == EP_INT )
       {
         assert(datasize == sizeof(int));
@@ -364 +364 @@
       }
 
-      else if(datatype == MPI_FLOAT )
+      else if(datatype == EP_FLOAT )
       {
         assert(datasize == sizeof(float));
@@ -371 +371 @@
       }
 
-      else if(datatype == MPI_DOUBLE )
+      else if(datatype == EP_DOUBLE )
       {
         assert(datasize == sizeof(double));
@@ -378 +378 @@
       }
 
-      else if(datatype == MPI_CHAR )
+      else if(datatype == EP_CHAR )
       {
         assert(datasize == sizeof(char));
@@ -385 +385 @@
       }
 
-      else if(datatype == MPI_LONG )
+      else if(datatype == EP_LONG )
       {
         assert(datasize == sizeof(long));
@@ -392 +392 @@
       }
 
-      else if(datatype == MPI_UNSIGNED_LONG )
+      else if(datatype == EP_UNSIGNED_LONG )
       {
         assert(datasize == sizeof(unsigned long));
@@ -399 +399 @@
       }
 
-      else if(datatype == MPI_LONG_LONG_INT )
+      else if(datatype == EP_LONG_LONG_INT )
       {
         assert(datasize == sizeof(long long int));

XIOS/dev/branch_openmp/extern/src_ep_dev/ep_scatter.cpp

@@ -73 +73 @@
     std::vector<int>ranks(ep_size);
 
-    if(mpi_rank == root_mpi_rank) MPI_Gather_local(&ep_rank, 1, MPI_INT, local_ranks.data(), root_ep_loc, comm);
-    else                          MPI_Gather_local(&ep_rank, 1, MPI_INT, local_ranks.data(), 0, comm);
+    if(mpi_rank == root_mpi_rank) MPI_Gather_local(&ep_rank, 1, EP_INT, local_ranks.data(), root_ep_loc, comm);
+    else                          MPI_Gather_local(&ep_rank, 1, EP_INT, local_ranks.data(), 0, comm);
 
 
@@ -91 +91 @@
         displs[i] = displs[i-1] + recvcounts[i-1];
 
-      ::MPI_Gatherv(local_ranks.data(), num_ep, to_mpi_type(MPI_INT), ranks.data(), recvcounts.data(), displs.data(), to_mpi_type(MPI_INT), root_mpi_rank, to_mpi_comm(comm->mpi_comm));
+      ::MPI_Gatherv(local_ranks.data(), num_ep, to_mpi_type(EP_INT), ranks.data(), recvcounts.data(), displs.data(), to_mpi_type(EP_INT), root_mpi_rank, to_mpi_comm(comm->mpi_comm));
     }
 
@@ -109 +109 @@
     {
       int local_sendcount = num_ep * count;
-      ::MPI_Gather(&local_sendcount, 1, to_mpi_type(MPI_INT), recvcounts.data(), 1, to_mpi_type(MPI_INT), root_mpi_rank, to_mpi_comm(comm->mpi_comm));
+      ::MPI_Gather(&local_sendcount, 1, to_mpi_type(EP_INT), recvcounts.data(), 1, to_mpi_type(EP_INT), root_mpi_rank, to_mpi_comm(comm->mpi_comm));
       
       if(is_root) for(int i=1; i<mpi_size; i++) displs[i] = displs[i-1] + recvcounts[i-1];

XIOS/dev/branch_openmp/extern/src_ep_dev/ep_scatterv.cpp

@@ -76 +76 @@
     std::vector<int>ranks(ep_size);
 
-    if(mpi_rank == root_mpi_rank) MPI_Gather_local(&ep_rank, 1, MPI_INT, local_ranks.data(), root_ep_loc, comm);
-    else                          MPI_Gather_local(&ep_rank, 1, MPI_INT, local_ranks.data(), 0, comm);
+    if(mpi_rank == root_mpi_rank) MPI_Gather_local(&ep_rank, 1, EP_INT, local_ranks.data(), root_ep_loc, comm);
+    else                          MPI_Gather_local(&ep_rank, 1, EP_INT, local_ranks.data(), 0, comm);
 
 
@@ -94 +94 @@
         my_displs[i] = my_displs[i-1] + recvcounts[i-1];
 
-      ::MPI_Gatherv(local_ranks.data(), num_ep, to_mpi_type(MPI_INT), ranks.data(), recvcounts.data(), my_displs.data(), to_mpi_type(MPI_INT), root_mpi_rank, to_mpi_comm(comm->mpi_comm));
+      ::MPI_Gatherv(local_ranks.data(), num_ep, to_mpi_type(EP_INT), ranks.data(), recvcounts.data(), my_displs.data(), to_mpi_type(EP_INT), root_mpi_rank, to_mpi_comm(comm->mpi_comm));
     }
 
@@ -112 +112 @@
     void* local_sendbuf;
     int local_sendcount;
-    if(mpi_rank == root_mpi_rank) MPI_Reduce_local(&recvcount, &local_sendcount, 1, MPI_INT, MPI_SUM, root_ep_loc, comm);
-    else                          MPI_Reduce_local(&recvcount, &local_sendcount, 1, MPI_INT, MPI_SUM, 0, comm);
+    if(mpi_rank == root_mpi_rank) MPI_Reduce_local(&recvcount, &local_sendcount, 1, EP_INT, EP_SUM, root_ep_loc, comm);
+    else                          MPI_Reduce_local(&recvcount, &local_sendcount, 1, EP_INT, EP_SUM, 0, comm);
 
     if(is_master) 
@@ -119 +119 @@
       local_sendbuf = new void*[datasize * local_sendcount];
  
-      ::MPI_Gather(&local_sendcount, 1, to_mpi_type(MPI_INT), recvcounts.data(), 1, to_mpi_type(MPI_INT), root_mpi_rank, to_mpi_comm(comm->mpi_comm));
+      ::MPI_Gather(&local_sendcount, 1, to_mpi_type(EP_INT), recvcounts.data(), 1, to_mpi_type(EP_INT), root_mpi_rank, to_mpi_comm(comm->mpi_comm));
 
       if(is_root) for(int i=1; i<mpi_size; i++) my_displs[i] = my_displs[i-1] + recvcounts[i-1];
@@ -129 +129 @@
     std::vector<int>local_displs(num_ep, 0);
 
-    MPI_Gather_local(&recvcount, 1, MPI_INT, local_sendcounts.data(), 0, comm);
-    MPI_Bcast_local(local_sendcounts.data(), num_ep, MPI_INT, 0, comm);
+    MPI_Gather_local(&recvcount, 1, EP_INT, local_sendcounts.data(), 0, comm);
+    MPI_Bcast_local(local_sendcounts.data(), num_ep, EP_INT, 0, comm);
     for(int i=1; i<num_ep; i++)
       local_displs[i] = local_displs[i-1] + local_sendcounts[i-1];

XIOS/dev/branch_openmp/extern/src_ep_dev/ep_split.cpp

@@ -54 +54 @@
     vector<int> all_color_loc(num_ep);
 
-    MPI_Allgather(&color, 1, MPI_INT, all_color.data(), 1, MPI_INT, comm);
-    MPI_Allgather_local(&color, 1, MPI_INT, all_color_loc.data(), comm);
+    MPI_Allgather(&color, 1, EP_INT, all_color.data(), 1, EP_INT, comm);
+    MPI_Allgather_local(&color, 1, EP_INT, all_color_loc.data(), comm);
 
     list<int> color_list(all_color.begin(), all_color.end());
@@ -106 +106 @@
     
 
-    MPI_Allgather(&key, 1, MPI_INT, all_key.data(),1, MPI_INT, comm);
-    MPI_Allgather_local(&key, 1, MPI_INT, all_key_loc.data(), comm);
+    MPI_Allgather(&key, 1, EP_INT, all_key.data(),1, EP_INT, comm);
+    MPI_Allgather_local(&key, 1, EP_INT, all_key_loc.data(), comm);
     
     vector<int> all_ep_rank(ep_size);
@@ -115 +115 @@
     vector<int> colored_ep_rank_loc[num_color];
     
-    MPI_Allgather(&ep_rank, 1, MPI_INT, all_ep_rank.data(),1, MPI_INT, comm);
-    MPI_Allgather_local(&ep_rank, 1, MPI_INT, all_ep_rank_loc.data(), comm);
+    MPI_Allgather(&ep_rank, 1, EP_INT, all_ep_rank.data(),1, EP_INT, comm);
+    MPI_Allgather_local(&ep_rank, 1, EP_INT, all_ep_rank_loc.data(), comm);
 
     for(int i=0; i<num_ep; i++)
@@ -270 +270 @@
         if(new_ep_rank_loc == 0) my_triple_vector_recv.resize(3*new_ep_size);
         
-        MPI_Gather_local(my_triple, 3, MPI_INT, my_triple_vector.data(), 0, *newcomm);
+        MPI_Gather_local(my_triple, 3, EP_INT, my_triple_vector.data(), 0, *newcomm);
         
         if(new_ep_rank_loc == 0)
@@ -277 +277 @@
           int *displs = new int[new_mpi_size];
           int new_num_epx3 = new_num_ep * 3;
-          ::MPI_Allgather(&new_num_epx3, 1, to_mpi_type(MPI_INT), recvcounts, 1, to_mpi_type(MPI_INT), to_mpi_comm((*newcomm)->mpi_comm));
+          ::MPI_Allgather(&new_num_epx3, 1, to_mpi_type(EP_INT), recvcounts, 1, to_mpi_type(EP_INT), to_mpi_comm((*newcomm)->mpi_comm));
           displs[0]=0;
           for(int i=1; i<new_mpi_size; i++)
             displs[i] = displs[i-1] + recvcounts[i-1];
              
-          ::MPI_Allgatherv(my_triple_vector.data(), 3*new_num_ep, to_mpi_type(MPI_INT), my_triple_vector_recv.data(), recvcounts, displs, to_mpi_type(MPI_INT), to_mpi_comm((*newcomm)->mpi_comm));
+          ::MPI_Allgatherv(my_triple_vector.data(), 3*new_num_ep, to_mpi_type(EP_INT), my_triple_vector_recv.data(), recvcounts, displs, to_mpi_type(EP_INT), to_mpi_comm((*newcomm)->mpi_comm));
           
           for(int i=0; i<new_ep_size; i++)
           {
-            (*newcomm)->ep_comm_ptr->comm_list[0]->ep_rank_map->insert(std::pair< int, std::pair<int,int> >(my_triple_vector_recv[3*i], my_triple_vector_recv[3*i+1], my_triple_vector_recv[3*i+2]));
+            //(*newcomm)->ep_comm_ptr->comm_list[0]->ep_rank_map->insert(std::pair< int, std::pair<int,int> >(my_triple_vector_recv[3*i], my_triple_vector_recv[3*i+1], my_triple_vector_recv[3*i+2]));
+            (*((*newcomm)->ep_comm_ptr->comm_list[0]->ep_rank_map)) [ my_triple_vector_recv[3*i] ] = std::make_pair(my_triple_vector_recv[3*i+1], my_triple_vector_recv[3*i+2]);
           }
           

XIOS/dev/branch_openmp/extern/src_ep_dev/ep_type.hpp

@@ -68 +68 @@
     
     MPI_Fint() {}
-    MPI_Fint(void* fint);
+    MPI_Fint(void* fint) {mpi_fint = fint;};
                                  
   };

XIOS/dev/branch_openmp/extern/src_ep_dev/ep_win.cpp

@@ -21 +21 @@
 
     int num_ep_max;
-    MPI_Allreduce(&num_ep, &num_ep_max, 1, MPI_INT, MPI_MAX, comm);
+    MPI_Allreduce(&num_ep, &num_ep_max, 1, EP_INT, EP_MAX, comm);
 
     assert(num_ep_max > 1);
@@ -88 +88 @@
     
     int num_ep_max;
-    MPI_Allreduce(&num_ep, &num_ep_max, 1, MPI_INT, MPI_MAX, (*win)->comm);
+    MPI_Allreduce(&num_ep, &num_ep_max, 1, EP_INT, EP_MAX, (*win)->comm);
 
     //printf("rank_loc = %d, thread_num = %d, num_ep_max = %d\n", rank_loc, omp_get_thread_num(), num_ep_max);
@@ -139 +139 @@
     
     int num_ep_max;
-    MPI_Allreduce(&num_ep, &num_ep_max, 1, MPI_INT, MPI_MAX, win->comm);
+    MPI_Allreduce(&num_ep, &num_ep_max, 1, EP_INT, EP_MAX, win->comm);
 
     if(num_ep == 1)