Ignore:
Timestamp:
11/15/17 12:14:34 (4 years ago)
Author:
yushan
Message:

dev_omp

Location:
XIOS/dev/branch_openmp/extern/remap/src
Files:
23 edited

Legend:

Unmodified
Added
Removed
  • XIOS/dev/branch_openmp/extern/remap/src/clipper.cpp

    r1205 r1328  
    42744274{ 
    42754275  //The equation of a line in general form (Ax + By + C = 0) 
    4276   //given 2 points (x,y) & (x,y) is ... 
    4277   //(y - y)x + (x - x)y + (y - y)x - (x - x)y = 0 
    4278   //A = (y - y); B = (x - x); C = (y - y)x - (x - x)y 
    4279   //perpendicular distance of point (x,y) = (Ax + By + C)/Sqrt(A + B) 
     4276  //given 2 points (x¹,y¹) & (x²,y²) is ... 
     4277  //(y¹ - y²)x + (x² - x¹)y + (y² - y¹)x¹ - (x² - x¹)y¹ = 0 
     4278  //A = (y¹ - y²); B = (x² - x¹); C = (y² - y¹)x¹ - (x² - x¹)y¹ 
     4279  //perpendicular distance of point (x³,y³) = (Ax³ + By³ + C)/Sqrt(A² + B²) 
    42804280  //see http://en.wikipedia.org/wiki/Perpendicular_distance 
    42814281  double A = double(ln1.Y - ln2.Y); 
  • XIOS/dev/branch_openmp/extern/remap/src/cputime.cpp

    r694 r1328  
    11#include "mpi.hpp" 
     2using namespace ep_lib; 
    23 
    34namespace sphereRemap { 
  • XIOS/dev/branch_openmp/extern/remap/src/elt.hpp

    r1172 r1328  
    5353struct Elt : Polyg 
    5454{ 
    55   Elt() {} 
    56   Elt(const double *bounds_lon, const double *bounds_lat, int max_num_vert) 
    57   { 
    58     int k = 0; 
    59     vertex[k++] = xyz(bounds_lon[0], bounds_lat[0]); 
    60     for (int i = 1; i < max_num_vert; i++) 
    61     { 
    62       vertex[k] = xyz(bounds_lon[i], bounds_lat[i]); 
    63       /* netCDF convention: if first vertex repeats element is finished (at least three vertices == triagle) */ 
    64       if (k >= 3 && squaredist(vertex[k], vertex[0]) < EPS*EPS)  
    65         break; 
    66       /* eliminate zero edges: move to next vertex only if it is different */ 
    67       if (squaredist(vertex[k], vertex[k-1]) > EPS*EPS) 
    68         k++; 
    69       //else cout << "Removed edge " << k << " due to zero length (coinciding endpoints)." << endl ; 
    70     } 
    71     n = k; 
    72     x = barycentre(vertex, n); 
    73   } 
     55        Elt() {} 
     56        Elt(const double *bounds_lon, const double *bounds_lat, int max_num_vert) 
     57        { 
     58                int k = 0; 
     59                vertex[k++] = xyz(bounds_lon[0], bounds_lat[0]); 
     60                for (int i = 1; i < max_num_vert; i++) 
     61                { 
     62                        vertex[k] = xyz(bounds_lon[i], bounds_lat[i]); 
     63                        /* netCDF convention: if first vertex repeats element is finished (at least three vertices == triagle) */ 
     64                        if (k >= 3 && squaredist(vertex[k], vertex[0]) < EPS*EPS)  
     65                                break; 
     66                        /* eliminate zero edges: move to next vertex only if it is different */ 
     67                        if (squaredist(vertex[k], vertex[k-1]) > EPS*EPS) 
     68                                k++; 
     69                        else 
     70                                /* cout << "Removed edge " << k << " due to zero length (coinciding endpoints)." << endl */ ; 
     71                } 
     72                n = k; 
     73                x = barycentre(vertex, n); 
     74        } 
    7475 
    7576        Elt& operator=(const Elt& rhs) 
     
    9596        } 
    9697 
    97   void delete_intersections() 
    98   { 
    99     for (list<Polyg*>::iterator it = this->is.begin(); it != this->is.end(); it++) 
    100     { 
    101       Polyg* poly = *it; 
    102       delete poly; 
    103     } 
    104   } 
     98        void delete_intersections() 
     99        { 
     100                for (list<Polyg*>::iterator it = this->is.begin(); it != this->is.end(); it++) 
     101                { 
     102                        Polyg* poly = *it; 
     103                        delete poly; 
     104                } 
     105        } 
    105106 
    106107  void insert_vertex(int i, const Coord& v) 
  • XIOS/dev/branch_openmp/extern/remap/src/gridRemap.cpp

    r1155 r1328  
    1111 
    1212CRemapGrid srcGrid; 
    13 #pragma omp threadprivate(srcGrid) 
    14  
    1513CRemapGrid tgtGrid; 
    16 #pragma omp threadprivate(tgtGrid) 
    1714 
    1815} 
  • XIOS/dev/branch_openmp/extern/remap/src/gridRemap.hpp

    r1220 r1328  
    1414Coord readPole(std::istream&); 
    1515 
    16  
     16extern CRemapGrid srcGrid; 
     17extern CRemapGrid tgtGrid; 
    1718 
    1819} 
  • XIOS/dev/branch_openmp/extern/remap/src/intersect.cpp

    r1220 r1328  
    1515namespace sphereRemap { 
    1616 
    17 extern CRemapGrid srcGrid; 
    18 #pragma omp threadprivate(srcGrid) 
    19  
    20 extern CRemapGrid tgtGrid; 
    21 #pragma omp threadprivate(tgtGrid) 
    22  
    2317using namespace std; 
    2418 
     
    2721int neighbour_idx(const Elt& a, const Elt& b) 
    2822{ 
    29   for (int i = 0; i < a.n; i++) 
    30   { 
    31     for (int j = 0; j < b.n; j++) 
    32     { 
    33       assert(squaredist(a.vertex[ i       ], b.vertex[ j       ]) > EPS*EPS || 
    34              squaredist(a.vertex[(i+1)%a.n], b.vertex[(j+1)%b.n]) > EPS*EPS); 
    35       if (   squaredist(a.vertex[ i       ], b.vertex[ j           ]) < 1e-13*1e-13 && 
    36              squaredist(a.vertex[(i+1)%a.n], b.vertex[(j+b.n-1)%b.n]) < 1e-13*1e-13) 
    37       { 
    38         return i; 
    39       } 
    40     } 
    41   } 
    42   return NOT_FOUND; 
     23        for (int i = 0; i < a.n; i++) 
     24        { 
     25                for (int j = 0; j < b.n; j++) 
     26                { 
     27                        assert(squaredist(a.vertex[ i       ], b.vertex[ j       ]) > EPS*EPS || 
     28                               squaredist(a.vertex[(i+1)%a.n], b.vertex[(j+1)%b.n]) > EPS*EPS); 
     29                        if (   squaredist(a.vertex[ i       ], b.vertex[ j           ]) < 1e-13*1e-13 && 
     30                               squaredist(a.vertex[(i+1)%a.n], b.vertex[(j+b.n-1)%b.n]) < 1e-13*1e-13) 
     31                        { 
     32                                return i; 
     33                        } 
     34                } 
     35        } 
     36        return NOT_FOUND; 
    4337} 
    4438 
     
    211205bool isNeighbour(Elt& a, const Elt& b) 
    212206{ 
    213   // return neighbour_idx(a, b) != NOT_FOUND; 
     207        // return neighbour_idx(a, b) != NOT_FOUND; 
    214208  return insertNeighbour(a,b,false) ; 
    215209} 
     
    219213void intersect(Elt *a, Elt *b) 
    220214{ 
    221   int na = a->n; /* vertices of a */ 
    222   int nb = b->n; /* vertices of b */ 
    223   Coord *c   = new Coord[na+nb]; 
    224   Coord *c2  = new Coord[na+nb]; 
    225   Coord *xc  = new Coord[na+nb]; 
    226   Coord *xc2 = new Coord[na+nb]; 
    227   Coord gc, gc2; 
    228   double *d = new double[na+nb]; 
    229   double *d2 = new double[na+nb]; 
    230   double are, are2; 
    231   Ipt ipt[NMAX*NMAX]; 
    232   Ipt ipt2[NMAX*NMAX]; 
    233   ptsec(a, b, ipt); 
    234   /* make ipt2 transpose of ipt */ 
    235   for (int ii = 0; ii < na; ii++) 
    236     for (int jj = 0; jj < nb; jj++) 
    237       ipt2[jj*na+ii] = ipt[ii*nb+jj]; 
    238   list<Sgm> iscot; 
    239   recense(a, b, ipt, iscot, 0); 
    240   recense(b, a, ipt2, iscot, 1); 
    241  
    242   int nseg = iscot.size(); 
    243   int nc = 0; 
    244   int nc2 = 0; 
    245   while (iscot.size() && nc < 2) 
    246     nc = assemble(iscot, c, d, xc); 
    247   while (iscot.size() && nc2 < 2) 
    248     nc2 = assemble(iscot, c2, d2, xc2); 
    249 //  assert(nseg == nc + nc2 || nseg == 1); // unused segment 
     215        int na = a->n; /* vertices of a */ 
     216        int nb = b->n; /* vertices of b */ 
     217        Coord *c   = new Coord[na+nb]; 
     218        Coord *c2  = new Coord[na+nb]; 
     219        Coord *xc  = new Coord[na+nb]; 
     220        Coord *xc2 = new Coord[na+nb]; 
     221        Coord gc, gc2; 
     222        double *d = new double[na+nb]; 
     223        double *d2 = new double[na+nb]; 
     224        double are, are2; 
     225        Ipt ipt[NMAX*NMAX]; 
     226        Ipt ipt2[NMAX*NMAX]; 
     227        ptsec(a, b, ipt); 
     228        /* make ipt2 transpose of ipt */ 
     229        for (int ii = 0; ii < na; ii++) 
     230                for (int jj = 0; jj < nb; jj++) 
     231                        ipt2[jj*na+ii] = ipt[ii*nb+jj]; 
     232        list<Sgm> iscot; 
     233        recense(a, b, ipt, iscot, 0); 
     234        recense(b, a, ipt2, iscot, 1); 
     235 
     236        int nseg = iscot.size(); 
     237        int nc = 0; 
     238        int nc2 = 0; 
     239        while (iscot.size() && nc < 2) 
     240                nc = assemble(iscot, c, d, xc); 
     241        while (iscot.size() && nc2 < 2) 
     242                nc2 = assemble(iscot, c2, d2, xc2); 
     243//      assert(nseg == nc + nc2 || nseg == 1); // unused segment 
    250244 
    251245        if (!(nseg == nc + nc2 || nseg == 1)) 
     
    266260 
    267261//    intersect_ym(a,b) ; 
    268   if (nc == 1) nc = 0; 
    269   if (nc2 == 1) nc2 = 0; 
    270   gc = barycentre(xc, nc); 
    271   gc2 = barycentre(xc2, nc2); 
    272   orient(nc, xc, c, d, gc); 
    273  
    274   Coord pole = srcGrid.pole; 
    275   if (pole == ORIGIN) pole = tgtGrid.pole; 
    276   const double MINBASE = 1e-11; 
    277   if (nc == 2) /* nc is the number of vertices of super mesh element */ 
    278   { 
    279     double base = arcdist(xc[0], xc[1]); 
     262        if (nc == 1) nc = 0; 
     263        if (nc2 == 1) nc2 = 0; 
     264        gc = barycentre(xc, nc); 
     265        gc2 = barycentre(xc2, nc2); 
     266        orient(nc, xc, c, d, gc); 
     267 
     268        Coord pole = srcGrid.pole; 
     269        if (pole == ORIGIN) pole = tgtGrid.pole; 
     270        const double MINBASE = 1e-11; 
     271        if (nc == 2) /* nc is the number of vertices of super mesh element */ 
     272        { 
     273                double base = arcdist(xc[0], xc[1]); 
    280274cerr << "DID ARRIVE " << base << xc[0] << xc[1] << endl; 
    281     gc = midpoint(gc, midpointSC(xc[0], xc[1])); 
    282     /* intersection area `are` must be zero here unless we have one great and one small circle */ 
    283     are = alun(base, fabs(scalarprod(xc[0], pole))); 
    284   } 
    285   else 
    286   { 
    287     are = airbar(nc, xc, c, d, pole, gc); 
    288   } 
    289   if (nc2 == 2) 
    290   { 
    291     double base = arcdist(xc2[0], xc2[1]); 
     275                gc = midpoint(gc, midpointSC(xc[0], xc[1])); 
     276                /* intersection area `are` must be zero here unless we have one great and one small circle */ 
     277                are = alun(base, fabs(scalarprod(xc[0], pole))); 
     278        } 
     279        else 
     280        { 
     281                are = airbar(nc, xc, c, d, pole, gc); 
     282        } 
     283        if (nc2 == 2) 
     284        { 
     285                double base = arcdist(xc2[0], xc2[1]); 
    292286cerr << "DID ARRIVE " << base << xc2[0] << xc2[1] << endl; 
    293     assert(base > MINBASE); 
    294     gc2 = midpoint(gc2, midpointSC(xc2[0], xc2[1])); 
    295     are2 = alun(base, fabs(scalarprod(xc2[0], pole))); // 0 
    296   } 
    297   else 
    298   { 
    299     are2 = airbar(nc2, xc2, c2, d2, pole, gc2); 
    300   } 
     287                assert(base > MINBASE); 
     288                gc2 = midpoint(gc2, midpointSC(xc2[0], xc2[1])); 
     289                are2 = alun(base, fabs(scalarprod(xc2[0], pole))); // 0 
     290        } 
     291        else 
     292        { 
     293                are2 = airbar(nc2, xc2, c2, d2, pole, gc2); 
     294        } 
    301295 
    302296//  double ym_area=intersect_ym(a,b) ; 
    303297 
    304298  if (nc > 1) 
    305   { 
    306     /* create one super mesh polygon that src and dest point to */ 
    307     Polyg *is = new Polyg; 
    308     is->x = gc; 
    309     is->area = are; 
    310     is->id = b->id; 
    311     is->src_id = b->src_id; 
    312     is->n = nc; 
    313     (a->is).push_back(is); 
    314     (b->is).push_back(is); 
     299        { 
     300                /* create one super mesh polygon that src and dest point to */ 
     301                Polyg *is = new Polyg; 
     302                is->x = gc; 
     303                is->area = are; 
     304                is->id = b->id; 
     305                is->src_id = b->src_id; 
     306                is->n = nc; 
     307                (a->is).push_back(is); 
     308                (b->is).push_back(is); 
    315309/* 
    316     if (  2*fabs(are-ym_area)/(are+ym_area) > 1.1 && ym_area>1e-8) 
     310          if (  2*fabs(are-ym_area)/(are+ym_area) > 1.1 && ym_area>1e-8) 
    317311    { 
    318312      cout<<"Big area difference : "<<are<<"  "<<ym_area<<endl ; 
     
    320314    } 
    321315*/ 
    322 //    cout<<"intersection : "<<are<<" "<< ym_area<<"  diff : "<<fabs(are-ym_area)<<"  ratio : "<<fabs(are-ym_area)/(0.5*(are+ym_area))<<endl ; 
    323   } 
    324   if (nc2 > 1) 
    325   { 
    326     Polyg *is = new Polyg; 
    327     is->x = gc2; 
    328     is->area = are2; 
    329     is->id = b->id; /* intersection holds id of corresponding source element (see Elt class definition for details about id) */ 
    330     is->src_id = b->src_id; 
    331     is->n = nc2; 
    332     (a->is).push_back(is); 
    333     (b->is).push_back(is); 
     316//              cout<<"intersection : "<<are<<" "<< ym_area<<"  diff : "<<fabs(are-ym_area)<<"  ratio : "<<fabs(are-ym_area)/(0.5*(are+ym_area))<<endl ; 
     317        } 
     318        if (nc2 > 1) 
     319        { 
     320                Polyg *is = new Polyg; 
     321                is->x = gc2; 
     322                is->area = are2; 
     323                is->id = b->id; /* intersection holds id of corresponding source element (see Elt class definition for details about id) */ 
     324                is->src_id = b->src_id; 
     325                is->n = nc2; 
     326                (a->is).push_back(is); 
     327                (b->is).push_back(is); 
    334328/* 
    335     if (  2*fabs(are-ym_area)/(are+ym_area) > 1.1 && ym_area>1e-8 ) 
     329    if (        2*fabs(are-ym_area)/(are+ym_area) > 1.1 && ym_area>1e-8 ) 
    336330    { 
    337331      cout<<"Big area difference : "<<are<<"  "<<ym_area<<endl ; 
     
    339333    } 
    340334*/ 
    341 //    cout<<"intersection : "<<are2<<" "<< ym_area<<"  diff : "<<fabs(are-ym_area)<<"  ratio : "<<fabs(are-ym_area)/(0.5*(are+ym_area))<<endl ; 
    342   } 
     335//              cout<<"intersection : "<<are2<<" "<< ym_area<<"  diff : "<<fabs(are-ym_area)<<"  ratio : "<<fabs(are-ym_area)/(0.5*(are+ym_area))<<endl ; 
     336        } 
    343337/* 
    344338  if (nc<=1 && nc2<=1) 
     
    350344  } 
    351345*/ 
    352   delete [] c; 
    353   delete [] c2; 
    354   delete [] xc; 
    355   delete [] xc2; 
    356   delete [] d; 
    357   delete [] d2; 
    358 } 
    359  
    360 } 
     346        delete [] c; 
     347        delete [] c2; 
     348        delete [] xc; 
     349        delete [] xc2; 
     350        delete [] d; 
     351        delete [] d2; 
     352} 
     353 
     354} 
  • XIOS/dev/branch_openmp/extern/remap/src/intersection_ym.cpp

    r1220 r1328  
    1313 
    1414namespace sphereRemap { 
    15  
    16 extern CRemapGrid srcGrid; 
    17 #pragma omp threadprivate(srcGrid) 
    18  
    19 extern CRemapGrid tgtGrid; 
    20 #pragma omp threadprivate(tgtGrid) 
    21  
    2215 
    2316using namespace std; 
  • XIOS/dev/branch_openmp/extern/remap/src/libmapper.cpp

    r1205 r1328  
    1414#include "mapper.hpp" 
    1515#include "cputime.hpp" // cputime 
    16 #include <stdio.h> 
     16 
     17using namespace ep_lib; 
    1718 
    1819using namespace sphereRemap ; 
     
    2122   and deallocated during the second step (computing the weights) */ 
    2223Mapper *mapper; 
    23 #pragma omp threadprivate(mapper) 
     24 
    2425 
    2526/** xxx_bounds_yyy is of length n_vert_per_cell_xxx*n_cell_xxx 
     
    3334                     int order, int* n_weights) 
    3435{ 
    35   assert(src_bounds_lon); 
    36   assert(src_bounds_lat); 
    37   assert(n_vert_per_cell_src >= 3); 
    38   assert(n_cell_src >= 4); 
    39   assert(dst_bounds_lon); 
    40   assert(dst_bounds_lat); 
    41   assert(n_vert_per_cell_dst >= 3); 
    42   assert(n_cell_dst >= 4); 
    43   assert(1 <= order && order <= 2); 
     36        assert(src_bounds_lon); 
     37        assert(src_bounds_lat); 
     38        assert(n_vert_per_cell_src >= 3); 
     39        assert(n_cell_src >= 4); 
     40        assert(dst_bounds_lon); 
     41        assert(dst_bounds_lat); 
     42        assert(n_vert_per_cell_dst >= 3); 
     43        assert(n_cell_dst >= 4); 
     44        assert(1 <= order && order <= 2); 
    4445 
    4546  mapper = new Mapper(MPI_COMM_WORLD); 
     
    8081        double tic = cputime(); 
    8182        mapper = new Mapper(MPI_COMM_WORLD); 
    82         mapper->setVerbosity(PROGRESS) ; 
     83  mapper->setVerbosity(PROGRESS) ; 
    8384        mapper->buildSSTree(src_msh, dst_msh); 
    8485        double tac = cputime(); 
     
    149150        char **argv = NULL; 
    150151        MPI_Init(&argc, &argv);*/ 
    151         //MPI_Init(NULL, NULL); 
    152         int provided; 
    153         MPI_Init_thread(NULL, NULL, 3, &provided); 
    154         assert(provided >= 3); 
     152        MPI_Init(NULL, NULL); 
    155153} 
    156154 
  • XIOS/dev/branch_openmp/extern/remap/src/mapper.cpp

    r1220 r1328  
    1212 
    1313#include "mapper.hpp" 
     14using namespace ep_lib; 
    1415 
    1516namespace sphereRemap { 
    16  
    17 extern CRemapGrid srcGrid; 
    18 #pragma omp threadprivate(srcGrid) 
    19  
    20 extern CRemapGrid tgtGrid; 
    21 #pragma omp threadprivate(tgtGrid) 
    22  
    2317 
    2418/* A subdivition of an array into N sub-arays 
     
    7266void Mapper::setTargetMesh(const double* boundsLon, const double* boundsLat, int nVertex, int nbCells, const double* pole, const long int* globalId) 
    7367{ 
    74     tgtGrid.pole = Coord(pole[0], pole[1], pole[2]); 
    75  
    76     int mpiRank, mpiSize; 
    77     MPI_Comm_rank(communicator, &mpiRank); 
    78     MPI_Comm_size(communicator, &mpiSize); 
    79  
    80     targetElements.reserve(nbCells); 
    81     targetMesh.reserve(nbCells); 
    82  
    83     targetGlobalId.resize(nbCells) ; 
    84     if (globalId==NULL) 
    85     { 
    86         long int offset ; 
    87         long int nb=nbCells ; 
    88         MPI_Scan(&nb,&offset,1,MPI_LONG,MPI_SUM,communicator) ; 
    89         offset=offset-nb ; 
    90         for(int i=0;i<nbCells;i++) targetGlobalId[i]=offset+i ; 
    91     } 
    92     else targetGlobalId.assign(globalId,globalId+nbCells); 
    93  
    94     for (int i = 0; i < nbCells; i++) 
    95     { 
    96         int offs = i*nVertex; 
    97         Elt elt(boundsLon + offs, boundsLat + offs, nVertex); 
    98         targetElements.push_back(elt); 
    99         targetMesh.push_back(Node(elt.x, cptRadius(elt), &sourceElements.back())); 
    100         cptEltGeom(targetElements[i], Coord(pole[0], pole[1], pole[2])); 
    101     } 
     68  tgtGrid.pole = Coord(pole[0], pole[1], pole[2]); 
     69 
     70  int mpiRank, mpiSize; 
     71  MPI_Comm_rank(communicator, &mpiRank); 
     72  MPI_Comm_size(communicator, &mpiSize); 
     73 
     74  targetElements.reserve(nbCells); 
     75  targetMesh.reserve(nbCells); 
     76 
     77  targetGlobalId.resize(nbCells) ; 
     78  if (globalId==NULL) 
     79  { 
     80    long int offset ; 
     81    long int nb=nbCells ; 
     82    MPI_Scan(&nb,&offset,1,MPI_LONG,MPI_SUM,communicator) ; 
     83    offset=offset-nb ; 
     84    for(int i=0;i<nbCells;i++) targetGlobalId[i]=offset+i ; 
     85  } 
     86  else targetGlobalId.assign(globalId,globalId+nbCells); 
     87 
     88  for (int i = 0; i < nbCells; i++) 
     89  { 
     90    int offs = i*nVertex; 
     91    Elt elt(boundsLon + offs, boundsLat + offs, nVertex); 
     92    targetElements.push_back(elt); 
     93    targetMesh.push_back(Node(elt.x, cptRadius(elt), &sourceElements.back())); 
     94    cptEltGeom(targetElements[i], Coord(pole[0], pole[1], pole[2])); 
     95  } 
    10296 
    10397 
     
    106100void Mapper::setSourceValue(const double* val) 
    107101{ 
    108     int size=sourceElements.size() ; 
    109     for(int i=0;i<size;++i) sourceElements[i].val=val[i] ; 
     102  int size=sourceElements.size() ; 
     103  for(int i=0;i<size;++i) sourceElements[i].val=val[i] ; 
    110104} 
    111105 
    112106void Mapper::getTargetValue(double* val) 
    113107{ 
    114     int size=targetElements.size() ; 
    115     for(int i=0;i<size;++i) val[i]=targetElements[i].val ; 
     108  int size=targetElements.size() ; 
     109  for(int i=0;i<size;++i) val[i]=targetElements[i].val ; 
    116110} 
    117111 
    118112vector<double> Mapper::computeWeights(int interpOrder, bool renormalize, bool quantity) 
    119113{ 
    120     vector<double> timings; 
    121     int mpiSize, mpiRank; 
    122     MPI_Comm_size(communicator, &mpiSize); 
    123     MPI_Comm_rank(communicator, &mpiRank); 
    124  
    125     this->buildSSTree(sourceMesh, targetMesh); 
    126  
    127     if (mpiRank == 0 && verbose) cout << "Computing intersections ..." << endl; 
    128     double tic = cputime(); 
    129     computeIntersection(&targetElements[0], targetElements.size()); 
    130     timings.push_back(cputime() - tic); 
    131  
    132     tic = cputime(); 
    133     if (interpOrder == 2) { 
    134         if (mpiRank == 0 && verbose) cout << "Computing grads ..." << endl; 
    135         buildMeshTopology(); 
    136         computeGrads(); 
    137     } 
    138     timings.push_back(cputime() - tic); 
    139  
    140     /* Prepare computation of weights */ 
    141     /* compute number of intersections which for the first order case 
    142        corresponds to the number of edges in the remap matrix */ 
    143     int nIntersections = 0; 
    144     for (int j = 0; j < targetElements.size(); j++) 
    145     { 
    146         Elt &elt = targetElements[j]; 
    147         for (list<Polyg*>::iterator it = elt.is.begin(); it != elt.is.end(); it++) 
    148             nIntersections++; 
    149     } 
    150     /* overallocate for NMAX neighbours for each elements */ 
    151     remapMatrix = new double[nIntersections*NMAX]; 
    152     srcAddress = new int[nIntersections*NMAX]; 
    153     srcRank = new int[nIntersections*NMAX]; 
    154     dstAddress = new int[nIntersections*NMAX]; 
    155     sourceWeightId =new long[nIntersections*NMAX]; 
    156     targetWeightId =new long[nIntersections*NMAX]; 
    157  
    158  
    159     if (mpiRank == 0 && verbose) cout << "Remapping..." << endl; 
    160     tic = cputime(); 
    161     nWeights = remap(&targetElements[0], targetElements.size(), interpOrder, renormalize, quantity); 
    162     timings.push_back(cputime() - tic); 
    163  
    164     for (int i = 0; i < targetElements.size(); i++) targetElements[i].delete_intersections(); 
    165  
    166     return timings; 
     114  vector<double> timings; 
     115  int mpiSize, mpiRank; 
     116  MPI_Comm_size(communicator, &mpiSize); 
     117  MPI_Comm_rank(communicator, &mpiRank); 
     118 
     119 
     120  this->buildSSTree(sourceMesh, targetMesh); 
     121 
     122  if (mpiRank == 0 && verbose) cout << "Computing intersections ..." << endl; 
     123  double tic = cputime(); 
     124  computeIntersection(&targetElements[0], targetElements.size()); 
     125  timings.push_back(cputime() - tic); 
     126 
     127        tic = cputime(); 
     128        if (interpOrder == 2) { 
     129                if (mpiRank == 0 && verbose) cout << "Computing grads ..." << endl; 
     130                buildMeshTopology(); 
     131                computeGrads(); 
     132        } 
     133        timings.push_back(cputime() - tic); 
     134 
     135        /* Prepare computation of weights */ 
     136        /* compute number of intersections which for the first order case 
     137           corresponds to the number of edges in the remap matrix */ 
     138        int nIntersections = 0; 
     139        for (int j = 0; j < targetElements.size(); j++) 
     140        { 
     141                Elt &elt = targetElements[j]; 
     142                for (list<Polyg*>::iterator it = elt.is.begin(); it != elt.is.end(); it++) 
     143                        nIntersections++; 
     144        } 
     145        /* overallocate for NMAX neighbours for each elements */ 
     146        remapMatrix = new double[nIntersections*NMAX]; 
     147        srcAddress = new int[nIntersections*NMAX]; 
     148        srcRank = new int[nIntersections*NMAX]; 
     149        dstAddress = new int[nIntersections*NMAX]; 
     150  sourceWeightId =new long[nIntersections*NMAX]; 
     151  targetWeightId =new long[nIntersections*NMAX]; 
     152 
     153 
     154        if (mpiRank == 0 && verbose) cout << "Remapping..." << endl; 
     155        tic = cputime(); 
     156        nWeights = remap(&targetElements[0], targetElements.size(), interpOrder, renormalize, quantity); 
     157        timings.push_back(cputime() - tic); 
     158 
     159  for (int i = 0; i < targetElements.size(); i++) targetElements[i].delete_intersections(); 
     160 
     161        return timings; 
    167162} 
    168163 
    169164/** 
    170   @param elements are cells of the target grid that are distributed over CPUs 
    171   indepentently of the distribution of the SS-tree. 
    172   @param nbElements is the size of the elements array. 
    173   @param order is the order of interpolaton (must be 1 or 2). 
    174   */ 
     165   @param elements are cells of the target grid that are distributed over CPUs 
     166          indepentently of the distribution of the SS-tree. 
     167   @param nbElements is the size of the elements array. 
     168   @param order is the order of interpolaton (must be 1 or 2). 
     169*/ 
    175170int Mapper::remap(Elt *elements, int nbElements, int order, bool renormalize, bool quantity) 
    176171{ 
    177     int mpiSize, mpiRank; 
    178     MPI_Comm_size(communicator, &mpiSize); 
    179     MPI_Comm_rank(communicator, &mpiRank); 
    180  
    181     /* create list of intersections (super mesh elements) for each rank */ 
    182     multimap<int, Polyg *> *elementList = new multimap<int, Polyg *>[mpiSize]; 
    183     for (int j = 0; j < nbElements; j++) 
    184     { 
    185         Elt& e = elements[j]; 
    186         for (list<Polyg *>::iterator it = e.is.begin(); it != e.is.end(); it++) 
    187             elementList[(*it)->id.rank].insert(pair<int, Polyg *>((*it)->id.ind, *it)); 
    188     } 
    189  
    190     int *nbSendElement = new int[mpiSize]; 
    191     int **sendElement = new int*[mpiSize]; /* indices of elements required from other rank */ 
    192     double **recvValue = new double*[mpiSize]; 
    193     double **recvArea = new double*[mpiSize]; 
    194     Coord **recvGrad = new Coord*[mpiSize]; 
    195     GloId **recvNeighIds = new GloId*[mpiSize]; /* ids of the of the source neighbours which also contribute through gradient */ 
    196     for (int rank = 0; rank < mpiSize; rank++) 
    197     { 
    198         /* get size for allocation */ 
    199         int last = -1; /* compares unequal to any index */ 
    200         int index = -1; /* increased to starting index 0 in first iteration */ 
    201         for (multimap<int, Polyg *>::iterator it = elementList[rank].begin(); it != elementList[rank].end(); ++it) 
    202         { 
    203             if (last != it->first) 
    204                 index++; 
    205             (it->second)->id.ind = index; 
    206             last = it->first; 
    207         } 
    208         nbSendElement[rank] = index + 1; 
    209  
    210         /* if size is non-zero allocate and collect indices of elements on other ranks that we intersect */ 
    211         if (nbSendElement[rank] > 0) 
    212         { 
    213             sendElement[rank] = new int[nbSendElement[rank]]; 
    214             recvValue[rank]   = new double[nbSendElement[rank]]; 
    215             recvArea[rank]    = new double[nbSendElement[rank]]; 
    216             if (order == 2) 
    217             { 
    218                 recvNeighIds[rank] = new GloId[nbSendElement[rank]*(NMAX+1)]; 
    219                 recvGrad[rank]    = new Coord[nbSendElement[rank]*(NMAX+1)]; 
    220             } 
    221             else 
    222                 recvNeighIds[rank] = new GloId[nbSendElement[rank]]; 
    223  
    224             last = -1; 
    225             index = -1; 
    226             for (multimap<int, Polyg *>::iterator it = elementList[rank].begin(); it != elementList[rank].end(); ++it) 
    227             { 
    228                 if (last != it->first) 
    229                     index++; 
    230                 sendElement[rank][index] = it->first; 
    231                 last = it->first; 
    232             } 
    233         } 
    234     } 
    235  
    236     /* communicate sizes of source elements to be sent (index lists and later values and gradients) */ 
    237     int *nbRecvElement = new int[mpiSize]; 
    238     MPI_Alltoall(nbSendElement, 1, MPI_INT, nbRecvElement, 1, MPI_INT, communicator); 
    239  
    240     /* communicate indices of source elements on other ranks whoes value and gradient we need (since intersection) */ 
    241     int nbSendRequest = 0; 
    242     int nbRecvRequest = 0; 
    243     int **recvElement = new int*[mpiSize]; 
    244     double **sendValue = new double*[mpiSize]; 
    245     double **sendArea = new double*[mpiSize]; 
    246     Coord **sendGrad = new Coord*[mpiSize]; 
    247     GloId **sendNeighIds = new GloId*[mpiSize]; 
    248     MPI_Request *sendRequest = new MPI_Request[4*mpiSize]; 
    249     MPI_Request *recvRequest = new MPI_Request[4*mpiSize]; 
    250     for (int rank = 0; rank < mpiSize; rank++) 
    251     { 
    252         if (nbSendElement[rank] > 0) 
    253         { 
    254             MPI_Issend(sendElement[rank], nbSendElement[rank], MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]); 
    255             nbSendRequest++; 
    256         } 
    257  
    258         if (nbRecvElement[rank] > 0) 
    259         { 
    260             recvElement[rank] = new int[nbRecvElement[rank]]; 
    261             sendValue[rank]   = new double[nbRecvElement[rank]]; 
    262             sendArea[rank]   = new double[nbRecvElement[rank]]; 
    263             if (order == 2) 
    264             { 
    265                 sendNeighIds[rank] = new GloId[nbRecvElement[rank]*(NMAX+1)]; 
    266                 sendGrad[rank]    = new Coord[nbRecvElement[rank]*(NMAX+1)]; 
    267             } 
    268             else 
    269             { 
    270                 sendNeighIds[rank] = new GloId[nbRecvElement[rank]]; 
    271             } 
    272             MPI_Irecv(recvElement[rank], nbRecvElement[rank], MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
    273             nbRecvRequest++; 
    274         } 
    275     } 
    276  
    277     MPI_Status *status = new MPI_Status[4*mpiSize]; 
    278  
    279     MPI_Waitall(nbSendRequest, sendRequest, status); 
    280     MPI_Waitall(nbRecvRequest, recvRequest, status); 
    281  
    282     /* for all indices that have been received from requesting ranks: pack values and gradients, then send */ 
    283     nbSendRequest = 0; 
    284     nbRecvRequest = 0; 
    285     for (int rank = 0; rank < mpiSize; rank++) 
    286     { 
    287         if (nbRecvElement[rank] > 0) 
    288         { 
    289             int jj = 0; // jj == j if no weight writing 
    290             for (int j = 0; j < nbRecvElement[rank]; j++) 
    291             { 
    292                 sendValue[rank][j] = sstree.localElements[recvElement[rank][j]].val; 
    293                 sendArea[rank][j] = sstree.localElements[recvElement[rank][j]].area; 
    294                 if (order == 2) 
    295                 { 
    296                     sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].grad; 
    297                     sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].src_id; 
    298                     jj++; 
    299                     for (int i = 0; i < NMAX; i++) 
    300                     { 
    301                         sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].gradNeigh[i]; 
    302                         sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].neighId[i]; 
    303                         jj++; 
    304                     } 
    305                 } 
    306                 else 
    307                     sendNeighIds[rank][j] = sstree.localElements[recvElement[rank][j]].src_id; 
    308             } 
    309             MPI_Issend(sendValue[rank], nbRecvElement[rank], MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]); 
    310             nbSendRequest++; 
    311             MPI_Issend(sendArea[rank], nbRecvElement[rank], MPI_DOUBLE, rank, 1, communicator, &sendRequest[nbSendRequest]); 
    312             nbSendRequest++; 
    313             if (order == 2) 
    314             { 
    315                 MPI_Issend(sendGrad[rank], 3*nbRecvElement[rank]*(NMAX+1), 
    316                         MPI_DOUBLE, rank, 2, communicator, &sendRequest[nbSendRequest]); 
    317                 nbSendRequest++; 
    318                 MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank]*(NMAX+1), MPI_INT, rank, 3, communicator, &sendRequest[nbSendRequest]); 
    319                 //ym  --> attention taille GloId 
    320                 nbSendRequest++; 
    321             } 
    322             else 
    323             { 
    324                 MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank], MPI_INT, rank, 4, communicator, &sendRequest[nbSendRequest]); 
    325                 //ym  --> attention taille GloId 
    326                 nbSendRequest++;                 
    327             } 
    328         } 
    329         if (nbSendElement[rank] > 0) 
    330         { 
    331             MPI_Irecv(recvValue[rank],  nbSendElement[rank], MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
    332             nbRecvRequest++; 
    333             MPI_Irecv(recvArea[rank],  nbSendElement[rank], MPI_DOUBLE, rank, 1, communicator, &recvRequest[nbRecvRequest]); 
    334             nbRecvRequest++; 
    335             if (order == 2) 
    336             { 
    337                 MPI_Irecv(recvGrad[rank], 3*nbSendElement[rank]*(NMAX+1), 
    338                         MPI_DOUBLE, rank, 2, communicator, &recvRequest[nbRecvRequest]); 
    339                 nbRecvRequest++; 
    340                 MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank]*(NMAX+1), MPI_INT, rank, 3, communicator, &recvRequest[nbRecvRequest]); 
    341                 //ym  --> attention taille GloId 
    342                 nbRecvRequest++; 
    343             } 
    344             else 
    345             { 
    346                 MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank], MPI_INT, rank, 4, communicator, &recvRequest[nbRecvRequest]); 
    347                 //ym  --> attention taille GloId 
    348                 nbRecvRequest++; 
    349             } 
    350         } 
    351     } 
    352      
    353     MPI_Waitall(nbRecvRequest, recvRequest, status); 
    354     MPI_Waitall(nbSendRequest, sendRequest, status);  
    355      
    356      
    357  
    358     /* now that all values and gradients are available use them to computed interpolated values on target 
    359        and also to compute weights */ 
    360     int i = 0; 
    361     for (int j = 0; j < nbElements; j++) 
    362     { 
    363         Elt& e = elements[j]; 
    364  
    365         /* since for the 2nd order case source grid elements can contribute to a destination grid element over several "paths" 
    366            (step1: gradient is computed using neighbours on same grid, step2: intersection uses several elements on other grid) 
    367            accumulate them so that there is only one final weight between two elements */ 
    368         map<GloId,double> wgt_map; 
    369  
    370         /* for destination element `e` loop over all intersetions/the corresponding source elements */ 
    371         for (list<Polyg *>::iterator it = e.is.begin(); it != e.is.end(); it++) 
    372         { 
    373             /* it is the intersection element, so it->x and it->area are barycentre and area of intersection element (super mesh) 
    374                but it->id is id of the source element that it intersects */ 
    375             int n1 = (*it)->id.ind; 
    376             int rank = (*it)->id.rank; 
    377             double fk = recvValue[rank][n1]; 
    378             double srcArea = recvArea[rank][n1]; 
    379             double w = (*it)->area; 
    380             if (quantity) w/=srcArea ; 
    381  
    382             /* first order: src value times weight (weight = supermesh area), later divide by target area */ 
    383             int kk = (order == 2) ? n1 * (NMAX + 1) : n1; 
    384             GloId neighID = recvNeighIds[rank][kk]; 
    385             wgt_map[neighID] += w; 
    386  
    387             if (order == 2) 
    388             { 
    389                 for (int k = 0; k < NMAX+1; k++) 
    390                 { 
    391                     int kk = n1 * (NMAX + 1) + k; 
    392                     GloId neighID = recvNeighIds[rank][kk]; 
    393                     if (neighID.ind != -1)  wgt_map[neighID] += w * scalarprod(recvGrad[rank][kk], (*it)->x); 
    394                 } 
    395  
    396             } 
    397         } 
    398  
    399         double renorm=0; 
    400         if (renormalize)  
    401             for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++) renorm+=it->second / e.area; 
    402         else renorm=1. ; 
    403  
    404         for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++) 
    405         { 
    406             if (quantity)  this->remapMatrix[i] = (it->second ) / renorm; 
    407             else this->remapMatrix[i] = (it->second / e.area) / renorm; 
    408             this->srcAddress[i] = it->first.ind; 
    409             this->srcRank[i] = it->first.rank; 
    410             this->dstAddress[i] = j; 
    411             this->sourceWeightId[i]= it->first.globalId ; 
    412             this->targetWeightId[i]= targetGlobalId[j] ; 
    413             i++; 
    414         } 
    415     } 
    416      
    417     /* free all memory allocated in this function */ 
    418     for (int rank = 0; rank < mpiSize; rank++) 
    419     { 
    420         if (nbSendElement[rank] > 0) 
    421         { 
    422             delete[] sendElement[rank]; 
    423             delete[] recvValue[rank]; 
    424             delete[] recvArea[rank]; 
    425             if (order == 2) 
    426             { 
    427                 delete[] recvGrad[rank]; 
    428             } 
    429             delete[] recvNeighIds[rank]; 
    430         } 
    431         if (nbRecvElement[rank] > 0) 
    432         { 
    433             delete[] recvElement[rank]; 
    434             delete[] sendValue[rank]; 
    435             delete[] sendArea[rank]; 
    436             if (order == 2) 
    437                 delete[] sendGrad[rank]; 
    438             delete[] sendNeighIds[rank]; 
    439         } 
    440     } 
    441     delete[] status; 
    442     delete[] sendRequest; 
    443     delete[] recvRequest; 
    444     delete[] elementList; 
    445     delete[] nbSendElement; 
    446     delete[] nbRecvElement; 
    447     delete[] sendElement; 
    448     delete[] recvElement; 
    449     delete[] sendValue; 
    450     delete[] recvValue; 
    451     delete[] sendGrad; 
    452     delete[] recvGrad; 
    453     delete[] sendNeighIds; 
    454     delete[] recvNeighIds; 
    455      
    456     return i; 
     172        int mpiSize, mpiRank; 
     173        MPI_Comm_size(communicator, &mpiSize); 
     174        MPI_Comm_rank(communicator, &mpiRank); 
     175 
     176        /* create list of intersections (super mesh elements) for each rank */ 
     177        multimap<int, Polyg *> *elementList = new multimap<int, Polyg *>[mpiSize]; 
     178        for (int j = 0; j < nbElements; j++) 
     179        { 
     180                Elt& e = elements[j]; 
     181                for (list<Polyg *>::iterator it = e.is.begin(); it != e.is.end(); it++) 
     182                        elementList[(*it)->id.rank].insert(pair<int, Polyg *>((*it)->id.ind, *it)); 
     183        } 
     184 
     185        int *nbSendElement = new int[mpiSize]; 
     186        int **sendElement = new int*[mpiSize]; /* indices of elements required from other rank */ 
     187        double **recvValue = new double*[mpiSize]; 
     188        double **recvArea = new double*[mpiSize]; 
     189        Coord **recvGrad = new Coord*[mpiSize]; 
     190        GloId **recvNeighIds = new GloId*[mpiSize]; /* ids of the of the source neighbours which also contribute through gradient */ 
     191        for (int rank = 0; rank < mpiSize; rank++) 
     192        { 
     193                /* get size for allocation */ 
     194                int last = -1; /* compares unequal to any index */ 
     195                int index = -1; /* increased to starting index 0 in first iteration */ 
     196                for (multimap<int, Polyg *>::iterator it = elementList[rank].begin(); it != elementList[rank].end(); ++it) 
     197                { 
     198                        if (last != it->first) 
     199                                index++; 
     200                        (it->second)->id.ind = index; 
     201                        last = it->first; 
     202                } 
     203                nbSendElement[rank] = index + 1; 
     204 
     205                /* if size is non-zero allocate and collect indices of elements on other ranks that we intersect */ 
     206                if (nbSendElement[rank] > 0) 
     207                { 
     208                        sendElement[rank] = new int[nbSendElement[rank]]; 
     209                        recvValue[rank]   = new double[nbSendElement[rank]]; 
     210                        recvArea[rank]    = new double[nbSendElement[rank]]; 
     211                        if (order == 2) 
     212                        { 
     213                                recvNeighIds[rank] = new GloId[nbSendElement[rank]*(NMAX+1)]; 
     214                                recvGrad[rank]    = new Coord[nbSendElement[rank]*(NMAX+1)]; 
     215                        } 
     216                        else 
     217                                recvNeighIds[rank] = new GloId[nbSendElement[rank]]; 
     218 
     219                        last = -1; 
     220                        index = -1; 
     221                        for (multimap<int, Polyg *>::iterator it = elementList[rank].begin(); it != elementList[rank].end(); ++it) 
     222                        { 
     223                                if (last != it->first) 
     224                                        index++; 
     225                                sendElement[rank][index] = it->first; 
     226                                last = it->first; 
     227                        } 
     228                } 
     229        } 
     230 
     231        /* communicate sizes of source elements to be sent (index lists and later values and gradients) */ 
     232        int *nbRecvElement = new int[mpiSize]; 
     233        MPI_Alltoall(nbSendElement, 1, MPI_INT, nbRecvElement, 1, MPI_INT, communicator); 
     234 
     235        /* communicate indices of source elements on other ranks whoes value and gradient we need (since intersection) */ 
     236        int nbSendRequest = 0; 
     237        int nbRecvRequest = 0; 
     238        int **recvElement = new int*[mpiSize]; 
     239        double **sendValue = new double*[mpiSize]; 
     240        double **sendArea = new double*[mpiSize]; 
     241        Coord **sendGrad = new Coord*[mpiSize]; 
     242        GloId **sendNeighIds = new GloId*[mpiSize]; 
     243        MPI_Request *sendRequest = new MPI_Request[4*mpiSize]; 
     244        MPI_Request *recvRequest = new MPI_Request[4*mpiSize]; 
     245        for (int rank = 0; rank < mpiSize; rank++) 
     246        { 
     247                if (nbSendElement[rank] > 0) 
     248                { 
     249                        MPI_Issend(sendElement[rank], nbSendElement[rank], MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]); 
     250                        nbSendRequest++; 
     251                } 
     252 
     253                if (nbRecvElement[rank] > 0) 
     254                { 
     255                        recvElement[rank] = new int[nbRecvElement[rank]]; 
     256                        sendValue[rank]   = new double[nbRecvElement[rank]]; 
     257                        sendArea[rank]   = new double[nbRecvElement[rank]]; 
     258                        if (order == 2) 
     259                        { 
     260                                sendNeighIds[rank] = new GloId[nbRecvElement[rank]*(NMAX+1)]; 
     261                                sendGrad[rank]    = new Coord[nbRecvElement[rank]*(NMAX+1)]; 
     262                        } 
     263                        else 
     264                        { 
     265                                sendNeighIds[rank] = new GloId[nbRecvElement[rank]]; 
     266                        } 
     267                        MPI_Irecv(recvElement[rank], nbRecvElement[rank], MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
     268                        nbRecvRequest++; 
     269                } 
     270        } 
     271 
     272  MPI_Status *status = new MPI_Status[4*mpiSize]; 
     273 
     274  MPI_Waitall(nbSendRequest, sendRequest, status); 
     275  MPI_Waitall(nbRecvRequest, recvRequest, status); 
     276 
     277        /* for all indices that have been received from requesting ranks: pack values and gradients, then send */ 
     278        nbSendRequest = 0; 
     279        nbRecvRequest = 0; 
     280        for (int rank = 0; rank < mpiSize; rank++) 
     281        { 
     282                if (nbRecvElement[rank] > 0) 
     283                { 
     284                        int jj = 0; // jj == j if no weight writing 
     285                        for (int j = 0; j < nbRecvElement[rank]; j++) 
     286                        { 
     287                                sendValue[rank][j] = sstree.localElements[recvElement[rank][j]].val; 
     288                                sendArea[rank][j] = sstree.localElements[recvElement[rank][j]].area; 
     289                                if (order == 2) 
     290                                { 
     291                                        sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].grad; 
     292                                        sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].src_id; 
     293                                        jj++; 
     294                                        for (int i = 0; i < NMAX; i++) 
     295                                        { 
     296                                                sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].gradNeigh[i]; 
     297            sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].neighId[i]; 
     298                                                jj++; 
     299                                        } 
     300                                } 
     301                                else 
     302                                        sendNeighIds[rank][j] = sstree.localElements[recvElement[rank][j]].src_id; 
     303                        } 
     304                        MPI_Issend(sendValue[rank],  nbRecvElement[rank], MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]); 
     305                        nbSendRequest++; 
     306                        MPI_Issend(sendArea[rank],  nbRecvElement[rank], MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]); 
     307                        nbSendRequest++; 
     308                        if (order == 2) 
     309                        { 
     310                                MPI_Issend(sendGrad[rank], 3*nbRecvElement[rank]*(NMAX+1), 
     311                                                                MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]); 
     312                                nbSendRequest++; 
     313                                MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank]*(NMAX+1), MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]); 
     314//ym  --> attention taille GloId 
     315                                nbSendRequest++; 
     316                        } 
     317                        else 
     318                        { 
     319                                MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank], MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]); 
     320//ym  --> attention taille GloId 
     321                                nbSendRequest++; 
     322                        } 
     323                } 
     324                if (nbSendElement[rank] > 0) 
     325                { 
     326                        MPI_Irecv(recvValue[rank],  nbSendElement[rank], MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
     327                        nbRecvRequest++; 
     328                        MPI_Irecv(recvArea[rank],  nbSendElement[rank], MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
     329                        nbRecvRequest++; 
     330                        if (order == 2) 
     331                        { 
     332                                MPI_Irecv(recvGrad[rank], 3*nbSendElement[rank]*(NMAX+1), 
     333                                                MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
     334                                nbRecvRequest++; 
     335                                MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank]*(NMAX+1), MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
     336//ym  --> attention taille GloId 
     337                                nbRecvRequest++; 
     338                        } 
     339                        else 
     340                        { 
     341                                MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank], MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
     342//ym  --> attention taille GloId 
     343                                nbRecvRequest++; 
     344                        } 
     345                } 
     346        } 
     347         
     348        MPI_Waitall(nbSendRequest, sendRequest, status); 
     349        MPI_Waitall(nbRecvRequest, recvRequest, status); 
     350         
     351 
     352        /* now that all values and gradients are available use them to computed interpolated values on target 
     353           and also to compute weights */ 
     354        int i = 0; 
     355        for (int j = 0; j < nbElements; j++) 
     356        { 
     357                Elt& e = elements[j]; 
     358 
     359                /* since for the 2nd order case source grid elements can contribute to a destination grid element over several "paths" 
     360                   (step1: gradient is computed using neighbours on same grid, step2: intersection uses several elements on other grid) 
     361                   accumulate them so that there is only one final weight between two elements */ 
     362                map<GloId,double> wgt_map; 
     363 
     364                /* for destination element `e` loop over all intersetions/the corresponding source elements */ 
     365                for (list<Polyg *>::iterator it = e.is.begin(); it != e.is.end(); it++) 
     366                { 
     367                        /* it is the intersection element, so it->x and it->area are barycentre and area of intersection element (super mesh) 
     368                        but it->id is id of the source element that it intersects */ 
     369                        int n1 = (*it)->id.ind; 
     370                        int rank = (*it)->id.rank; 
     371                        double fk = recvValue[rank][n1]; 
     372                        double srcArea = recvArea[rank][n1]; 
     373                        double w = (*it)->area; 
     374      if (quantity) w/=srcArea ; 
     375 
     376                        /* first order: src value times weight (weight = supermesh area), later divide by target area */ 
     377                        int kk = (order == 2) ? n1 * (NMAX + 1) : n1; 
     378                        GloId neighID = recvNeighIds[rank][kk]; 
     379                        wgt_map[neighID] += w; 
     380 
     381                        if (order == 2) 
     382                        { 
     383                                for (int k = 0; k < NMAX+1; k++) 
     384                                { 
     385                                        int kk = n1 * (NMAX + 1) + k; 
     386                                        GloId neighID = recvNeighIds[rank][kk]; 
     387                                        if (neighID.ind != -1)  wgt_map[neighID] += w * scalarprod(recvGrad[rank][kk], (*it)->x); 
     388                                } 
     389 
     390                        } 
     391                } 
     392 
     393    double renorm=0; 
     394    if (renormalize)  
     395      for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++) renorm+=it->second / e.area; 
     396    else renorm=1. ; 
     397 
     398    for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++) 
     399                { 
     400      if (quantity)  this->remapMatrix[i] = (it->second ) / renorm; 
     401                        else this->remapMatrix[i] = (it->second / e.area) / renorm; 
     402                        this->srcAddress[i] = it->first.ind; 
     403                        this->srcRank[i] = it->first.rank; 
     404                        this->dstAddress[i] = j; 
     405      this->sourceWeightId[i]= it->first.globalId ; 
     406      this->targetWeightId[i]= targetGlobalId[j] ; 
     407                        i++; 
     408                } 
     409        } 
     410 
     411        /* free all memory allocated in this function */ 
     412        for (int rank = 0; rank < mpiSize; rank++) 
     413        { 
     414                if (nbSendElement[rank] > 0) 
     415                { 
     416                        delete[] sendElement[rank]; 
     417                        delete[] recvValue[rank]; 
     418                        delete[] recvArea[rank]; 
     419                        if (order == 2) 
     420                        { 
     421                                delete[] recvGrad[rank]; 
     422                        } 
     423                        delete[] recvNeighIds[rank]; 
     424                } 
     425                if (nbRecvElement[rank] > 0) 
     426                { 
     427                        delete[] recvElement[rank]; 
     428                        delete[] sendValue[rank]; 
     429                        delete[] sendArea[rank]; 
     430                        if (order == 2) 
     431                                delete[] sendGrad[rank]; 
     432                        delete[] sendNeighIds[rank]; 
     433                } 
     434        } 
     435        delete[] status; 
     436        delete[] sendRequest; 
     437        delete[] recvRequest; 
     438        delete[] elementList; 
     439        delete[] nbSendElement; 
     440        delete[] nbRecvElement; 
     441        delete[] sendElement; 
     442        delete[] recvElement; 
     443        delete[] sendValue; 
     444        delete[] recvValue; 
     445        delete[] sendGrad; 
     446        delete[] recvGrad; 
     447        delete[] sendNeighIds; 
     448        delete[] recvNeighIds; 
     449        return i; 
    457450} 
    458451 
    459452void Mapper::computeGrads() 
    460453{ 
    461     /* array of pointers to collect local elements and elements received from other cpu */ 
    462     vector<Elt*> globalElements(sstree.nbLocalElements + nbNeighbourElements); 
    463     int index = 0; 
    464     for (int i = 0; i < sstree.nbLocalElements; i++, index++) 
    465         globalElements[index] = &(sstree.localElements[i]); 
    466     for (int i = 0; i < nbNeighbourElements; i++, index++) 
    467         globalElements[index] = &neighbourElements[i]; 
    468  
    469     update_baryc(sstree.localElements, sstree.nbLocalElements); 
    470     computeGradients(&globalElements[0], sstree.nbLocalElements); 
     454        /* array of pointers to collect local elements and elements received from other cpu */ 
     455        vector<Elt*> globalElements(sstree.nbLocalElements + nbNeighbourElements); 
     456        int index = 0; 
     457        for (int i = 0; i < sstree.nbLocalElements; i++, index++) 
     458                globalElements[index] = &(sstree.localElements[i]); 
     459        for (int i = 0; i < nbNeighbourElements; i++, index++) 
     460                globalElements[index] = &neighbourElements[i]; 
     461 
     462        update_baryc(sstree.localElements, sstree.nbLocalElements); 
     463        computeGradients(&globalElements[0], sstree.nbLocalElements); 
    471464} 
    472465 
    473466/** for each element of the source grid, finds all the neighbouring elements that share an edge 
    474   (filling array neighbourElements). This is used later to compute gradients */ 
     467    (filling array neighbourElements). This is used later to compute gradients */ 
    475468void Mapper::buildMeshTopology() 
    476469{ 
    477     int mpiSize, mpiRank; 
    478     MPI_Comm_size(communicator, &mpiSize); 
    479     MPI_Comm_rank(communicator, &mpiRank); 
    480  
    481     vector<Node> *routingList = new vector<Node>[mpiSize]; 
    482     vector<vector<int> > routes(sstree.localTree.leafs.size()); 
    483  
    484     sstree.routeIntersections(routes, sstree.localTree.leafs); 
    485  
    486     for (int i = 0; i < routes.size(); ++i) 
    487         for (int k = 0; k < routes[i].size(); ++k) 
    488             routingList[routes[i][k]].push_back(sstree.localTree.leafs[i]); 
    489     routingList[mpiRank].clear(); 
    490  
    491  
    492     CMPIRouting mpiRoute(communicator); 
    493     mpiRoute.init(routes); 
    494     int nRecv = mpiRoute.getTotalSourceElement(); 
    495  
    496     int *nbSendNode = new int[mpiSize]; 
    497     int *nbRecvNode = new int[mpiSize]; 
    498     int *sendMessageSize = new int[mpiSize]; 
    499     int *recvMessageSize = new int[mpiSize]; 
    500  
    501     for (int rank = 0; rank < mpiSize; rank++) 
    502     { 
    503         nbSendNode[rank] = routingList[rank].size(); 
    504         sendMessageSize[rank] = 0; 
    505         for (size_t j = 0; j < routingList[rank].size(); j++) 
    506         { 
    507             Elt *elt = (Elt *) (routingList[rank][j].data); 
    508             sendMessageSize[rank] += packedPolygonSize(*elt); 
    509         } 
    510     } 
    511  
    512     MPI_Alltoall(nbSendNode, 1, MPI_INT, nbRecvNode, 1, MPI_INT, communicator); 
    513     MPI_Alltoall(sendMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator); 
    514  
    515     char **sendBuffer = new char*[mpiSize]; 
    516     char **recvBuffer = new char*[mpiSize]; 
    517     int *pos = new int[mpiSize]; 
    518  
    519     for (int rank = 0; rank < mpiSize; rank++) 
    520     { 
    521         if (nbSendNode[rank] > 0) sendBuffer[rank] = new char[sendMessageSize[rank]]; 
    522         if (nbRecvNode[rank] > 0) recvBuffer[rank] = new char[recvMessageSize[rank]]; 
    523     } 
    524  
    525     for (int rank = 0; rank < mpiSize; rank++) 
    526     { 
    527         pos[rank] = 0; 
    528         for (size_t j = 0; j < routingList[rank].size(); j++) 
    529         { 
    530             Elt *elt = (Elt *) (routingList[rank][j].data); 
    531             packPolygon(*elt, sendBuffer[rank], pos[rank]); 
    532         } 
    533     } 
    534     delete [] routingList; 
    535  
    536  
    537     int nbSendRequest = 0; 
    538     int nbRecvRequest = 0; 
    539     MPI_Request *sendRequest = new MPI_Request[mpiSize]; 
    540     MPI_Request *recvRequest = new MPI_Request[mpiSize]; 
    541     MPI_Status  *status      = new MPI_Status[mpiSize]; 
    542  
    543     for (int rank = 0; rank < mpiSize; rank++) 
    544     { 
    545         if (nbSendNode[rank] > 0) 
    546         { 
    547             MPI_Issend(sendBuffer[rank], sendMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]); 
    548             nbSendRequest++; 
    549         } 
    550         if (nbRecvNode[rank] > 0) 
    551         { 
    552             MPI_Irecv(recvBuffer[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
    553             nbRecvRequest++; 
    554         } 
    555     } 
    556      
    557     MPI_Waitall(nbRecvRequest, recvRequest, status); 
    558     MPI_Waitall(nbSendRequest, sendRequest, status); 
    559  
    560     for (int rank = 0; rank < mpiSize; rank++) 
    561         if (nbSendNode[rank] > 0) delete [] sendBuffer[rank]; 
    562     delete [] sendBuffer; 
    563  
    564     char **sendBuffer2 = new char*[mpiSize]; 
    565     char **recvBuffer2 = new char*[mpiSize]; 
    566  
    567     for (int rank = 0; rank < mpiSize; rank++) 
    568     { 
    569         nbSendNode[rank] = 0; 
    570         sendMessageSize[rank] = 0; 
    571  
    572         if (nbRecvNode[rank] > 0) 
    573         { 
    574             set<NodePtr> neighbourList; 
    575             pos[rank] = 0; 
    576             for (int j = 0; j < nbRecvNode[rank]; j++) 
    577             { 
    578                 Elt elt; 
    579                 unpackPolygon(elt, recvBuffer[rank], pos[rank]); 
    580                 Node node(elt.x, cptRadius(elt), &elt); 
    581                 findNeighbour(sstree.localTree.root, &node, neighbourList); 
    582             } 
    583             nbSendNode[rank] = neighbourList.size(); 
    584             for (set<NodePtr>::iterator it = neighbourList.begin(); it != neighbourList.end(); it++) 
    585             { 
    586                 Elt *elt = (Elt *) ((*it)->data); 
    587                 sendMessageSize[rank] += packedPolygonSize(*elt); 
    588             } 
    589  
    590             sendBuffer2[rank] = new char[sendMessageSize[rank]]; 
    591             pos[rank] = 0; 
    592  
    593             for (set<NodePtr>::iterator it = neighbourList.begin(); it != neighbourList.end(); it++) 
    594             { 
    595                 Elt *elt = (Elt *) ((*it)->data); 
    596                 packPolygon(*elt, sendBuffer2[rank], pos[rank]); 
    597             } 
    598         } 
    599     } 
    600     for (int rank = 0; rank < mpiSize; rank++) 
    601         if (nbRecvNode[rank] > 0) delete [] recvBuffer[rank]; 
    602     delete [] recvBuffer; 
    603  
    604  
    605     MPI_Barrier(communicator); 
    606     MPI_Alltoall(nbSendNode, 1, MPI_INT, nbRecvNode, 1, MPI_INT, communicator); 
    607     MPI_Alltoall(sendMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator); 
    608  
    609     for (int rank = 0; rank < mpiSize; rank++) 
    610         if (nbRecvNode[rank] > 0) recvBuffer2[rank] = new char[recvMessageSize[rank]]; 
    611  
    612     nbSendRequest = 0; 
    613     nbRecvRequest = 0; 
    614  
    615     for (int rank = 0; rank < mpiSize; rank++) 
    616     { 
    617         if (nbSendNode[rank] > 0) 
    618         { 
    619             MPI_Issend(sendBuffer2[rank], sendMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]); 
    620             nbSendRequest++; 
    621         } 
    622         if (nbRecvNode[rank] > 0) 
    623         { 
    624             MPI_Irecv(recvBuffer2[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
    625             nbRecvRequest++; 
    626         } 
    627     } 
    628  
    629     MPI_Waitall(nbRecvRequest, recvRequest, status); 
    630     MPI_Waitall(nbSendRequest, sendRequest, status); 
    631  
    632     int nbNeighbourNodes = 0; 
    633     for (int rank = 0; rank < mpiSize; rank++) 
    634         nbNeighbourNodes += nbRecvNode[rank]; 
    635  
    636     neighbourElements = new Elt[nbNeighbourNodes]; 
    637     nbNeighbourElements = nbNeighbourNodes; 
    638  
    639     int index = 0; 
    640     for (int rank = 0; rank < mpiSize; rank++) 
    641     { 
    642         pos[rank] = 0; 
    643         for (int j = 0; j < nbRecvNode[rank]; j++) 
    644         { 
    645             unpackPolygon(neighbourElements[index], recvBuffer2[rank], pos[rank]); 
    646             neighbourElements[index].id.ind = sstree.localTree.leafs.size() + index; 
    647             index++; 
    648         } 
    649     } 
    650     for (int rank = 0; rank < mpiSize; rank++) 
    651     { 
    652         if (nbRecvNode[rank] > 0) delete [] recvBuffer2[rank]; 
    653         if (nbSendNode[rank] > 0) delete [] sendBuffer2[rank]; 
    654     } 
    655     delete [] recvBuffer2; 
    656     delete [] sendBuffer2; 
    657     delete [] sendMessageSize; 
    658     delete [] recvMessageSize; 
    659     delete [] nbSendNode; 
    660     delete [] nbRecvNode; 
    661     delete [] sendRequest; 
    662     delete [] recvRequest; 
    663     delete [] status; 
    664     delete [] pos; 
    665  
    666     /* re-compute on received elements to avoid having to send this information */ 
    667     neighbourNodes.resize(nbNeighbourNodes); 
    668     setCirclesAndLinks(neighbourElements, neighbourNodes); 
    669     cptAllEltsGeom(neighbourElements, nbNeighbourNodes, srcGrid.pole); 
    670  
    671     /* the local SS tree must include nodes from other cpus if they are potential 
    672        intersector of a local node */ 
    673     sstree.localTree.insertNodes(neighbourNodes); 
    674  
    675     /* for every local element, 
    676        use the SS-tree to find all elements (including neighbourElements) 
    677        who are potential neighbours because their circles intersect, 
    678        then check all canditates for common edges to build up connectivity information 
    679        */ 
    680     for (int j = 0; j < sstree.localTree.leafs.size(); j++) 
    681     { 
    682         Node& node = sstree.localTree.leafs[j]; 
    683  
    684         /* find all leafs whoes circles that intersect node's circle and save into node->intersectors */ 
    685         node.search(sstree.localTree.root); 
    686  
    687         Elt *elt = (Elt *)(node.data); 
    688  
    689         for (int i = 0; i < elt->n; i++) elt->neighbour[i] = NOT_FOUND; 
    690  
    691         /* for element `elt` loop through all nodes in the SS-tree 
    692            whoes circles intersect with the circle around `elt` (the SS intersectors) 
    693            and check if they are neighbours in the sense that the two elements share an edge. 
    694            If they do, save this information for elt */ 
    695         for (list<NodePtr>::iterator it = (node.intersectors).begin(); it != (node.intersectors).end(); ++it) 
    696         { 
    697             Elt *elt2 = (Elt *)((*it)->data); 
    698             set_neighbour(*elt, *elt2); 
    699         } 
    700  
    701     } 
     470        int mpiSize, mpiRank; 
     471        MPI_Comm_size(communicator, &mpiSize); 
     472        MPI_Comm_rank(communicator, &mpiRank); 
     473 
     474        vector<Node> *routingList = new vector<Node>[mpiSize]; 
     475        vector<vector<int> > routes(sstree.localTree.leafs.size()); 
     476 
     477        sstree.routeIntersections(routes, sstree.localTree.leafs); 
     478 
     479        for (int i = 0; i < routes.size(); ++i) 
     480                for (int k = 0; k < routes[i].size(); ++k) 
     481                        routingList[routes[i][k]].push_back(sstree.localTree.leafs[i]); 
     482        routingList[mpiRank].clear(); 
     483 
     484 
     485        CMPIRouting mpiRoute(communicator); 
     486        mpiRoute.init(routes); 
     487        int nRecv = mpiRoute.getTotalSourceElement(); 
     488 
     489        int *nbSendNode = new int[mpiSize]; 
     490        int *nbRecvNode = new int[mpiSize]; 
     491        int *sendMessageSize = new int[mpiSize]; 
     492        int *recvMessageSize = new int[mpiSize]; 
     493 
     494        for (int rank = 0; rank < mpiSize; rank++) 
     495        { 
     496                nbSendNode[rank] = routingList[rank].size(); 
     497                sendMessageSize[rank] = 0; 
     498                for (size_t j = 0; j < routingList[rank].size(); j++) 
     499                { 
     500                        Elt *elt = (Elt *) (routingList[rank][j].data); 
     501                        sendMessageSize[rank] += packedPolygonSize(*elt); 
     502                } 
     503        } 
     504 
     505        MPI_Alltoall(nbSendNode, 1, MPI_INT, nbRecvNode, 1, MPI_INT, communicator); 
     506        MPI_Alltoall(sendMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator); 
     507 
     508        char **sendBuffer = new char*[mpiSize]; 
     509        char **recvBuffer = new char*[mpiSize]; 
     510        int *pos = new int[mpiSize]; 
     511 
     512        for (int rank = 0; rank < mpiSize; rank++) 
     513        { 
     514                if (nbSendNode[rank] > 0) sendBuffer[rank] = new char[sendMessageSize[rank]]; 
     515                if (nbRecvNode[rank] > 0) recvBuffer[rank] = new char[recvMessageSize[rank]]; 
     516        } 
     517 
     518        for (int rank = 0; rank < mpiSize; rank++) 
     519        { 
     520                pos[rank] = 0; 
     521                for (size_t j = 0; j < routingList[rank].size(); j++) 
     522                { 
     523                        Elt *elt = (Elt *) (routingList[rank][j].data); 
     524                        packPolygon(*elt, sendBuffer[rank], pos[rank]); 
     525                } 
     526        } 
     527        delete [] routingList; 
     528 
     529 
     530        int nbSendRequest = 0; 
     531        int nbRecvRequest = 0; 
     532        MPI_Request *sendRequest = new MPI_Request[mpiSize]; 
     533        MPI_Request *recvRequest = new MPI_Request[mpiSize]; 
     534        MPI_Status  *status      = new MPI_Status[mpiSize]; 
     535 
     536        for (int rank = 0; rank < mpiSize; rank++) 
     537        { 
     538                if (nbSendNode[rank] > 0) 
     539                { 
     540                        MPI_Issend(sendBuffer[rank], sendMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]); 
     541                        nbSendRequest++; 
     542                } 
     543                if (nbRecvNode[rank] > 0) 
     544                { 
     545                        MPI_Irecv(recvBuffer[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
     546                        nbRecvRequest++; 
     547                } 
     548        } 
     549 
     550  MPI_Waitall(nbRecvRequest, recvRequest, status); 
     551  MPI_Waitall(nbSendRequest, sendRequest, status); 
     552 
     553        for (int rank = 0; rank < mpiSize; rank++) 
     554                if (nbSendNode[rank] > 0) delete [] sendBuffer[rank]; 
     555        delete [] sendBuffer; 
     556 
     557        char **sendBuffer2 = new char*[mpiSize]; 
     558        char **recvBuffer2 = new char*[mpiSize]; 
     559 
     560        for (int rank = 0; rank < mpiSize; rank++) 
     561        { 
     562                nbSendNode[rank] = 0; 
     563                sendMessageSize[rank] = 0; 
     564 
     565                if (nbRecvNode[rank] > 0) 
     566                { 
     567                        set<NodePtr> neighbourList; 
     568                        pos[rank] = 0; 
     569                        for (int j = 0; j < nbRecvNode[rank]; j++) 
     570                        { 
     571                                Elt elt; 
     572                                unpackPolygon(elt, recvBuffer[rank], pos[rank]); 
     573                                Node node(elt.x, cptRadius(elt), &elt); 
     574                                findNeighbour(sstree.localTree.root, &node, neighbourList); 
     575                        } 
     576                        nbSendNode[rank] = neighbourList.size(); 
     577                        for (set<NodePtr>::iterator it = neighbourList.begin(); it != neighbourList.end(); it++) 
     578                        { 
     579                                Elt *elt = (Elt *) ((*it)->data); 
     580                                sendMessageSize[rank] += packedPolygonSize(*elt); 
     581                        } 
     582 
     583                        sendBuffer2[rank] = new char[sendMessageSize[rank]]; 
     584                        pos[rank] = 0; 
     585 
     586                        for (set<NodePtr>::iterator it = neighbourList.begin(); it != neighbourList.end(); it++) 
     587                        { 
     588                                Elt *elt = (Elt *) ((*it)->data); 
     589                                packPolygon(*elt, sendBuffer2[rank], pos[rank]); 
     590                        } 
     591                } 
     592        } 
     593        for (int rank = 0; rank < mpiSize; rank++) 
     594                if (nbRecvNode[rank] > 0) delete [] recvBuffer[rank]; 
     595        delete [] recvBuffer; 
     596 
     597 
     598        MPI_Barrier(communicator); 
     599        MPI_Alltoall(nbSendNode, 1, MPI_INT, nbRecvNode, 1, MPI_INT, communicator); 
     600        MPI_Alltoall(sendMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator); 
     601 
     602        for (int rank = 0; rank < mpiSize; rank++) 
     603                if (nbRecvNode[rank] > 0) recvBuffer2[rank] = new char[recvMessageSize[rank]]; 
     604 
     605        nbSendRequest = 0; 
     606        nbRecvRequest = 0; 
     607 
     608        for (int rank = 0; rank < mpiSize; rank++) 
     609        { 
     610                if (nbSendNode[rank] > 0) 
     611                { 
     612                        MPI_Issend(sendBuffer2[rank], sendMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]); 
     613                        nbSendRequest++; 
     614                } 
     615                if (nbRecvNode[rank] > 0) 
     616                { 
     617                        MPI_Irecv(recvBuffer2[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
     618                        nbRecvRequest++; 
     619                } 
     620        } 
     621 
     622        MPI_Waitall(nbRecvRequest, recvRequest, status); 
     623        MPI_Waitall(nbSendRequest, sendRequest, status); 
     624 
     625        int nbNeighbourNodes = 0; 
     626        for (int rank = 0; rank < mpiSize; rank++) 
     627                nbNeighbourNodes += nbRecvNode[rank]; 
     628 
     629        neighbourElements = new Elt[nbNeighbourNodes]; 
     630        nbNeighbourElements = nbNeighbourNodes; 
     631 
     632        int index = 0; 
     633        for (int rank = 0; rank < mpiSize; rank++) 
     634        { 
     635                pos[rank] = 0; 
     636                for (int j = 0; j < nbRecvNode[rank]; j++) 
     637                { 
     638                        unpackPolygon(neighbourElements[index], recvBuffer2[rank], pos[rank]); 
     639                        neighbourElements[index].id.ind = sstree.localTree.leafs.size() + index; 
     640                        index++; 
     641                } 
     642        } 
     643        for (int rank = 0; rank < mpiSize; rank++) 
     644        { 
     645                if (nbRecvNode[rank] > 0) delete [] recvBuffer2[rank]; 
     646                if (nbSendNode[rank] > 0) delete [] sendBuffer2[rank]; 
     647        } 
     648        delete [] recvBuffer2; 
     649        delete [] sendBuffer2; 
     650        delete [] sendMessageSize; 
     651        delete [] recvMessageSize; 
     652        delete [] nbSendNode; 
     653        delete [] nbRecvNode; 
     654        delete [] sendRequest; 
     655        delete [] recvRequest; 
     656        delete [] status; 
     657        delete [] pos; 
     658 
     659        /* re-compute on received elements to avoid having to send this information */ 
     660        neighbourNodes.resize(nbNeighbourNodes); 
     661        setCirclesAndLinks(neighbourElements, neighbourNodes); 
     662        cptAllEltsGeom(neighbourElements, nbNeighbourNodes, srcGrid.pole); 
     663 
     664        /* the local SS tree must include nodes from other cpus if they are potential 
     665           intersector of a local node */ 
     666        sstree.localTree.insertNodes(neighbourNodes); 
     667 
     668        /* for every local element, 
     669           use the SS-tree to find all elements (including neighbourElements) 
     670           who are potential neighbours because their circles intersect, 
     671           then check all canditates for common edges to build up connectivity information 
     672        */ 
     673        for (int j = 0; j < sstree.localTree.leafs.size(); j++) 
     674        { 
     675                Node& node = sstree.localTree.leafs[j]; 
     676 
     677                /* find all leafs whoes circles that intersect node's circle and save into node->intersectors */ 
     678                node.search(sstree.localTree.root); 
     679 
     680                Elt *elt = (Elt *)(node.data); 
     681 
     682                for (int i = 0; i < elt->n; i++) elt->neighbour[i] = NOT_FOUND; 
     683 
     684                /* for element `elt` loop through all nodes in the SS-tree 
     685                   whoes circles intersect with the circle around `elt` (the SS intersectors) 
     686                   and check if they are neighbours in the sense that the two elements share an edge. 
     687                   If they do, save this information for elt */ 
     688                for (list<NodePtr>::iterator it = (node.intersectors).begin(); it != (node.intersectors).end(); ++it) 
     689                { 
     690                        Elt *elt2 = (Elt *)((*it)->data); 
     691                        set_neighbour(*elt, *elt2); 
     692                } 
     693 
     694        } 
    702695} 
    703696 
     
    705698void Mapper::computeIntersection(Elt *elements, int nbElements) 
    706699{ 
    707     int mpiSize, mpiRank; 
    708     MPI_Comm_size(communicator, &mpiSize); 
    709     MPI_Comm_rank(communicator, &mpiRank); 
    710  
    711     MPI_Barrier(communicator); 
    712  
    713     vector<Node> *routingList = new vector<Node>[mpiSize]; 
    714  
    715     vector<Node> routeNodes;  routeNodes.reserve(nbElements); 
    716     for (int j = 0; j < nbElements; j++) 
    717     { 
    718         elements[j].id.ind = j; 
    719         elements[j].id.rank = mpiRank; 
    720         routeNodes.push_back(Node(elements[j].x, cptRadius(elements[j]), &elements[j])); 
    721     } 
    722  
    723     vector<vector<int> > routes(routeNodes.size()); 
    724     sstree.routeIntersections(routes, routeNodes); 
    725     for (int i = 0; i < routes.size(); ++i) 
    726         for (int k = 0; k < routes[i].size(); ++k) 
    727             routingList[routes[i][k]].push_back(routeNodes[i]); 
    728  
    729     if (verbose >= 2) 
    730     { 
    731         cout << " --> rank  " << mpiRank << " nbElements " << nbElements << " : "; 
    732         for (int rank = 0; rank < mpiSize; rank++) 
    733             cout << routingList[rank].size() << "   "; 
    734         cout << endl; 
    735     } 
    736     MPI_Barrier(communicator); 
    737  
    738     int *nbSendNode = new int[mpiSize]; 
    739     int *nbRecvNode = new int[mpiSize]; 
    740     int *sentMessageSize = new int[mpiSize]; 
    741     int *recvMessageSize = new int[mpiSize]; 
    742  
    743     for (int rank = 0; rank < mpiSize; rank++) 
    744     { 
    745         nbSendNode[rank] = routingList[rank].size(); 
    746         sentMessageSize[rank] = 0; 
    747         for (size_t j = 0; j < routingList[rank].size(); j++) 
    748         { 
    749             Elt *elt = (Elt *) (routingList[rank][j].data); 
    750             sentMessageSize[rank] += packedPolygonSize(*elt); 
    751         } 
    752     } 
    753  
    754     MPI_Alltoall(nbSendNode, 1, MPI_INT, nbRecvNode, 1, MPI_INT, communicator); 
    755     MPI_Alltoall(sentMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator); 
    756  
    757     int total = 0; 
    758  
    759     for (int rank = 0; rank < mpiSize; rank++) 
    760     { 
    761         total = total + nbRecvNode[rank]; 
    762     } 
    763  
    764     if (verbose >= 2) cout << "---> rank " << mpiRank << " : compute intersection : total received nodes  " << total << endl; 
    765     char **sendBuffer = new char*[mpiSize]; 
    766     char **recvBuffer = new char*[mpiSize]; 
    767     int *pos = new int[mpiSize]; 
    768  
    769     for (int rank = 0; rank < mpiSize; rank++) 
    770     { 
    771         if (nbSendNode[rank] > 0) sendBuffer[rank] = new char[sentMessageSize[rank]]; 
    772         if (nbRecvNode[rank] > 0) recvBuffer[rank] = new char[recvMessageSize[rank]]; 
    773     } 
    774  
    775     for (int rank = 0; rank < mpiSize; rank++) 
    776     { 
    777         pos[rank] = 0; 
    778         for (size_t j = 0; j < routingList[rank].size(); j++) 
    779         { 
    780             Elt* elt = (Elt *) (routingList[rank][j].data); 
    781             packPolygon(*elt, sendBuffer[rank], pos[rank]); 
    782         } 
    783     } 
    784     delete [] routingList; 
    785  
    786     int nbSendRequest = 0; 
    787     int nbRecvRequest = 0; 
    788     MPI_Request *sendRequest = new MPI_Request[mpiSize]; 
    789     MPI_Request *recvRequest = new MPI_Request[mpiSize]; 
    790     MPI_Status   *status = new MPI_Status[mpiSize]; 
    791  
    792     for (int rank = 0; rank < mpiSize; rank++) 
    793     { 
    794         if (nbSendNode[rank] > 0) 
    795         { 
    796             MPI_Issend(sendBuffer[rank], sentMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]); 
    797             nbSendRequest++; 
    798         } 
    799         if (nbRecvNode[rank] > 0) 
    800         { 
    801             MPI_Irecv(recvBuffer[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
    802             nbRecvRequest++; 
    803         } 
    804     } 
    805  
    806     MPI_Waitall(nbRecvRequest, recvRequest, status); 
    807     MPI_Waitall(nbSendRequest, sendRequest, status); 
    808  
    809     char **sendBuffer2 = new char*[mpiSize]; 
    810     char **recvBuffer2 = new char*[mpiSize]; 
    811  
    812     double tic = cputime(); 
    813     for (int rank = 0; rank < mpiSize; rank++) 
    814     { 
    815         sentMessageSize[rank] = 0; 
    816  
    817         if (nbRecvNode[rank] > 0) 
    818         { 
    819             Elt *recvElt = new Elt[nbRecvNode[rank]]; 
    820             pos[rank] = 0; 
    821             for (int j = 0; j < nbRecvNode[rank]; j++) 
    822             { 
    823                 unpackPolygon(recvElt[j], recvBuffer[rank], pos[rank]); 
    824                 cptEltGeom(recvElt[j], tgtGrid.pole); 
    825                 Node recvNode(recvElt[j].x, cptRadius(recvElt[j]), &recvElt[j]); 
    826                 recvNode.search(sstree.localTree.root); 
    827                 /* for a node holding an element of the target, loop throught candidates for intersecting source */ 
    828                 for (list<NodePtr>::iterator it = (recvNode.intersectors).begin(); it != (recvNode.intersectors).end(); ++it) 
    829                 { 
    830                     Elt *elt2 = (Elt *) ((*it)->data); 
    831                     /* recvElt is target, elt2 is source */ 
    832                     intersect(&recvElt[j], elt2); 
    833                     //intersect_ym(&recvElt[j], elt2); 
    834                 } 
    835                 if (recvElt[j].is.size() > 0) sentMessageSize[rank] += packIntersectionSize(recvElt[j]); 
    836  
    837                 // here recvNode goes out of scope 
    838             } 
    839  
    840             if (sentMessageSize[rank] > 0) 
    841             { 
    842                 sentMessageSize[rank] += sizeof(int); 
    843                 sendBuffer2[rank] = new char[sentMessageSize[rank]]; 
    844                 *((int *) sendBuffer2[rank]) = 0; 
    845                 pos[rank] = sizeof(int); 
    846                 for (int j = 0; j < nbRecvNode[rank]; j++) 
    847                 { 
    848                     packIntersection(recvElt[j], sendBuffer2[rank], pos[rank]); 
    849                     //FIXME should be deleted: recvElt[j].delete_intersections(); // intersection areas have been packed to buffer and won't be used any more 
    850                 } 
    851             } 
    852             delete [] recvElt; 
    853         } 
    854     } 
    855     delete [] pos; 
    856  
    857     for (int rank = 0; rank < mpiSize; rank++) 
    858     { 
    859         if (nbSendNode[rank] > 0) delete [] sendBuffer[rank]; 
    860         if (nbRecvNode[rank] > 0) delete [] recvBuffer[rank]; 
    861         nbSendNode[rank] = 0; 
    862     } 
    863  
    864     if (verbose >= 2) cout << "Rank " << mpiRank << "  Compute (internal) intersection " << cputime() - tic << " s" << endl; 
    865     MPI_Alltoall(sentMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator); 
    866  
    867     for (int rank = 0; rank < mpiSize; rank++) 
    868         if (recvMessageSize[rank] > 0) 
    869             recvBuffer2[rank] = new char[recvMessageSize[rank]]; 
    870  
    871     nbSendRequest = 0; 
    872     nbRecvRequest = 0; 
    873  
    874     for (int rank = 0; rank < mpiSize; rank++) 
    875     { 
    876         if (sentMessageSize[rank] > 0) 
    877         { 
    878             MPI_Issend(sendBuffer2[rank], sentMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]); 
    879             nbSendRequest++; 
    880         } 
    881         if (recvMessageSize[rank] > 0) 
    882         { 
    883             MPI_Irecv(recvBuffer2[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
    884             nbRecvRequest++; 
    885         } 
    886     } 
    887      
    888     MPI_Waitall(nbRecvRequest, recvRequest, status); 
    889     MPI_Waitall(nbSendRequest, sendRequest, status); 
    890  
    891     delete [] sendRequest; 
    892     delete [] recvRequest; 
    893     delete [] status; 
    894     for (int rank = 0; rank < mpiSize; rank++) 
    895     { 
    896         if (nbRecvNode[rank] > 0) 
    897         { 
    898             if (sentMessageSize[rank] > 0) 
    899                 delete [] sendBuffer2[rank]; 
    900         } 
    901  
    902         if (recvMessageSize[rank] > 0) 
    903         { 
    904             unpackIntersection(elements, recvBuffer2[rank]); 
    905             delete [] recvBuffer2[rank]; 
    906         } 
    907     } 
    908     delete [] sendBuffer2; 
    909     delete [] recvBuffer2; 
    910     delete [] sendBuffer; 
    911     delete [] recvBuffer; 
    912  
    913     delete [] nbSendNode; 
    914     delete [] nbRecvNode; 
    915     delete [] sentMessageSize; 
    916     delete [] recvMessageSize; 
     700        int mpiSize, mpiRank; 
     701        MPI_Comm_size(communicator, &mpiSize); 
     702        MPI_Comm_rank(communicator, &mpiRank); 
     703 
     704        MPI_Barrier(communicator); 
     705 
     706        vector<Node> *routingList = new vector<Node>[mpiSize]; 
     707 
     708        vector<Node> routeNodes;  routeNodes.reserve(nbElements); 
     709        for (int j = 0; j < nbElements; j++) 
     710        { 
     711                elements[j].id.ind = j; 
     712                elements[j].id.rank = mpiRank; 
     713                routeNodes.push_back(Node(elements[j].x, cptRadius(elements[j]), &elements[j])); 
     714        } 
     715 
     716        vector<vector<int> > routes(routeNodes.size()); 
     717        sstree.routeIntersections(routes, routeNodes); 
     718        for (int i = 0; i < routes.size(); ++i) 
     719                for (int k = 0; k < routes[i].size(); ++k) 
     720                        routingList[routes[i][k]].push_back(routeNodes[i]); 
     721 
     722        if (verbose >= 2) 
     723        { 
     724                cout << " --> rank  " << mpiRank << " nbElements " << nbElements << " : "; 
     725                for (int rank = 0; rank < mpiSize; rank++) 
     726                        cout << routingList[rank].size() << "   "; 
     727                cout << endl; 
     728        } 
     729        MPI_Barrier(communicator); 
     730 
     731        int *nbSendNode = new int[mpiSize]; 
     732        int *nbRecvNode = new int[mpiSize]; 
     733        int *sentMessageSize = new int[mpiSize]; 
     734        int *recvMessageSize = new int[mpiSize]; 
     735 
     736        for (int rank = 0; rank < mpiSize; rank++) 
     737        { 
     738                nbSendNode[rank] = routingList[rank].size(); 
     739                sentMessageSize[rank] = 0; 
     740                for (size_t j = 0; j < routingList[rank].size(); j++) 
     741                { 
     742                        Elt *elt = (Elt *) (routingList[rank][j].data); 
     743                        sentMessageSize[rank] += packedPolygonSize(*elt); 
     744                } 
     745        } 
     746 
     747        MPI_Alltoall(nbSendNode, 1, MPI_INT, nbRecvNode, 1, MPI_INT, communicator); 
     748        MPI_Alltoall(sentMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator); 
     749 
     750        int total = 0; 
     751 
     752        for (int rank = 0; rank < mpiSize; rank++) 
     753        { 
     754                total = total + nbRecvNode[rank]; 
     755        } 
     756 
     757        if (verbose >= 2) cout << "---> rank " << mpiRank << " : compute intersection : total received nodes  " << total << endl; 
     758        char **sendBuffer = new char*[mpiSize]; 
     759        char **recvBuffer = new char*[mpiSize]; 
     760        int *pos = new int[mpiSize]; 
     761 
     762        for (int rank = 0; rank < mpiSize; rank++) 
     763        { 
     764                if (nbSendNode[rank] > 0) sendBuffer[rank] = new char[sentMessageSize[rank]]; 
     765                if (nbRecvNode[rank] > 0) recvBuffer[rank] = new char[recvMessageSize[rank]]; 
     766        } 
     767 
     768        for (int rank = 0; rank < mpiSize; rank++) 
     769        { 
     770                pos[rank] = 0; 
     771                for (size_t j = 0; j < routingList[rank].size(); j++) 
     772                { 
     773                        Elt* elt = (Elt *) (routingList[rank][j].data); 
     774                        packPolygon(*elt, sendBuffer[rank], pos[rank]); 
     775                } 
     776        } 
     777        delete [] routingList; 
     778 
     779        int nbSendRequest = 0; 
     780        int nbRecvRequest = 0; 
     781        MPI_Request *sendRequest = new MPI_Request[mpiSize]; 
     782        MPI_Request *recvRequest = new MPI_Request[mpiSize]; 
     783        MPI_Status   *status = new MPI_Status[mpiSize]; 
     784 
     785        for (int rank = 0; rank < mpiSize; rank++) 
     786        { 
     787                if (nbSendNode[rank] > 0) 
     788                { 
     789                        MPI_Issend(sendBuffer[rank], sentMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]); 
     790                        nbSendRequest++; 
     791                } 
     792                if (nbRecvNode[rank] > 0) 
     793                { 
     794                        MPI_Irecv(recvBuffer[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
     795                        nbRecvRequest++; 
     796                } 
     797        } 
     798 
     799        MPI_Waitall(nbRecvRequest, recvRequest, status); 
     800        MPI_Waitall(nbSendRequest, sendRequest, status); 
     801 
     802        char **sendBuffer2 = new char*[mpiSize]; 
     803        char **recvBuffer2 = new char*[mpiSize]; 
     804 
     805        double tic = cputime(); 
     806        for (int rank = 0; rank < mpiSize; rank++) 
     807        { 
     808                sentMessageSize[rank] = 0; 
     809 
     810                if (nbRecvNode[rank] > 0) 
     811                { 
     812                        Elt *recvElt = new Elt[nbRecvNode[rank]]; 
     813                        pos[rank] = 0; 
     814                        for (int j = 0; j < nbRecvNode[rank]; j++) 
     815                        { 
     816                                unpackPolygon(recvElt[j], recvBuffer[rank], pos[rank]); 
     817                                cptEltGeom(recvElt[j], tgtGrid.pole); 
     818                                Node recvNode(recvElt[j].x, cptRadius(recvElt[j]), &recvElt[j]); 
     819                                recvNode.search(sstree.localTree.root); 
     820                                /* for a node holding an element of the target, loop throught candidates for intersecting source */ 
     821                                for (list<NodePtr>::iterator it = (recvNode.intersectors).begin(); it != (recvNode.intersectors).end(); ++it) 
     822                                { 
     823                                        Elt *elt2 = (Elt *) ((*it)->data); 
     824                                        /* recvElt is target, elt2 is source */ 
     825//                                      intersect(&recvElt[j], elt2); 
     826                                        intersect_ym(&recvElt[j], elt2); 
     827                                } 
     828                                if (recvElt[j].is.size() > 0) sentMessageSize[rank] += packIntersectionSize(recvElt[j]); 
     829 
     830                                // here recvNode goes out of scope 
     831                        } 
     832 
     833                        if (sentMessageSize[rank] > 0) 
     834                        { 
     835                                sentMessageSize[rank] += sizeof(int); 
     836                                sendBuffer2[rank] = new char[sentMessageSize[rank]]; 
     837                                *((int *) sendBuffer2[rank]) = 0; 
     838                                pos[rank] = sizeof(int); 
     839                                for (int j = 0; j < nbRecvNode[rank]; j++) 
     840                                { 
     841                                        packIntersection(recvElt[j], sendBuffer2[rank], pos[rank]); 
     842                                        //FIXME should be deleted: recvElt[j].delete_intersections(); // intersection areas have been packed to buffer and won't be used any more 
     843                                } 
     844                        } 
     845                        delete [] recvElt; 
     846                } 
     847        } 
     848        delete [] pos; 
     849 
     850        for (int rank = 0; rank < mpiSize; rank++) 
     851        { 
     852                if (nbSendNode[rank] > 0) delete [] sendBuffer[rank]; 
     853                if (nbRecvNode[rank] > 0) delete [] recvBuffer[rank]; 
     854                nbSendNode[rank] = 0; 
     855        } 
     856 
     857        if (verbose >= 2) cout << "Rank " << mpiRank << "  Compute (internal) intersection " << cputime() - tic << " s" << endl; 
     858        MPI_Alltoall(sentMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator); 
     859 
     860        for (int rank = 0; rank < mpiSize; rank++) 
     861                if (recvMessageSize[rank] > 0) 
     862                        recvBuffer2[rank] = new char[recvMessageSize[rank]]; 
     863 
     864        nbSendRequest = 0; 
     865        nbRecvRequest = 0; 
     866 
     867        for (int rank = 0; rank < mpiSize; rank++) 
     868        { 
     869                if (sentMessageSize[rank] > 0) 
     870                { 
     871                        MPI_Issend(sendBuffer2[rank], sentMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]); 
     872                        nbSendRequest++; 
     873                } 
     874                if (recvMessageSize[rank] > 0) 
     875                { 
     876                        MPI_Irecv(recvBuffer2[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
     877                        nbRecvRequest++; 
     878                } 
     879        } 
     880 
     881  MPI_Waitall(nbRecvRequest, recvRequest, status); 
     882  MPI_Waitall(nbSendRequest, sendRequest, status); 
     883 
     884        delete [] sendRequest; 
     885        delete [] recvRequest; 
     886        delete [] status; 
     887        for (int rank = 0; rank < mpiSize; rank++) 
     888        { 
     889                if (nbRecvNode[rank] > 0) 
     890                { 
     891                        if (sentMessageSize[rank] > 0) 
     892                                delete [] sendBuffer2[rank]; 
     893                } 
     894 
     895                if (recvMessageSize[rank] > 0) 
     896                { 
     897                        unpackIntersection(elements, recvBuffer2[rank]); 
     898                        delete [] recvBuffer2[rank]; 
     899                } 
     900        } 
     901        delete [] sendBuffer2; 
     902        delete [] recvBuffer2; 
     903        delete [] sendBuffer; 
     904        delete [] recvBuffer; 
     905 
     906        delete [] nbSendNode; 
     907        delete [] nbRecvNode; 
     908        delete [] sentMessageSize; 
     909        delete [] recvMessageSize; 
    917910} 
    918911 
    919912Mapper::~Mapper() 
    920913{ 
    921     delete [] remapMatrix; 
    922     delete [] srcAddress; 
    923     delete [] srcRank; 
    924     delete [] dstAddress; 
    925     if (neighbourElements) delete [] neighbourElements; 
    926 } 
    927  
    928 } 
     914        delete [] remapMatrix; 
     915        delete [] srcAddress; 
     916        delete [] srcRank; 
     917        delete [] dstAddress; 
     918        if (neighbourElements) delete [] neighbourElements; 
     919} 
     920 
     921} 
  • XIOS/dev/branch_openmp/extern/remap/src/mapper.hpp

    r1134 r1328  
    33#include "parallel_tree.hpp" 
    44#include "mpi.hpp" 
    5  
    6 #ifdef _usingEP 
    7 #include "ep_declaration.hpp" 
    8 #endif 
    95 
    106namespace sphereRemap { 
     
    2218{ 
    2319public: 
    24        Mapper(ep_lib::MPI_Comm comm=MPI_COMM_WORLD) : communicator(comm), verbose(SILENT), neighbourElements(NULL), sstree(comm) {} 
     20       Mapper(ep_lib::MPI_Comm comm) : communicator(comm), verbose(SILENT), neighbourElements(NULL), sstree(comm) {} 
    2521       ~Mapper(); 
    2622       void setVerbosity(verbosity v) {verbose=v ;} 
  • XIOS/dev/branch_openmp/extern/remap/src/mpi_cascade.cpp

    r688 r1328  
    11#include "mpi_cascade.hpp" 
    22#include <iostream> 
     3using namespace ep_lib; 
    34 
    45namespace sphereRemap { 
  • XIOS/dev/branch_openmp/extern/remap/src/mpi_cascade.hpp

    r694 r1328  
    1212{ 
    1313public: 
    14         CCascadeLevel(MPI_Comm comm) : comm(comm) 
    15         { 
    16                 MPI_Comm_size(comm, &size); 
    17                 MPI_Comm_rank(comm, &rank); 
    18         } 
    19         int colour() const { return rank % group_size; }; 
    20         int key() const { return p_colour() + rank/(p_grp_size*group_size)*p_grp_size; } 
     14  CCascadeLevel(ep_lib::MPI_Comm comm) : comm(comm) 
     15  { 
     16    ep_lib::MPI_Comm_size(comm, &size); 
     17    ep_lib::MPI_Comm_rank(comm, &rank); 
     18  } 
     19  int colour() const { return rank % group_size; }; 
     20  int key() const { return p_colour() + rank/(p_grp_size*group_size)*p_grp_size; } 
    2121 
    22         // perpendicular group 
    23         int p_colour() const { return (rank%group_size + rank/group_size) % p_grp_size; } 
    24         int p_key() const { return colour() + rank/(p_grp_size*group_size)*group_size; } 
     22  // perpendicular group 
     23  int p_colour() const { return (rank%group_size + rank/group_size) % p_grp_size; } 
     24  int p_key() const { return colour() + rank/(p_grp_size*group_size)*group_size; } 
    2525 
    26         MPI_Comm comm, pg_comm; 
    27         int rank; 
    28         int size; 
    29         int group_size; // group_size and p_grp_size are interchanged?? FIXME 
    30         int p_grp_size; 
     26  ep_lib::MPI_Comm comm, pg_comm; 
     27  int rank; 
     28  int size; 
     29  int group_size; // group_size and p_grp_size are interchanged?? FIXME 
     30  int p_grp_size; 
    3131}; 
    3232 
     
    3434{ 
    3535public: 
    36         //  
    37         CMPICascade(int nodes_per_level, MPI_Comm comm); 
     36  CMPICascade(int nodes_per_level, ep_lib::MPI_Comm comm); 
    3837 
    39         int num_levels; 
    40         std::vector<CCascadeLevel> level; 
     38  int num_levels; 
     39  std::vector<CCascadeLevel> level; 
    4140}; 
    4241 
  • XIOS/dev/branch_openmp/extern/remap/src/mpi_routing.cpp

    r1289 r1328  
    55#include "timerRemap.hpp" 
    66#include <iostream> 
    7 #ifdef _usingEP 
    8 #include "ep_declaration.hpp" 
    9 #endif 
     7using namespace ep_lib; 
    108 
    119namespace sphereRemap { 
     
    125123        CTimer::get("CMPIRouting::init(reduce_scatter)").print(); 
    126124 
    127         MPI_Info info_null; 
    128  
    129         // MPI_Alloc_mem(nbTarget *sizeof(int), info_null, &targetRank); 
    130         // MPI_Alloc_mem(nbSource *sizeof(int), info_null, &sourceRank); 
    131125        MPI_Alloc_mem(nbTarget *sizeof(int), MPI_INFO_NULL, &targetRank); 
    132126        MPI_Alloc_mem(nbSource *sizeof(int), MPI_INFO_NULL, &sourceRank); 
     
    158152        { 
    159153                #ifdef _usingEP 
    160                 MPI_Irecv(&sourceRank[i], 1, MPI_INT, -1, 0, communicator, &request[indexRequest]); 
     154                MPI_Irecv(&sourceRank[i], 1, MPI_INT, -2, 0, communicator, &request[indexRequest]); 
    161155                #else 
    162156                MPI_Irecv(&sourceRank[i], 1, MPI_INT, MPI_ANY_SOURCE, 0, communicator, &request[indexRequest]); 
     
    182176        { 
    183177                #ifdef _usingEP 
    184                 MPI_Irecv(&sourceRank[i], 1, MPI_INT, -1, 0, communicator, &request[indexRequest]); 
     178                MPI_Irecv(&sourceRank[i], 1, MPI_INT, -2, 0, communicator, &request[indexRequest]); 
    185179                #else 
    186180                MPI_Irecv(&sourceRank[i], 1, MPI_INT, MPI_ANY_SOURCE, 0, communicator, &request[indexRequest]); 
  • XIOS/dev/branch_openmp/extern/remap/src/mpi_routing.hpp

    r694 r1328  
    1111{ 
    1212 
    13         MPI_Comm communicator; 
     13        ep_lib::MPI_Comm communicator; 
    1414        int mpiRank; 
    1515        int mpiSize; 
     
    2929 
    3030public: 
    31         CMPIRouting(MPI_Comm comm); 
     31        CMPIRouting(ep_lib::MPI_Comm comm); 
    3232        ~CMPIRouting(); 
    3333        template<typename T> void init(const std::vector<T>& route, CMPICascade *cascade = NULL); 
     
    4444template <typename T> 
    4545void alltoalls_known(const std::vector<std::vector<T> >& send, std::vector<std::vector<T> >& recv, 
    46                      const std::vector<int>& ranks, MPI_Comm communicator); 
     46                     const std::vector<int>& ranks, ep_lib::MPI_Comm communicator); 
    4747 
    4848template <typename T> 
    4949void alltoalls_unknown(const std::vector<std::vector<T> >& send, std::vector<std::vector<T> >& recv, 
    50                        const std::vector<int>& ranks, MPI_Comm communicator); 
     50                       const std::vector<int>& ranks, ep_lib::MPI_Comm communicator); 
    5151} 
    5252#endif 
  • XIOS/dev/branch_openmp/extern/remap/src/node.cpp

    r1205 r1328  
    254254NodePtr insert(NodePtr thIs, NodePtr node) 
    255255{ 
    256   int la = thIs->level; // node to be inserted 
    257   int lb = node->level; // node where insertation 
    258   assert(la < lb); // node to be inserted must have lower level then parent 
    259   //if (thIs->parent) assert(find_in_tree1(thIs) == true); 
    260   NodePtr q = NULL; 
    261   NodePtr chd = NULL; 
    262   node->move(thIs); 
    263   if (la == lb - 1) 
    264   { 
     256        int la = thIs->level; // node to be inserted 
     257        int lb = node->level; // node where insertation 
     258        assert(la < lb); // node to be inserted must have lower level then parent 
     259        //if (thIs->parent) assert(find_in_tree1(thIs) == true); 
     260        NodePtr q = NULL; 
     261        NodePtr chd = NULL; 
     262        node->move(thIs); 
     263        if (la == lb - 1) 
     264        { 
    265265    node->child.push_back(thIs); 
    266     thIs->parent = node; 
    267     if (node->child.size() > MAX_NODE_SZ &&  node->tree->canSplit() ) // with us as additional child `node` is now too large :( 
    268     return (node->reinserted || node->parent == NULL) ? split(node) : reinsert(node); 
    269   } 
    270   else // la < lb - 1 
    271   { 
    272     chd = thIs->closest(node->child); 
    273     q = insert(thIs, chd); 
    274   } 
    275   if ((node->updateCount + 1) % UPDATE_EVERY == 0) 
    276     node->update(); 
    277   else 
    278   { 
    279     if (q) node->remove(q); 
    280     node->inflate(chd); 
    281   } 
     266                thIs->parent = node; 
     267                if (node->child.size() > MAX_NODE_SZ &&  node->tree->canSplit() ) // with us as additional child `node` is now too large :( 
     268                        return (node->reinserted || node->parent == NULL) ? split(node) : reinsert(node); 
     269        } 
     270        else // la < lb - 1 
     271        { 
     272                chd = thIs->closest(node->child); 
     273                q = insert(thIs, chd); 
     274        } 
     275        if ((node->updateCount + 1) % UPDATE_EVERY == 0) 
     276                node->update(); 
     277        else 
     278        { 
     279                if (q) node->remove(q); 
     280                node->inflate(chd); 
     281        } 
    282282 
    283283  return q; 
  • XIOS/dev/branch_openmp/extern/remap/src/node.hpp

    r1153 r1328  
    1515struct Circle 
    1616{ 
    17   Coord centre; 
    18   double radius; 
     17        Coord centre; 
     18        double radius; 
    1919}; 
    2020 
     
    116116struct Node 
    117117{ 
    118   int level; /* FIXME leafs are 0 and root is max level? */ 
    119   int leafCount; /* number of leafs that are descendants of this node (the elements in it's cycle) */ 
    120   Coord centre; 
    121   double radius; 
    122   NodePtr parent, ref; 
    123   std::vector<NodePtr> child; 
    124   std::list<NodePtr> intersectors; 
    125   bool reinserted; 
    126   int updateCount;  // double var; 
    127   CBasicTree* tree; 
    128   void *data; 
    129   int route; 
     118        int level; /* FIXME leafs are 0 and root is max level? */ 
     119        int leafCount; /* number of leafs that are descendants of this node (the elements in it's cycle) */ 
     120        Coord centre; 
     121        double radius; 
     122        NodePtr parent, ref; 
     123        std::vector<NodePtr> child; 
     124        std::list<NodePtr> intersectors; 
     125        bool reinserted; 
     126        int updateCount;  // double var; 
     127        CBasicTree* tree; 
     128        void *data; 
     129        int route; 
    130130  bool toDelete ; 
    131131 
    132   Node() : level(0), leafCount(1), centre(ORIGIN), radius(0), reinserted(false), updateCount(0), toDelete(false) {} 
    133   Node(const Coord& centre, double radius, void *data) 
    134     : level(0), leafCount(1), centre(centre), radius(radius), reinserted(false), updateCount(0), data(data), toDelete(false) {} 
     132        Node() : level(0), leafCount(1), centre(ORIGIN), radius(0), reinserted(false), updateCount(0), toDelete(false) {} 
     133        Node(const Coord& centre, double radius, void *data) 
     134                : level(0), leafCount(1), centre(centre), radius(radius), reinserted(false), updateCount(0), data(data), toDelete(false) {} 
    135135 
    136136//#ifdef DEBUG 
     
    178178//#endif 
    179179 
    180   void move(const NodePtr node); 
    181   void remove(const NodePtr node); 
    182   void inflate(const NodePtr node); 
    183   void update(); 
     180        void move(const NodePtr node); 
     181        void remove(const NodePtr node); 
     182        void inflate(const NodePtr node); 
     183        void update(); 
    184184  void output(std::ostream& flux, int level, int color) ; 
    185   NodePtr closest(std::vector<NodePtr>& list, int n = CLOSEST); 
    186   NodePtr farthest(std::vector<NodePtr>& list); 
    187   void findClosest(int level, NodePtr src, double& minDist, NodePtr &closest); 
    188  
    189   void search(NodePtr node); 
    190   bool centreInside(Node &node); 
    191   bool intersects(NodePtr node); 
    192   bool isInside(Node &node); 
    193   int incluCheck(); 
     185        NodePtr closest(std::vector<NodePtr>& list, int n = CLOSEST); 
     186        NodePtr farthest(std::vector<NodePtr>& list); 
     187        void findClosest(int level, NodePtr src, double& minDist, NodePtr &closest); 
     188 
     189        void search(NodePtr node); 
     190        bool centreInside(Node &node); 
     191        bool intersects(NodePtr node); 
     192        bool isInside(Node &node); 
     193        int incluCheck(); 
    194194  void checkParent(void) ; 
    195   void printChildren(); 
    196   void assignRoute(std::vector<int>::iterator& rank, int level); 
    197   void assignCircleAndPropagateUp(Coord *centres, double *radia, int level); 
    198   void printLevel(int level); 
    199   void routeNode(NodePtr node, int level); 
    200   void routingIntersecting(std::vector<Node>* routingList, NodePtr node); 
    201   void routeIntersection(std::vector<int>& routes, NodePtr node); 
     195        void printChildren(); 
     196        void assignRoute(std::vector<int>::iterator& rank, int level); 
     197        void assignCircleAndPropagateUp(Coord *centres, double *radia, int level); 
     198        void printLevel(int level); 
     199        void routeNode(NodePtr node, int level); 
     200        void routingIntersecting(std::vector<Node>* routingList, NodePtr node); 
     201        void routeIntersection(std::vector<int>& routes, NodePtr node); 
    202202  void getNodeLevel(int level,std::list<NodePtr>& NodeList) ; 
    203203  bool removeDeletedNodes(int assignLevel) ; 
  • XIOS/dev/branch_openmp/extern/remap/src/parallel_tree.cpp

    r1295 r1328  
    1212 
    1313#include "parallel_tree.hpp" 
     14using namespace ep_lib; 
    1415 
    1516namespace sphereRemap { 
    16  
    17 extern CRemapGrid srcGrid; 
    18 #pragma omp threadprivate(srcGrid) 
    19  
    20 extern CRemapGrid tgtGrid; 
    21 #pragma omp threadprivate(tgtGrid) 
    2217 
    2318static const int assignLevel = 2; 
     
    292287{ 
    293288 
    294   int assignLevel = 2; 
    295   int nbSampleNodes = 2*ipow(MAX_NODE_SZ + 1, assignLevel); 
     289        int assignLevel = 2; 
     290        int nbSampleNodes = 2*ipow(MAX_NODE_SZ + 1, assignLevel); 
    296291 
    297292 
     
    305300  MPI_Comm_size(communicator,&commSize) ; 
    306301   
    307   // make multiple of two 
    308   nbSampleNodes /= 2; 
    309   nbSampleNodes *= 2; 
    310   //assert( nbTot > nbSampleNodes*commSize) ; 
     302        // make multiple of two 
     303        nbSampleNodes /= 2; 
     304        nbSampleNodes *= 2; 
     305//  assert( nbTot > nbSampleNodes*commSize) ; 
    311306     
    312307  int nbSampleNodes1 = nbSampleNodes * (nb1*commSize)/(1.*nbTot) ; 
     
    314309   
    315310 
    316   //assert(node.size() > nbSampleNodes); 
    317   //assert(node2.size() > nbSampleNodes); 
    318   //assert(node.size() + node2.size() > nbSampleNodes); 
    319   vector<Node> sampleNodes; sampleNodes.reserve(nbSampleNodes1+nbSampleNodes2); 
    320  
    321   vector<int> randomArray1(node.size()); 
    322   randomizeArray(randomArray1); 
    323   vector<int> randomArray2(node2.size()); 
    324   randomizeArray(randomArray2); 
    325  
    326   for (int i = 0; i <nbSampleNodes1; i++) sampleNodes.push_back(Node(node[randomArray1[i%nb1]].centre,  node[randomArray1[i%nb1]].radius, NULL)); 
    327   for (int i = 0; i <nbSampleNodes2; i++) sampleNodes.push_back(Node(node2[randomArray2[i%nb2]].centre, node2[randomArray2[i%nb2]].radius, NULL)); 
     311//      assert(node.size() > nbSampleNodes); 
     312//      assert(node2.size() > nbSampleNodes); 
     313//      assert(node.size() + node2.size() > nbSampleNodes); 
     314        vector<Node> sampleNodes; sampleNodes.reserve(nbSampleNodes1+nbSampleNodes2); 
     315 
     316        vector<int> randomArray1(node.size()); 
     317        randomizeArray(randomArray1); 
     318        vector<int> randomArray2(node2.size()); 
     319        randomizeArray(randomArray2); 
     320 
     321        for (int i = 0; i <nbSampleNodes1; i++) sampleNodes.push_back(Node(node[randomArray1[i%nb1]].centre,  node[randomArray1[i%nb1]].radius, NULL)); 
     322        for (int i = 0; i <nbSampleNodes2; i++) sampleNodes.push_back(Node(node2[randomArray2[i%nb2]].centre, node2[randomArray2[i%nb2]].radius, NULL)); 
    328323 
    329324        CTimer::get("buildParallelSampleTree").resume(); 
     
    336331        CTimer::get("parallelRouteNode").resume(); 
    337332        vector<int> route(node.size()); 
     333        cout<<"node.size = "<<node.size()<<endl; 
    338334        routeNodes(route /*out*/, node); 
    339335        CTimer::get("parallelRouteNode").suspend(); 
  • XIOS/dev/branch_openmp/extern/remap/src/parallel_tree.hpp

    r1134 r1328  
    66#include "mpi_cascade.hpp" 
    77#include "mpi.hpp" 
    8 #ifdef _usingEP 
    9 #include "ep_declaration.hpp" 
    10 #endif 
    11  
    128namespace sphereRemap { 
    139 
     
    1511{ 
    1612public: 
    17         CParallelTree(ep_lib::MPI_Comm comm); 
    18         ~CParallelTree(); 
     13  CParallelTree(ep_lib::MPI_Comm comm); 
     14  ~CParallelTree(); 
    1915 
    20         void build(vector<Node>& node, vector<Node>& node2); 
     16  void build(vector<Node>& node, vector<Node>& node2); 
    2117 
    22         void routeNodes(vector<int>& route, vector<Node>& nodes, int level = 0); 
    23         void routeIntersections(vector<vector<int> >& route, vector<Node>& nodes, int level = 0); 
     18  void routeNodes(vector<int>& route, vector<Node>& nodes, int level = 0); 
     19  void routeIntersections(vector<vector<int> >& route, vector<Node>& nodes, int level = 0); 
    2420 
    25         int nbLocalElements; 
    26         Elt* localElements; 
     21  int nbLocalElements; 
     22  Elt* localElements; 
    2723 
    28         CTree localTree; 
     24  CTree localTree; 
    2925 
    3026private: 
    31         void updateCirclesForRouting(Coord rootCentre, double rootRadius, int level = 0); 
    32         void buildSampleTreeCascade(vector<Node>& sampleNodes, int level = 0); 
    33         void buildLocalTree(const vector<Node>& node, const vector<int>& route); 
    34         void buildRouteTree(); 
     27  void updateCirclesForRouting(Coord rootCentre, double rootRadius, int level = 0); 
     28  void buildSampleTreeCascade(vector<Node>& sampleNodes, int level = 0); 
     29  void buildLocalTree(const vector<Node>& node, const vector<int>& route); 
     30  void buildRouteTree(); 
    3531 
    36         //CSampleTree sampleTree; 
    37         vector<CSampleTree> treeCascade; // first for sample tree, then for routing tree 
    38         CMPICascade cascade; 
    39         ep_lib::MPI_Comm communicator ; 
     32  //CSampleTree sampleTree; 
     33  vector<CSampleTree> treeCascade; // first for sample tree, then for routing tree 
     34  CMPICascade cascade; 
     35  ep_lib::MPI_Comm communicator ; 
    4036 
    4137}; 
  • XIOS/dev/branch_openmp/extern/remap/src/polyg.cpp

    r1289 r1328  
    77 
    88#include "polyg.hpp" 
    9 #include <stdio.h> 
    109 
    1110namespace sphereRemap { 
     
    4645Coord barycentre(const Coord *x, int n) 
    4746{ 
    48  
    4947        if (n == 0) return ORIGIN; 
    5048        Coord bc = ORIGIN; 
    5149        for (int i = 0; i < n; i++) 
    52         { 
    5350                bc = bc + x[i]; 
    54         }        
    5551        /* both distances can be equal down to roundoff when norm(bc) < mashineepsilon  
    5652           which can occur when weighted with tiny area */ 
    57    
    58          
    59   //assert(squaredist(bc, proj(bc)) <= squaredist(bc, proj(bc * (-1.0))));       
     53 
     54        assert(squaredist(bc, proj(bc)) <= squaredist(bc, proj(bc * (-1.0)))); 
     55        //if (squaredist(bc, proj(bc)) > squaredist(bc, proj(bc * (-1.0)))) return proj(bc * (-1.0)); 
    6056 
    6157        return proj(bc); 
     
    165161{ 
    166162        if (N < 3) 
    167                 return 0; /* polygons with less than three vertices have zero area */ 
     163          return 0; /* polygons with less than three vertices have zero area */ 
    168164        Coord t[3]; 
    169165        t[0] = barycentre(x, N); 
     
    177173                t[1] = x[i]; 
    178174                t[2] = x[ii]; 
    179     double sc=scalarprod(crossprod(t[1] - t[0], t[2] - t[0]), t[0]) ; 
     175                double sc=scalarprod(crossprod(t[1] - t[0], t[2] - t[0]), t[0]) ; 
    180176                assert(sc >= -1e-10); // Error: tri a l'env (wrong orientation) 
    181177                double area_gc = triarea(t[0], t[1], t[2]); 
    182                 //if(area_gc<=0) printf("area_gc = %e\n", area_gc); 
    183178                double area_sc_gc_moon = 0; 
    184179                if (d[i]) /* handle small circle case */ 
     
    188183                        char sgl = (mext > 0) ? -1 : 1; 
    189184                        area_sc_gc_moon = sgl * alun(arcdist(t[1], t[2]), fabs(scalarprod(t[1], pole))); 
    190                         //if(area_sc_gc_moon<=0) printf("area_sc_gc_moon = %e\n", area_sc_gc_moon); 
    191185                        gg_exact = gg_exact + sc_gc_moon_normalintegral(t[1], t[2], pole); 
    192186                } 
    193187                area += area_gc + area_sc_gc_moon; /* for "spherical circle segment" sum triangular part (at) and "small moon" and => account for small circle */ 
    194188                g[i] = barycentre(t, 3) * (area_gc + area_sc_gc_moon); 
    195                 //printf("g[%d] = (%e,%e,%e) * (%e+%e) = (%e,%e,%e) norm = %e\n", i, barycentre(t, 3).x, barycentre(t, 3).y, barycentre(t, 3).z, area_gc,  area_sc_gc_moon, g[i].x, g[i].y, g[i].z, norm(g[i])); 
    196189        } 
    197190        gg = barycentre(g, N); 
  • XIOS/dev/branch_openmp/extern/remap/src/timerRemap.cpp

    r1146 r1328  
    44#include <map> 
    55#include <iostream> 
     6using namespace ep_lib; 
    67 
    78namespace sphereRemap { 
     
    910using namespace std; 
    1011 
    11 map<string,CTimer*> *CTimer::allTimer = 0; 
     12//map<string,CTimer*> CTimer::allTimer; 
     13map<string,CTimer*> *CTimer::allTimer_ptr = 0; 
    1214 
    1315CTimer::CTimer(const string& name_) : name(name_) 
     
    5557CTimer& CTimer::get(const string name) 
    5658{ 
    57         if(allTimer == 0) allTimer = new map<string,CTimer*>; 
    5859        map<string,CTimer*>::iterator it; 
    59         it=(*allTimer).find(name); 
    60         if (it==(*allTimer).end()) it=(*allTimer).insert(pair<string,CTimer*>(name,new CTimer(name))).first; 
     60        if(allTimer_ptr == 0) allTimer_ptr = new map<string,CTimer*>; 
     61        //it=allTimer.find(name); 
     62        it=allTimer_ptr->find(name); 
     63        //if (it==allTimer.end()) it=allTimer.insert(pair<string,CTimer*>(name,new CTimer(name))).first; 
     64        if (it==allTimer_ptr->end()) it=allTimer_ptr->insert(pair<string,CTimer*>(name,new CTimer(name))).first; 
    6165        return *(it->second); 
    6266} 
  • XIOS/dev/branch_openmp/extern/remap/src/timerRemap.hpp

    r1146 r1328  
    2727    void print(void); 
    2828    //static map<string,CTimer*> allTimer; 
    29     static map<string,CTimer*> *allTimer; 
    30     #pragma omp threadprivate(allTimer) 
    31      
     29    static map<string,CTimer*> *allTimer_ptr; 
    3230    static double getTime(void); 
    3331    static CTimer& get(string name); 
  • XIOS/dev/branch_openmp/extern/remap/src/tree.cpp

    r1172 r1328  
    2929void CBasicTree::routeNodes(vector<int>& route, vector<Node>& nodes, int assignLevel) 
    3030{ 
    31   for (int i = 0; i < nodes.size(); i++) 
    32   { 
    33     root->routeNode(&nodes[i], assignLevel); 
    34     route[i] = nodes[i].route; 
    35   } 
     31        for (int i = 0; i < nodes.size(); i++) 
     32        { 
     33                root->routeNode(&nodes[i], assignLevel); 
     34                route[i] = nodes[i].route; 
     35        } 
    3636} 
    3737 
    3838void CBasicTree::routeIntersections(vector<vector<int> >& routes, vector<Node>& nodes) 
    3939{ 
    40   for (int i = 0; i < nodes.size(); i++) 
    41     root->routeIntersection(routes[i], &nodes[i]); 
     40        for (int i = 0; i < nodes.size(); i++) 
     41                root->routeIntersection(routes[i], &nodes[i]); 
    4242} 
    4343 
    4444void CBasicTree::build(vector<Node>& nodes) 
    4545{ 
    46   newRoot(1); 
    47   insertNodes(nodes); 
     46        newRoot(1); 
     47        insertNodes(nodes); 
    4848} 
    4949 
    5050void CBasicTree::output(ostream& flux, int level) 
    5151{ 
    52   root->output(flux,level,0) ; 
     52        root->output(flux,level,0) ; 
    5353} 
    5454void CBasicTree::slim(int nbIts) 
    5555{ 
    56   for (int i = 0; i < nbIts; i++) 
    57   { 
    58     for (int level = root->level - 1; level > 0; level--) 
    59     { 
    60       slim2(root, level); 
    61       ri = 0; 
    62       emptyPool(); 
    63     } 
    64  
    65     for (int level = 2; level < root->level; level++) 
    66     { 
    67       slim2(root, level); 
    68       ri = 0; 
    69       emptyPool(); 
    70     } 
    71   } 
     56        for (int i = 0; i < nbIts; i++) 
     57        { 
     58                for (int level = root->level - 1; level > 0; level--) 
     59                { 
     60                        slim2(root, level); 
     61                        ri = 0; 
     62                        emptyPool(); 
     63                } 
     64 
     65                for (int level = 2; level < root->level; level++) 
     66                { 
     67                        slim2(root, level); 
     68                        ri = 0; 
     69                        emptyPool(); 
     70                } 
     71        } 
    7272} 
    7373 
     
    7676void CBasicTree::insertNode(NodePtr node) 
    7777{ 
    78   node->tree = this; 
    79   increaseLevelSize(0); 
    80   push_back(node); 
    81  
    82   NodePtr q; 
    83   while (pool.size()) 
    84   { 
    85     q = pool.front(); 
    86     pool.pop_front(); 
    87     q = insert(q, root); 
    88     if (ri) 
    89     { 
    90       delete q; 
    91       ri = 0; 
    92     } 
    93   } 
     78        node->tree = this; 
     79        increaseLevelSize(0); 
     80        push_back(node); 
     81 
     82        NodePtr q; 
     83        while (pool.size()) 
     84        { 
     85                q = pool.front(); 
     86                pool.pop_front(); 
     87                q = insert(q, root); 
     88                if (ri) 
     89                { 
     90                        delete q; 
     91                        ri = 0; 
     92                } 
     93        } 
    9494} 
    9595 
    9696void CBasicTree::emptyPool(void) 
    9797{ 
    98   while (pool.size()) 
    99   { 
    100     NodePtr q = pool.front(); 
    101     pool.pop_front(); 
    102     q = insert(q, root); 
    103     if (ri) 
    104     { 
    105       delete q; 
    106       ri = 0; 
    107     } 
    108   } 
     98        while (pool.size()) 
     99        { 
     100                NodePtr q = pool.front(); 
     101                pool.pop_front(); 
     102                q = insert(q, root); 
     103                if (ri) 
     104                { 
     105                        delete q; 
     106                        ri = 0; 
     107                } 
     108        } 
    109109} 
    110110 
     
    142142        root->parent = 0; 
    143143        root->leafCount = 0; 
    144         // initialize root node on the sphere 
    145         root->centre.x=1 ;  
    146         root->centre.y=0 ;  
    147         root->centre.z=0 ;  
     144// initialize root node on the sphere 
     145  root->centre.x=1 ; root->centre.y=0 ; root->centre.z=0 ;  
    148146        root->radius = 0.; 
    149147        root->reinserted = false; 
  • XIOS/dev/branch_openmp/extern/remap/src/tree.hpp

    r1172 r1328  
    1414class CBasicTree 
    1515{ 
    16   public: 
     16public: 
    1717 
    18   NodePtr root; /* The main tree is stored as Nodes which can be reached through traversal starting here */ 
    19   NodePtr ref; // FIXME this reference, set by a node is odd, try to remove 
    20   int ri; /** this is set to one by a node in case of reinsertion */ 
    21   vector<int> levelSize; /** e.g. levelSize[0] == leafs.size() */ 
    22   vector<Node> leafs; /** leafs are stored in vector for easy access and rest of the tree nodes as separate allocations, only reachable through tree traversal */ 
     18        NodePtr root; /* The main tree is stored as Nodes which can be reached through traversal starting here */ 
     19        NodePtr ref; // FIXME this reference, set by a node is odd, try to remove 
     20        int ri; /** this is set to one by a node in case of reinsertion */ 
     21        vector<int> levelSize; /** e.g. levelSize[0] == leafs.size() */ 
     22        vector<Node> leafs; /** leafs are stored in vector for easy access and rest of the tree nodes as separate allocations, only reachable through tree traversal */ 
    2323 
    24   CBasicTree() : ri(0), levelSize(MAX_LEVEL_SIZE), root(NULL), isAssignedLevel(false), okSplit(true), isActiveOkSplit(false) {}  
    25   ~CBasicTree();  
    26   void build(vector<Node>& nodes); 
    27   void slim(int nbIts = 1); 
    28   virtual void insertNodes(vector<Node>& node) = 0; 
     24        CBasicTree() : ri(0), levelSize(MAX_LEVEL_SIZE), root(NULL), isAssignedLevel(false), okSplit(true), isActiveOkSplit(false) {}  
     25        ~CBasicTree();  
     26        void build(vector<Node>& nodes); 
     27        void slim(int nbIts = 1); 
     28        virtual void insertNodes(vector<Node>& node) = 0; 
    2929 
    30   void routeNodes(vector<int>& route, vector<Node>& nodes, int assignLevel); 
    31   void routeIntersections(vector<vector<int> >& route, vector<Node>& nodes); 
     30        void routeNodes(vector<int>& route, vector<Node>& nodes, int assignLevel); 
     31        void routeIntersections(vector<vector<int> >& route, vector<Node>& nodes); 
    3232 
    33   void push_back(NodePtr node); 
    34   void push_front(NodePtr node); 
    35   void increaseLevelSize(int level); 
    36   void decreaseLevelSize(int level); 
    37   void newRoot(int level); 
    38   void insertNode(NodePtr node); 
     33        void push_back(NodePtr node); 
     34        void push_front(NodePtr node); 
     35        void increaseLevelSize(int level); 
     36        void decreaseLevelSize(int level); 
     37        void newRoot(int level); 
     38        void insertNode(NodePtr node); 
    3939  void output(ostream& flux, int level) ; 
    4040 
    41   int keepNodes; 
     41        int keepNodes; 
    4242  bool isAssignedLevel ;  
    4343  int assignLevel; 
     
    5050 
    5151   
    52   private: 
    53   deque<NodePtr > pool; 
     52private: 
     53        deque<NodePtr > pool; 
    5454         
    5555  bool okSplit ; 
    5656  
    57   protected: 
     57protected: 
    5858  void emptyPool(); 
    59   CBasicTree(int keepNodes_, int assignLevel_) : ri(0), levelSize(MAX_LEVEL_SIZE), root(NULL), keepNodes(keepNodes_), assignLevel(assignLevel_), isAssignedLevel(true),  
    60                                                  okSplit(true), isActiveOkSplit(false) {}  
     59  CBasicTree(int keepNodes_, int assignLevel_) : ri(0), levelSize(MAX_LEVEL_SIZE), root(NULL), keepNodes(keepNodes_), assignLevel(assignLevel_), isAssignedLevel(true), okSplit(true), isActiveOkSplit(false) {}  
    6160}; 
    6261 
    6362class CTree : public CBasicTree 
    6463{ 
    65   public: 
    66   void insertNodes(vector<Node>& nodes); 
     64public: 
     65        void insertNodes(vector<Node>& nodes); 
    6766}; 
    6867 
     
    7069{ 
    7170 
    72   public: 
    73   CSampleTree(int keepNodes_, int assignLevel_) : CBasicTree(keepNodes_,assignLevel_) {} 
     71public: 
     72        CSampleTree(int keepNodes_, int assignLevel_) : CBasicTree(keepNodes_,assignLevel_) {} 
    7473  void slimAssignedLevel() ; 
    7574  void removeExtraNode(void) ; 
    76   void insertNodes(vector<Node>& nodes); 
     75        void insertNodes(vector<Node>& nodes); 
    7776}; 
    7877 
Note: See TracChangeset for help on using the changeset viewer.