Changeset 1614


Ignore:
Timestamp:
11/26/18 10:16:06 (6 years ago)
Author:
ymipsl
Message:

Interpolation : enhancement : you have now the possibility to give computed area to the remaper for better global conservation.

YM

Location:
XIOS/trunk/extern/remap/src
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • XIOS/trunk/extern/remap/src/elt.hpp

    r1158 r1614  
    4848        int n; /* number of vertices */ 
    4949        double area; 
     50  double given_area ; 
    5051        Coord x; /* barycentre */ 
    5152}; 
    … …  
    8081                n    = rhs.n; 
    8182                area = rhs.area; 
     83                given_area = rhs.given_area; 
    8284                x    = rhs.x; 
    8385                val  = rhs.val; 

  • XIOS/trunk/extern/remap/src/libmapper.cpp

    r694 r1614  
    4141        assert(n_cell_dst >= 4); 
    4242        assert(1 <= order && order <= 2); 
    43  
     43  double* src_area=NULL ; 
     44  double* dst_area=NULL ; 
    4445  mapper = new Mapper(MPI_COMM_WORLD); 
    4546  mapper->setVerbosity(PROGRESS) ; 
    46   mapper->setSourceMesh(src_bounds_lon, src_bounds_lat, n_vert_per_cell_src, n_cell_src, src_pole ) ; 
    47   mapper->setTargetMesh(dst_bounds_lon, dst_bounds_lat, n_vert_per_cell_dst, n_cell_dst, dst_pole ) ; 
     47  mapper->setSourceMesh(src_bounds_lon, src_bounds_lat, src_area, n_vert_per_cell_src, n_cell_src, src_pole ) ; 
     48  mapper->setTargetMesh(dst_bounds_lon, dst_bounds_lat, dst_area, n_vert_per_cell_dst, n_cell_dst, dst_pole ) ; 
    4849 
    4950/* 

  • XIOS/trunk/extern/remap/src/mapper.cpp

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

  • XIOS/trunk/extern/remap/src/mapper.hpp

    r1158 r1614  
    2222       void setVerbosity(verbosity v) {verbose=v ;} 
    2323 
    24        void setSourceMesh(const double* boundsLon, const double* boundsLat, int nVertex, int nbCells, const double* pole, const long int* globalId=NULL) ; 
    25        void setTargetMesh(const double* boundsLon, const double* boundsLat, int nVertex, int nbCells, const double* pole, const long int* globalId=NULL) ; 
     24       void setSourceMesh(const double* boundsLon, const double* boundsLat, const double* area, int nVertex, int nbCells, const double* pole, const long int* globalId=NULL) ; 
     25       void setTargetMesh(const double* boundsLon, const double* boundsLat, const double* area, int nVertex, int nbCells, const double* pole, const long int* globalId=NULL) ; 
    2626       void setSourceValue(const double* val) ; 
    2727       void getTargetValue(double* val) ; 

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

    r1579 r1614  
    218218int packedPolygonSize(const Elt& e) 
    219219{ 
    220   return sizeof(e.id) + sizeof(e.src_id) + sizeof(e.x) + sizeof(e.val) + 
     220  return sizeof(e.id) + sizeof(e.src_id) + sizeof(e.x) + sizeof(e.val)  + sizeof(e.given_area)+ 
    221221         sizeof(e.n) + e.n*(sizeof(double)+sizeof(Coord)); 
    222222} 
    … …  
    235235  pos += sizeof(e.val); 
    236236 
     237  *((double*) &(buffer[pos])) = e.given_area; 
     238  pos += sizeof(e.val); 
     239 
    237240  *((int *) & (buffer[pos])) = e.n; 
    238241  pos += sizeof(e.n); 
    … …  
    262265  pos += sizeof(double); 
    263266 
     267  e.given_area = *((double *) & (buffer[pos])); 
     268  pos += sizeof(double); 
     269 
    264270  e.n = *((int *) & (buffer[pos])); 
    265271  pos += sizeof(int); 
    … …  
    291297    *((double *) &(buffer[pos])) = e.area; 
    292298    pos += sizeof(double); 
     299 
    293300 
    294301    *((GloId *) &(buffer[pos])) = (*it)->id; 
    … …  
    322329    pos += sizeof(double); 
    323330 
     331 
    324332    Polyg *polygon = new Polyg; 
    325333 
Note: See TracChangeset for help on using the changeset viewer.