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

dev_omp

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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} 
Note: See TracChangeset for help on using the changeset viewer.