Ignore:
Timestamp:
01/31/19 12:12:52 (5 years ago)
Author:
yushan
Message:

branch merged with trunk @1645. arch file (ep&mpi) added for ADA

File:
1 edited

Legend:

Unmodified
Added
Removed

  • XIOS/dev/dev_trunk_omp/extern/remap/src/mapper.cpp

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