Changeset 1619


Ignore:
Timestamp:
11/30/18 17:37:26 (9 months ago)
Author:
yushan
Message:

branch_openmp merged and NOT tested with trunk r1618

Location:
XIOS/dev/dev_trunk_omp
Files:
12 edited

Legend:

Unmodified
Added
Removed
  • XIOS/dev/dev_trunk_omp/extern/remap/src/elt.hpp

    r1158 r1619  
    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/dev/dev_trunk_omp/extern/remap/src/libmapper.cpp

    r1602 r1619  
    4343        assert(n_cell_dst >= 4); 
    4444        assert(1 <= order && order <= 2); 
    45  
     45  double* src_area=NULL ; 
     46  double* dst_area=NULL ; 
    4647  mapper = new Mapper(MPI_COMM_WORLD); 
    4748  mapper->setVerbosity(PROGRESS) ; 
    48   mapper->setSourceMesh(src_bounds_lon, src_bounds_lat, n_vert_per_cell_src, n_cell_src, src_pole ) ; 
    49   mapper->setTargetMesh(dst_bounds_lon, dst_bounds_lat, n_vert_per_cell_dst, n_cell_dst, dst_pole ) ; 
     49  mapper->setSourceMesh(src_bounds_lon, src_bounds_lat, src_area, n_vert_per_cell_src, n_cell_src, src_pole ) ; 
     50  mapper->setTargetMesh(dst_bounds_lon, dst_bounds_lat, dst_area, n_vert_per_cell_dst, n_cell_dst, dst_pole ) ; 
    5051 
    5152/* 
  • XIOS/dev/dev_trunk_omp/extern/remap/src/mapper.cpp

    r1602 r1619  
    2929void cptOffsetsFromLengths(const int *lengths, int *offsets, int sz) 
    3030{ 
    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) 
     31  offsets[0] = 0; 
     32  for (int i = 1; i < sz; i++) 
     33    offsets[i] = offsets[i-1] + lengths[i-1]; 
     34} 
     35 
     36 
     37void Mapper::setSourceMesh(const double* boundsLon, const double* boundsLat, const double* area, int nVertex, int nbCells, const double* pole, const long int* globalId) 
    3838{ 
    3939  srcGrid.pole = Coord(pole[0], pole[1], pole[2]); 
    4040 
    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); 
     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); 
    4747  sourceGlobalId.resize(nbCells) ; 
    4848 
     
    5757  else sourceGlobalId.assign(globalId,globalId+nbCells); 
    5858 
    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; 
     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; 
    6565    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) 
     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    if (area!=NULL) sourceElements[i].given_area=area[i] ; 
     70    else sourceElements[i].given_area=sourceElements[i].area ; 
     71  } 
     72 
     73} 
     74 
     75void Mapper::setTargetMesh(const double* boundsLon, const double* boundsLat, const double* area, int nVertex, int nbCells, const double* pole, const long int* globalId) 
    7476{ 
    7577  tgtGrid.pole = Coord(pole[0], pole[1], pole[2]); 
    7678 
    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); 
     79  int mpiRank, mpiSize; 
     80  MPI_Comm_rank(communicator, &mpiRank); 
     81  MPI_Comm_size(communicator, &mpiSize); 
     82 
     83  targetElements.reserve(nbCells); 
     84  targetMesh.reserve(nbCells); 
    8385 
    8486  targetGlobalId.resize(nbCells) ; 
     
    9395  else targetGlobalId.assign(globalId,globalId+nbCells); 
    9496 
    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         } 
     97  for (int i = 0; i < nbCells; i++) 
     98  { 
     99    int offs = i*nVertex; 
     100    Elt elt(boundsLon + offs, boundsLat + offs, nVertex); 
     101    targetElements.push_back(elt); 
     102    targetMesh.push_back(Node(elt.x, cptRadius(elt), &sourceElements.back())); 
     103    cptEltGeom(targetElements[i], Coord(pole[0], pole[1], pole[2])); 
     104    if (area!=NULL) targetElements[i].given_area=area[i] ; 
     105    else targetElements[i].given_area=targetElements[i].area ; 
     106  } 
    103107 
    104108 
     
    119123vector<double> Mapper::computeWeights(int interpOrder, bool renormalize, bool quantity) 
    120124{ 
    121         vector<double> timings; 
    122         int mpiSize, mpiRank; 
    123         MPI_Comm_size(communicator, &mpiSize); 
    124         MPI_Comm_rank(communicator, &mpiRank); 
     125  vector<double> timings; 
     126  int mpiSize, mpiRank; 
     127  MPI_Comm_size(communicator, &mpiSize); 
     128  MPI_Comm_rank(communicator, &mpiRank); 
    125129 
    126130  this->buildSSTree(sourceMesh, targetMesh); 
    127131 
    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 
     132  if (mpiRank == 0 && verbose) cout << "Computing intersections ..." << endl; 
     133  double tic = cputime(); 
     134  computeIntersection(&targetElements[0], targetElements.size()); 
     135  timings.push_back(cputime() - tic); 
     136 
     137  tic = cputime(); 
     138  if (interpOrder == 2) { 
     139    if (mpiRank == 0 && verbose) cout << "Computing grads ..." << endl; 
     140    buildMeshTopology(); 
     141    computeGrads(); 
     142  } 
     143  timings.push_back(cputime() - tic); 
     144 
     145  /* Prepare computation of weights */ 
     146  /* compute number of intersections which for the first order case 
    143147           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]; 
     148  int nIntersections = 0; 
     149  for (int j = 0; j < targetElements.size(); j++) 
     150  { 
     151    Elt &elt = targetElements[j]; 
     152    for (list<Polyg*>::iterator it = elt.is.begin(); it != elt.is.end(); it++) 
     153      nIntersections++; 
     154  } 
     155  /* overallocate for NMAX neighbours for each elements */ 
     156  remapMatrix = new double[nIntersections*NMAX]; 
     157  srcAddress = new int[nIntersections*NMAX]; 
     158  srcRank = new int[nIntersections*NMAX]; 
     159  dstAddress = new int[nIntersections*NMAX]; 
    156160  sourceWeightId =new long[nIntersections*NMAX]; 
    157161  targetWeightId =new long[nIntersections*NMAX]; 
    158162 
    159163 
    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); 
     164  if (mpiRank == 0 && verbose) cout << "Remapping..." << endl; 
     165  tic = cputime(); 
     166  nWeights = remap(&targetElements[0], targetElements.size(), interpOrder, renormalize, quantity); 
     167  timings.push_back(cputime() - tic); 
    164168 
    165169  for (int i = 0; i < targetElements.size(); i++) targetElements[i].delete_intersections(); 
    166170 
    167         return timings; 
     171  return timings; 
    168172} 
    169173 
     
    176180int Mapper::remap(Elt *elements, int nbElements, int order, bool renormalize, bool quantity) 
    177181{ 
    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); 
     182  int mpiSize, mpiRank; 
     183  MPI_Comm_size(communicator, &mpiSize); 
     184  MPI_Comm_rank(communicator, &mpiRank); 
     185 
     186  /* create list of intersections (super mesh elements) for each rank */ 
     187  multimap<int, Polyg *> *elementList = new multimap<int, Polyg *>[mpiSize]; 
     188  for (int j = 0; j < nbElements; j++) 
     189  { 
     190    Elt& e = elements[j]; 
     191    for (list<Polyg *>::iterator it = e.is.begin(); it != e.is.end(); it++) 
     192      elementList[(*it)->id.rank].insert(pair<int, Polyg *>((*it)->id.ind, *it)); 
     193  } 
     194 
     195  int *nbSendElement = new int[mpiSize]; 
     196  int **sendElement = new int*[mpiSize]; /* indices of elements required from other rank */ 
     197  double **recvValue = new double*[mpiSize]; 
     198  double **recvArea = new double*[mpiSize]; 
     199  double **recvGivenArea = new double*[mpiSize]; 
     200  Coord **recvGrad = new Coord*[mpiSize]; 
     201  GloId **recvNeighIds = new GloId*[mpiSize]; /* ids of the of the source neighbours which also contribute through gradient */ 
     202  for (int rank = 0; rank < mpiSize; rank++) 
     203  { 
     204    /* get size for allocation */ 
     205    int last = -1; /* compares unequal to any index */ 
     206    int index = -1; /* increased to starting index 0 in first iteration */ 
     207    for (multimap<int, Polyg *>::iterator it = elementList[rank].begin(); it != elementList[rank].end(); ++it) 
     208    { 
     209      if (last != it->first) 
     210        index++; 
     211      (it->second)->id.ind = index; 
     212      last = it->first; 
     213    } 
     214    nbSendElement[rank] = index + 1; 
     215 
     216    /* if size is non-zero allocate and collect indices of elements on other ranks that we intersect */ 
     217    if (nbSendElement[rank] > 0) 
     218    { 
     219      sendElement[rank] = new int[nbSendElement[rank]]; 
     220      recvValue[rank]   = new double[nbSendElement[rank]]; 
     221      recvArea[rank]    = new double[nbSendElement[rank]]; 
     222      recvGivenArea[rank] = new double[nbSendElement[rank]]; 
     223      if (order == 2) 
     224      { 
     225        recvNeighIds[rank] = new GloId[nbSendElement[rank]*(NMAX+1)]; 
     226        recvGrad[rank]    = new Coord[nbSendElement[rank]*(NMAX+1)]; 
     227      } 
     228      else 
     229        recvNeighIds[rank] = new GloId[nbSendElement[rank]]; 
     230 
     231      last = -1; 
     232      index = -1; 
     233      for (multimap<int, Polyg *>::iterator it = elementList[rank].begin(); it != elementList[rank].end(); ++it) 
     234      { 
     235        if (last != it->first) 
     236          index++; 
     237        sendElement[rank][index] = it->first; 
     238        last = it->first; 
     239      } 
     240    } 
     241  } 
     242 
     243  /* communicate sizes of source elements to be sent (index lists and later values and gradients) */ 
     244  int *nbRecvElement = new int[mpiSize]; 
     245  MPI_Alltoall(nbSendElement, 1, MPI_INT, nbRecvElement, 1, MPI_INT, communicator); 
     246 
     247  /* communicate indices of source elements on other ranks whoes value and gradient we need (since intersection) */ 
     248  int nbSendRequest = 0; 
     249  int nbRecvRequest = 0; 
     250  int **recvElement = new int*[mpiSize]; 
     251  double **sendValue = new double*[mpiSize]; 
     252  double **sendArea = new double*[mpiSize]; 
     253  double **sendGivenArea = new double*[mpiSize]; 
     254  Coord **sendGrad = new Coord*[mpiSize]; 
     255  GloId **sendNeighIds = new GloId*[mpiSize]; 
     256  MPI_Request *sendRequest = new MPI_Request[5*mpiSize]; 
     257  MPI_Request *recvRequest = new MPI_Request[5*mpiSize]; 
     258  for (int rank = 0; rank < mpiSize; rank++) 
     259  { 
     260    if (nbSendElement[rank] > 0) 
     261    { 
     262      MPI_Issend(sendElement[rank], nbSendElement[rank], MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]); 
     263      nbSendRequest++; 
     264    } 
     265 
     266    if (nbRecvElement[rank] > 0) 
     267    { 
     268      recvElement[rank] = new int[nbRecvElement[rank]]; 
     269      sendValue[rank]   = new double[nbRecvElement[rank]]; 
     270      sendArea[rank]   = new double[nbRecvElement[rank]]; 
     271      sendGivenArea[rank] = new double[nbRecvElement[rank]]; 
     272      if (order == 2) 
     273      { 
     274        sendNeighIds[rank] = new GloId[nbRecvElement[rank]*(NMAX+1)]; 
     275        sendGrad[rank]    = new Coord[nbRecvElement[rank]*(NMAX+1)]; 
     276      } 
     277      else 
     278      { 
     279        sendNeighIds[rank] = new GloId[nbRecvElement[rank]]; 
     280      } 
     281      MPI_Irecv(recvElement[rank], nbRecvElement[rank], MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
     282      nbRecvRequest++; 
     283    } 
     284  } 
     285  MPI_Status *status = new MPI_Status[5*mpiSize]; 
     286   
     287  MPI_Waitall(nbSendRequest, sendRequest, status); 
    280288        MPI_Waitall(nbRecvRequest, recvRequest, status); 
    281289 
    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; 
     290  /* for all indices that have been received from requesting ranks: pack values and gradients, then send */ 
     291  nbSendRequest = 0; 
     292  nbRecvRequest = 0; 
     293  for (int rank = 0; rank < mpiSize; rank++) 
     294  { 
     295    if (nbRecvElement[rank] > 0) 
     296    { 
     297      int jj = 0; // jj == j if no weight writing 
     298      for (int j = 0; j < nbRecvElement[rank]; j++) 
     299      { 
     300        sendValue[rank][j] = sstree.localElements[recvElement[rank][j]].val; 
     301        sendArea[rank][j] = sstree.localElements[recvElement[rank][j]].area; 
     302        sendGivenArea[rank][j] = sstree.localElements[recvElement[rank][j]].given_area; 
     303        if (order == 2) 
     304        { 
     305          sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].grad; 
    297306//          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]; 
     307          sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].src_id; 
     308          jj++; 
     309          for (int i = 0; i < NMAX; i++) 
     310          { 
     311            sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].gradNeigh[i]; 
    303312//            cout<<"grad  "<<jj<<"  "<<sendGrad[rank][jj]<<" "<<sstree.localElements[recvElement[rank][j]].grad<<endl ; 
    304313            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]); 
     314            jj++; 
     315          } 
     316        } 
     317        else 
     318          sendNeighIds[rank][j] = sstree.localElements[recvElement[rank][j]].src_id; 
     319      } 
     320      MPI_Issend(sendValue[rank],  nbRecvElement[rank], MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]); 
     321      nbSendRequest++; 
     322      MPI_Issend(sendArea[rank],  nbRecvElement[rank], MPI_DOUBLE, rank, 1, communicator, &sendRequest[nbSendRequest]); 
     323      nbSendRequest++; 
     324      MPI_Issend(sendGivenArea[rank],  nbRecvElement[rank], MPI_DOUBLE, rank, 4, communicator, &sendRequest[nbSendRequest]); 
     325      nbSendRequest++; 
     326      if (order == 2) 
     327      { 
     328        MPI_Issend(sendGrad[rank], 3*nbRecvElement[rank]*(NMAX+1), MPI_DOUBLE, rank, 2, communicator, &sendRequest[nbSendRequest]); 
     329        nbSendRequest++; 
     330        MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank]*(NMAX+1), MPI_INT, rank, 3, communicator, &sendRequest[nbSendRequest]); 
    321331//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]); 
     332        nbSendRequest++; 
     333      } 
     334      else 
     335      { 
     336        MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank], MPI_INT, rank, 5, communicator, &sendRequest[nbSendRequest]); 
    327337//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]); 
     338        nbSendRequest++; 
     339      } 
     340    } 
     341    if (nbSendElement[rank] > 0) 
     342    { 
     343      MPI_Irecv(recvValue[rank],  nbSendElement[rank], MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
     344      nbRecvRequest++; 
     345      MPI_Irecv(recvArea[rank],  nbSendElement[rank], MPI_DOUBLE, rank, 1, communicator, &recvRequest[nbRecvRequest]); 
     346      nbRecvRequest++; 
     347      MPI_Irecv(recvGivenArea[rank],  nbSendElement[rank], MPI_DOUBLE, rank, 4, communicator, &recvRequest[nbRecvRequest]); 
     348      nbRecvRequest++; 
     349      if (order == 2) 
     350      { 
     351        MPI_Irecv(recvGrad[rank], 3*nbSendElement[rank]*(NMAX+1), 
     352            MPI_DOUBLE, rank, 2, communicator, &recvRequest[nbRecvRequest]); 
     353        nbRecvRequest++; 
     354        MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank]*(NMAX+1), MPI_INT, rank, 3, communicator, &recvRequest[nbRecvRequest]); 
    343355//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]); 
     356        nbRecvRequest++; 
     357      } 
     358      else 
     359      { 
     360        MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank], MPI_INT, rank, 5, communicator, &recvRequest[nbRecvRequest]); 
    349361//ym  --> attention taille GloId 
    350                                 nbRecvRequest++; 
    351                         } 
    352                 } 
    353         } 
     362        nbRecvRequest++; 
     363      } 
     364    } 
     365  } 
    354366         
    355367        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; 
     368  MPI_Waitall(nbRecvRequest, recvRequest, status); 
     369   
     370 
     371  /* now that all values and gradients are available use them to computed interpolated values on target 
     372     and also to compute weights */ 
     373  int i = 0; 
     374  for (int j = 0; j < nbElements; j++) 
     375  { 
     376    Elt& e = elements[j]; 
     377 
     378    /* since for the 2nd order case source grid elements can contribute to a destination grid element over several "paths" 
     379       (step1: gradient is computed using neighbours on same grid, step2: intersection uses several elements on other grid) 
     380       accumulate them so that there is only one final weight between two elements */ 
     381    map<GloId,double> wgt_map; 
     382 
     383    /* for destination element `e` loop over all intersetions/the corresponding source elements */ 
     384    for (list<Polyg *>::iterator it = e.is.begin(); it != e.is.end(); it++) 
     385    { 
     386      /* it is the intersection element, so it->x and it->area are barycentre and area of intersection element (super mesh) 
     387      but it->id is id of the source element that it intersects */ 
     388      int n1 = (*it)->id.ind; 
     389      int rank = (*it)->id.rank; 
     390      double fk = recvValue[rank][n1]; 
     391      double srcArea = recvArea[rank][n1]; 
     392      double srcGivenArea = recvGivenArea[rank][n1]; 
     393      double w = (*it)->area; 
    381394      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                 } 
     395      else w=w*srcGivenArea/srcArea*e.area/e.given_area ; 
     396 
     397      /* first order: src value times weight (weight = supermesh area), later divide by target area */ 
     398      int kk = (order == 2) ? n1 * (NMAX + 1) : n1; 
     399      GloId neighID = recvNeighIds[rank][kk]; 
     400      wgt_map[neighID] += w; 
     401 
     402      if (order == 2) 
     403      { 
     404        for (int k = 0; k < NMAX+1; k++) 
     405        { 
     406          int kk = n1 * (NMAX + 1) + k; 
     407          GloId neighID = recvNeighIds[rank][kk]; 
     408          if (neighID.ind != -1)  wgt_map[neighID] += w * scalarprod(recvGrad[rank][kk], (*it)->x); 
     409        } 
     410 
     411      } 
     412    } 
    399413 
    400414    double renorm=0; 
    401415    if (renormalize)  
    402       for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++) renorm+=it->second / e.area; 
     416    { 
     417      if (quantity) for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++) renorm+=it->second ; 
     418      else for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++) renorm+=it->second / e.area; 
     419    } 
    403420    else renorm=1. ; 
    404421 
    405422    for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++) 
    406                 { 
     423    { 
    407424      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; 
     425      else this->remapMatrix[i] = (it->second / e.area) / renorm; 
     426      this->srcAddress[i] = it->first.ind; 
     427      this->srcRank[i] = it->first.rank; 
     428      this->dstAddress[i] = j; 
    412429      this->sourceWeightId[i]= it->first.globalId ; 
    413430      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; 
     431      i++; 
     432    } 
     433  } 
     434 
     435  /* free all memory allocated in this function */ 
     436  for (int rank = 0; rank < mpiSize; rank++) 
     437  { 
     438    if (nbSendElement[rank] > 0) 
     439    { 
     440      delete[] sendElement[rank]; 
     441      delete[] recvValue[rank]; 
     442      delete[] recvArea[rank]; 
     443      delete[] recvGivenArea[rank]; 
     444      if (order == 2) 
     445      { 
     446        delete[] recvGrad[rank]; 
     447      } 
     448      delete[] recvNeighIds[rank]; 
     449    } 
     450    if (nbRecvElement[rank] > 0) 
     451    { 
     452      delete[] recvElement[rank]; 
     453      delete[] sendValue[rank]; 
     454      delete[] sendArea[rank]; 
     455      delete[] sendGivenArea[rank]; 
     456      if (order == 2) 
     457        delete[] sendGrad[rank]; 
     458      delete[] sendNeighIds[rank]; 
     459    } 
     460  } 
     461  delete[] status; 
     462  delete[] sendRequest; 
     463  delete[] recvRequest; 
     464  delete[] elementList; 
     465  delete[] nbSendElement; 
     466  delete[] nbRecvElement; 
     467  delete[] sendElement; 
     468  delete[] recvElement; 
     469  delete[] sendValue; 
     470  delete[] recvValue; 
     471  delete[] sendGrad; 
     472  delete[] recvGrad; 
     473  delete[] sendNeighIds; 
     474  delete[] recvNeighIds; 
     475  return i; 
    457476} 
    458477 
    459478void Mapper::computeGrads() 
    460479{ 
    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); 
     480  /* array of pointers to collect local elements and elements received from other cpu */ 
     481  vector<Elt*> globalElements(sstree.nbLocalElements + nbNeighbourElements); 
     482  int index = 0; 
     483  for (int i = 0; i < sstree.nbLocalElements; i++, index++) 
     484    globalElements[index] = &(sstree.localElements[i]); 
     485  for (int i = 0; i < nbNeighbourElements; i++, index++) 
     486    globalElements[index] = &neighbourElements[i]; 
     487 
     488  update_baryc(sstree.localElements, sstree.nbLocalElements); 
     489  computeGradients(&globalElements[0], sstree.nbLocalElements); 
    471490} 
    472491 
     
    475494void Mapper::buildMeshTopology() 
    476495{ 
    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(); 
     496  int mpiSize, mpiRank; 
     497  MPI_Comm_size(communicator, &mpiSize); 
     498  MPI_Comm_rank(communicator, &mpiRank); 
     499 
     500  vector<Node> *routingList = new vector<Node>[mpiSize]; 
     501  vector<vector<int> > routes(sstree.localTree.leafs.size()); 
     502 
     503  sstree.routeIntersections(routes, sstree.localTree.leafs); 
     504 
     505  for (int i = 0; i < routes.size(); ++i) 
     506    for (int k = 0; k < routes[i].size(); ++k) 
     507      routingList[routes[i][k]].push_back(sstree.localTree.leafs[i]); 
     508  routingList[mpiRank].clear(); 
     509 
     510 
     511  CMPIRouting mpiRoute(communicator); 
     512  mpiRoute.init(routes); 
     513  int nRecv = mpiRoute.getTotalSourceElement(); 
    495514// cout << mpiRank << " NRECV " << nRecv << "(" << routes.size() << ")"<< endl; 
    496515 
    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 
     516  int *nbSendNode = new int[mpiSize]; 
     517  int *nbRecvNode = new int[mpiSize]; 
     518  int *sendMessageSize = new int[mpiSize]; 
     519  int *recvMessageSize = new int[mpiSize]; 
     520 
     521  for (int rank = 0; rank < mpiSize; rank++) 
     522  { 
     523    nbSendNode[rank] = routingList[rank].size(); 
     524    sendMessageSize[rank] = 0; 
     525    for (size_t j = 0; j < routingList[rank].size(); j++) 
     526    { 
     527      Elt *elt = (Elt *) (routingList[rank][j].data); 
     528      sendMessageSize[rank] += packedPolygonSize(*elt); 
     529    } 
     530  } 
     531 
     532  MPI_Alltoall(nbSendNode, 1, MPI_INT, nbRecvNode, 1, MPI_INT, communicator); 
     533  MPI_Alltoall(sendMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator); 
     534 
     535  char **sendBuffer = new char*[mpiSize]; 
     536  char **recvBuffer = new char*[mpiSize]; 
     537  int *pos = new int[mpiSize]; 
     538 
     539  for (int rank = 0; rank < mpiSize; rank++) 
     540  { 
     541    if (nbSendNode[rank] > 0) sendBuffer[rank] = new char[sendMessageSize[rank]]; 
     542    if (nbRecvNode[rank] > 0) recvBuffer[rank] = new char[recvMessageSize[rank]]; 
     543  } 
     544 
     545  for (int rank = 0; rank < mpiSize; rank++) 
     546  { 
     547    pos[rank] = 0; 
     548    for (size_t j = 0; j < routingList[rank].size(); j++) 
     549    { 
     550      Elt *elt = (Elt *) (routingList[rank][j].data); 
     551      packPolygon(*elt, sendBuffer[rank], pos[rank]); 
     552    } 
     553  } 
     554  delete [] routingList; 
     555 
     556 
     557  int nbSendRequest = 0; 
     558  int nbRecvRequest = 0; 
     559  MPI_Request *sendRequest = new MPI_Request[mpiSize]; 
     560  MPI_Request *recvRequest = new MPI_Request[mpiSize]; 
     561  MPI_Status  *status      = new MPI_Status[mpiSize]; 
     562 
     563  for (int rank = 0; rank < mpiSize; rank++) 
     564  { 
     565    if (nbSendNode[rank] > 0) 
     566    { 
     567      MPI_Issend(sendBuffer[rank], sendMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]); 
     568      nbSendRequest++; 
     569    } 
     570    if (nbRecvNode[rank] > 0) 
     571    { 
     572      MPI_Irecv(recvBuffer[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
     573      nbRecvRequest++; 
     574    } 
     575  } 
     576 
     577  MPI_Waitall(nbRecvRequest, recvRequest, status); 
     578  MPI_Waitall(nbSendRequest, sendRequest, status); 
     579 
     580  for (int rank = 0; rank < mpiSize; rank++) 
     581    if (nbSendNode[rank] > 0) delete [] sendBuffer[rank]; 
     582  delete [] sendBuffer; 
     583 
     584  char **sendBuffer2 = new char*[mpiSize]; 
     585  char **recvBuffer2 = new char*[mpiSize]; 
     586 
     587  for (int rank = 0; rank < mpiSize; rank++) 
     588  { 
     589    nbSendNode[rank] = 0; 
     590    sendMessageSize[rank] = 0; 
     591 
     592    if (nbRecvNode[rank] > 0) 
     593    { 
     594      set<NodePtr> neighbourList; 
     595      pos[rank] = 0; 
     596      for (int j = 0; j < nbRecvNode[rank]; j++) 
     597      { 
     598        Elt elt; 
     599        unpackPolygon(elt, recvBuffer[rank], pos[rank]); 
     600        Node node(elt.x, cptRadius(elt), &elt); 
     601        findNeighbour(sstree.localTree.root, &node, neighbourList); 
     602      } 
     603      nbSendNode[rank] = neighbourList.size(); 
     604      for (set<NodePtr>::iterator it = neighbourList.begin(); it != neighbourList.end(); it++) 
     605      { 
     606        Elt *elt = (Elt *) ((*it)->data); 
     607        sendMessageSize[rank] += packedPolygonSize(*elt); 
     608      } 
     609 
     610      sendBuffer2[rank] = new char[sendMessageSize[rank]]; 
     611      pos[rank] = 0; 
     612 
     613      for (set<NodePtr>::iterator it = neighbourList.begin(); it != neighbourList.end(); it++) 
     614      { 
     615        Elt *elt = (Elt *) ((*it)->data); 
     616        packPolygon(*elt, sendBuffer2[rank], pos[rank]); 
     617      } 
     618    } 
     619  } 
     620  for (int rank = 0; rank < mpiSize; rank++) 
     621    if (nbRecvNode[rank] > 0) delete [] recvBuffer[rank]; 
     622  delete [] recvBuffer; 
     623 
     624 
     625  MPI_Barrier(communicator); 
     626  MPI_Alltoall(nbSendNode, 1, MPI_INT, nbRecvNode, 1, MPI_INT, communicator); 
     627  MPI_Alltoall(sendMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator); 
     628 
     629  for (int rank = 0; rank < mpiSize; rank++) 
     630    if (nbRecvNode[rank] > 0) recvBuffer2[rank] = new char[recvMessageSize[rank]]; 
     631 
     632  nbSendRequest = 0; 
     633  nbRecvRequest = 0; 
     634 
     635  for (int rank = 0; rank < mpiSize; rank++) 
     636  { 
     637    if (nbSendNode[rank] > 0) 
     638    { 
     639      MPI_Issend(sendBuffer2[rank], sendMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]); 
     640      nbSendRequest++; 
     641    } 
     642    if (nbRecvNode[rank] > 0) 
     643    { 
     644      MPI_Irecv(recvBuffer2[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
     645      nbRecvRequest++; 
     646    } 
     647  } 
     648 
     649  MPI_Waitall(nbRecvRequest, recvRequest, status); 
     650  MPI_Waitall(nbSendRequest, sendRequest, status); 
     651 
     652  int nbNeighbourNodes = 0; 
     653  for (int rank = 0; rank < mpiSize; rank++) 
     654    nbNeighbourNodes += nbRecvNode[rank]; 
     655 
     656  neighbourElements = new Elt[nbNeighbourNodes]; 
     657  nbNeighbourElements = nbNeighbourNodes; 
     658 
     659  int index = 0; 
     660  for (int rank = 0; rank < mpiSize; rank++) 
     661  { 
     662    pos[rank] = 0; 
     663    for (int j = 0; j < nbRecvNode[rank]; j++) 
     664    { 
     665      unpackPolygon(neighbourElements[index], recvBuffer2[rank], pos[rank]); 
     666      neighbourElements[index].id.ind = sstree.localTree.leafs.size() + index; 
     667      index++; 
     668    } 
     669  } 
     670  for (int rank = 0; rank < mpiSize; rank++) 
     671  { 
     672    if (nbRecvNode[rank] > 0) delete [] recvBuffer2[rank]; 
     673    if (nbSendNode[rank] > 0) delete [] sendBuffer2[rank]; 
     674  } 
     675  delete [] recvBuffer2; 
     676  delete [] sendBuffer2; 
     677  delete [] sendMessageSize; 
     678  delete [] recvMessageSize; 
     679  delete [] nbSendNode; 
     680  delete [] nbRecvNode; 
     681  delete [] sendRequest; 
     682  delete [] recvRequest; 
     683  delete [] status; 
     684  delete [] pos; 
     685 
     686  /* re-compute on received elements to avoid having to send this information */ 
     687  neighbourNodes.resize(nbNeighbourNodes); 
     688  setCirclesAndLinks(neighbourElements, neighbourNodes); 
     689  cptAllEltsGeom(neighbourElements, nbNeighbourNodes, srcGrid.pole); 
     690 
     691  /* the local SS tree must include nodes from other cpus if they are potential 
    673692           intersector of a local node */ 
    674         sstree.localTree.insertNodes(neighbourNodes); 
    675  
    676         /* for every local element, 
     693  sstree.localTree.insertNodes(neighbourNodes); 
     694 
     695  /* for every local element, 
    677696           use the SS-tree to find all elements (including neighbourElements) 
    678697           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 
     698     then check all canditates for common edges to build up connectivity information 
     699  */ 
     700  for (int j = 0; j < sstree.localTree.leafs.size(); j++) 
     701  { 
     702    Node& node = sstree.localTree.leafs[j]; 
     703 
     704    /* find all leafs whoes circles that intersect node's circle and save into node->intersectors */ 
     705    node.search(sstree.localTree.root); 
     706 
     707    Elt *elt = (Elt *)(node.data); 
     708 
     709    for (int i = 0; i < elt->n; i++) elt->neighbour[i] = NOT_FOUND; 
     710 
     711    /* for element `elt` loop through all nodes in the SS-tree 
    693712                   whoes circles intersect with the circle around `elt` (the SS intersectors) 
    694713                   and check if they are neighbours in the sense that the two elements share an edge. 
    695714                   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                 } 
     715    for (list<NodePtr>::iterator it = (node.intersectors).begin(); it != (node.intersectors).end(); ++it) 
     716    { 
     717      Elt *elt2 = (Elt *)((*it)->data); 
     718      set_neighbour(*elt, *elt2); 
     719    } 
    701720 
    702721/* 
    703                 for (int i = 0; i < elt->n; i++) 
    704                 { 
    705                         if (elt->neighbour[i] == NOT_FOUND) 
    706                                 error_exit("neighbour not found"); 
    707                 } 
     722    for (int i = 0; i < elt->n; i++) 
     723    { 
     724      if (elt->neighbour[i] == NOT_FOUND) 
     725        error_exit("neighbour not found"); 
     726    } 
    708727*/ 
    709         } 
     728  } 
    710729} 
    711730 
     
    713732void Mapper::computeIntersection(Elt *elements, int nbElements) 
    714733{ 
    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; 
     734  int mpiSize, mpiRank; 
     735  MPI_Comm_size(communicator, &mpiSize); 
     736  MPI_Comm_rank(communicator, &mpiRank); 
     737 
     738  MPI_Barrier(communicator); 
     739 
     740  vector<Node> *routingList = new vector<Node>[mpiSize]; 
     741 
     742  vector<Node> routeNodes;  routeNodes.reserve(nbElements); 
     743  for (int j = 0; j < nbElements; j++) 
     744  { 
     745    elements[j].id.ind = j; 
     746    elements[j].id.rank = mpiRank; 
     747    routeNodes.push_back(Node(elements[j].x, cptRadius(elements[j]), &elements[j])); 
     748  } 
     749 
     750  vector<vector<int> > routes(routeNodes.size()); 
     751  sstree.routeIntersections(routes, routeNodes); 
     752  for (int i = 0; i < routes.size(); ++i) 
     753    for (int k = 0; k < routes[i].size(); ++k) 
     754      routingList[routes[i][k]].push_back(routeNodes[i]); 
     755 
     756  if (verbose >= 2) 
     757  { 
     758    cout << " --> rank  " << mpiRank << " nbElements " << nbElements << " : "; 
     759    for (int rank = 0; rank < mpiSize; rank++) 
     760      cout << routingList[rank].size() << "   "; 
     761    cout << endl; 
     762  } 
     763  MPI_Barrier(communicator); 
     764 
     765  int *nbSendNode = new int[mpiSize]; 
     766  int *nbRecvNode = new int[mpiSize]; 
     767  int *sentMessageSize = new int[mpiSize]; 
     768  int *recvMessageSize = new int[mpiSize]; 
     769 
     770  for (int rank = 0; rank < mpiSize; rank++) 
     771  { 
     772    nbSendNode[rank] = routingList[rank].size(); 
     773    sentMessageSize[rank] = 0; 
     774    for (size_t j = 0; j < routingList[rank].size(); j++) 
     775    { 
     776      Elt *elt = (Elt *) (routingList[rank][j].data); 
     777      sentMessageSize[rank] += packedPolygonSize(*elt); 
     778    } 
     779  } 
     780 
     781  MPI_Alltoall(nbSendNode, 1, MPI_INT, nbRecvNode, 1, MPI_INT, communicator); 
     782  MPI_Alltoall(sentMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator); 
     783 
     784  int total = 0; 
     785 
     786  for (int rank = 0; rank < mpiSize; rank++) 
     787  { 
     788    total = total + nbRecvNode[rank]; 
     789  } 
     790 
     791  if (verbose >= 2) cout << "---> rank " << mpiRank << " : compute intersection : total received nodes  " << total << endl; 
     792  char **sendBuffer = new char*[mpiSize]; 
     793  char **recvBuffer = new char*[mpiSize]; 
     794  int *pos = new int[mpiSize]; 
     795 
     796  for (int rank = 0; rank < mpiSize; rank++) 
     797  { 
     798    if (nbSendNode[rank] > 0) sendBuffer[rank] = new char[sentMessageSize[rank]]; 
     799    if (nbRecvNode[rank] > 0) recvBuffer[rank] = new char[recvMessageSize[rank]]; 
     800  } 
     801 
     802  for (int rank = 0; rank < mpiSize; rank++) 
     803  { 
     804    pos[rank] = 0; 
     805    for (size_t j = 0; j < routingList[rank].size(); j++) 
     806    { 
     807      Elt* elt = (Elt *) (routingList[rank][j].data); 
     808      packPolygon(*elt, sendBuffer[rank], pos[rank]); 
     809    } 
     810  } 
     811  delete [] routingList; 
     812 
     813  int nbSendRequest = 0; 
     814  int nbRecvRequest = 0; 
     815  MPI_Request *sendRequest = new MPI_Request[mpiSize]; 
     816  MPI_Request *recvRequest = new MPI_Request[mpiSize]; 
     817  MPI_Status   *status = new MPI_Status[mpiSize]; 
     818 
     819  for (int rank = 0; rank < mpiSize; rank++) 
     820  { 
     821    if (nbSendNode[rank] > 0) 
     822    { 
     823      MPI_Issend(sendBuffer[rank], sentMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]); 
     824      nbSendRequest++; 
     825    } 
     826    if (nbRecvNode[rank] > 0) 
     827    { 
     828      MPI_Irecv(recvBuffer[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
     829      nbRecvRequest++; 
     830    } 
     831  } 
     832 
     833  MPI_Waitall(nbRecvRequest, recvRequest, status); 
     834  MPI_Waitall(nbSendRequest, sendRequest, status); 
     835  char **sendBuffer2 = new char*[mpiSize]; 
     836  char **recvBuffer2 = new char*[mpiSize]; 
     837 
     838  double tic = cputime(); 
     839  for (int rank = 0; rank < mpiSize; rank++) 
     840  { 
     841    sentMessageSize[rank] = 0; 
     842 
     843    if (nbRecvNode[rank] > 0) 
     844    { 
     845      Elt *recvElt = new Elt[nbRecvNode[rank]]; 
     846      pos[rank] = 0; 
     847      for (int j = 0; j < nbRecvNode[rank]; j++) 
     848      { 
     849        unpackPolygon(recvElt[j], recvBuffer[rank], pos[rank]); 
     850        cptEltGeom(recvElt[j], tgtGrid.pole); 
     851        Node recvNode(recvElt[j].x, cptRadius(recvElt[j]), &recvElt[j]); 
     852        recvNode.search(sstree.localTree.root); 
     853        /* for a node holding an element of the target, loop throught candidates for intersecting source */ 
     854        for (list<NodePtr>::iterator it = (recvNode.intersectors).begin(); it != (recvNode.intersectors).end(); ++it) 
     855        { 
     856          Elt *elt2 = (Elt *) ((*it)->data); 
     857          /* recvElt is target, elt2 is source */ 
     858//          intersect(&recvElt[j], elt2); 
     859          intersect_ym(&recvElt[j], elt2); 
     860        } 
     861 
     862        if (recvElt[j].is.size() > 0) sentMessageSize[rank] += packIntersectionSize(recvElt[j]); 
     863 
     864        // here recvNode goes out of scope 
     865      } 
     866 
     867      if (sentMessageSize[rank] > 0) 
     868      { 
     869        sentMessageSize[rank] += sizeof(int); 
     870        sendBuffer2[rank] = new char[sentMessageSize[rank]]; 
     871        *((int *) sendBuffer2[rank]) = 0; 
     872        pos[rank] = sizeof(int); 
     873        for (int j = 0; j < nbRecvNode[rank]; j++) 
     874        { 
     875          packIntersection(recvElt[j], sendBuffer2[rank], pos[rank]); 
     876          //FIXME should be deleted: recvElt[j].delete_intersections(); // intersection areas have been packed to buffer and won't be used any more 
     877        } 
     878      } 
     879      delete [] recvElt; 
     880 
     881    } 
     882  } 
     883  delete [] pos; 
     884 
     885  for (int rank = 0; rank < mpiSize; rank++) 
     886  { 
     887    if (nbSendNode[rank] > 0) delete [] sendBuffer[rank]; 
     888    if (nbRecvNode[rank] > 0) delete [] recvBuffer[rank]; 
     889    nbSendNode[rank] = 0; 
     890  } 
     891 
     892  if (verbose >= 2) cout << "Rank " << mpiRank << "  Compute (internal) intersection " << cputime() - tic << " s" << endl; 
     893  MPI_Alltoall(sentMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator); 
     894 
     895  for (int rank = 0; rank < mpiSize; rank++) 
     896    if (recvMessageSize[rank] > 0) 
     897      recvBuffer2[rank] = new char[recvMessageSize[rank]]; 
     898 
     899  nbSendRequest = 0; 
     900  nbRecvRequest = 0; 
     901 
     902  for (int rank = 0; rank < mpiSize; rank++) 
     903  { 
     904    if (sentMessageSize[rank] > 0) 
     905    { 
     906      MPI_Issend(sendBuffer2[rank], sentMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]); 
     907      nbSendRequest++; 
     908    } 
     909    if (recvMessageSize[rank] > 0) 
     910    { 
     911      MPI_Irecv(recvBuffer2[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
     912      nbRecvRequest++; 
     913    } 
     914  } 
     915 
     916  MPI_Waitall(nbRecvRequest, recvRequest, status); 
     917  MPI_Waitall(nbSendRequest, sendRequest, status); 
     918 
     919  delete [] sendRequest; 
     920  delete [] recvRequest; 
     921  delete [] status; 
     922  for (int rank = 0; rank < mpiSize; rank++) 
     923  { 
     924    if (nbRecvNode[rank] > 0) 
     925    { 
     926      if (sentMessageSize[rank] > 0) 
     927        delete [] sendBuffer2[rank]; 
     928    } 
     929 
     930    if (recvMessageSize[rank] > 0) 
     931    { 
     932      unpackIntersection(elements, recvBuffer2[rank]); 
     933      delete [] recvBuffer2[rank]; 
     934    } 
     935  } 
     936  delete [] sendBuffer2; 
     937  delete [] recvBuffer2; 
     938  delete [] sendBuffer; 
     939  delete [] recvBuffer; 
     940 
     941  delete [] nbSendNode; 
     942  delete [] nbRecvNode; 
     943  delete [] sentMessageSize; 
     944  delete [] recvMessageSize; 
    926945} 
    927946 
    928947Mapper::~Mapper() 
    929948{ 
    930         delete [] remapMatrix; 
    931         delete [] srcAddress; 
    932         delete [] srcRank; 
    933         delete [] dstAddress; 
    934         if (neighbourElements) delete [] neighbourElements; 
    935 } 
    936  
    937 } 
     949  delete [] remapMatrix; 
     950  delete [] srcAddress; 
     951  delete [] srcRank; 
     952  delete [] dstAddress; 
     953  if (neighbourElements) delete [] neighbourElements; 
     954} 
     955 
     956} 
  • XIOS/dev/dev_trunk_omp/extern/remap/src/mapper.hpp

    r1602 r1619  
    2323       void setVerbosity(verbosity v) {verbose=v ;} 
    2424 
    25        void setSourceMesh(const double* boundsLon, const double* boundsLat, int nVertex, int nbCells, const double* pole, const long int* globalId=NULL) ; 
    26        void setTargetMesh(const double* boundsLon, const double* boundsLat, int nVertex, int nbCells, const double* pole, const long int* globalId=NULL) ; 
     25       void setSourceMesh(const double* boundsLon, const double* boundsLat, const double* area, int nVertex, int nbCells, const double* pole, const long int* globalId=NULL) ; 
     26       void setTargetMesh(const double* boundsLon, const double* boundsLat, const double* area, int nVertex, int nbCells, const double* pole, const long int* globalId=NULL) ; 
    2727       void setSourceValue(const double* val) ; 
    2828       void getTargetValue(double* val) ; 
  • XIOS/dev/dev_trunk_omp/extern/remap/src/meshutil.cpp

    r1581 r1619  
    22#include "elt.hpp" 
    33#include "polyg.hpp" 
     4#include "intersection_ym.hpp" 
     5#include "earcut.hpp" 
     6#include <vector> 
    47 
    58namespace sphereRemap { 
    69 
    710using namespace std; 
     11 
     12double computePolygoneArea(Elt& a, const Coord &pole) 
     13{ 
     14  using N = uint32_t; 
     15  using Point = array<double, 2>; 
     16  vector<Point> vect_points; 
     17  vector< vector<Point> > polyline; 
     18   
     19  vector<Coord> dstPolygon ; 
     20  createGreatCirclePolygon(a, pole, dstPolygon) ; 
     21 
     22  int na=dstPolygon.size() ; 
     23  Coord *a_gno   = new Coord[na]; 
     24   
     25  Coord OC=barycentre(a.vertex,a.n) ; 
     26  Coord Oz=OC ; 
     27  Coord Ox=crossprod(Coord(0,0,1),Oz) ; 
     28// choose Ox not too small to avoid rounding error 
     29  if (norm(Ox)< 0.1) Ox=crossprod(Coord(0,1,0),Oz) ; 
     30  Ox=Ox*(1./norm(Ox)) ; 
     31  Coord Oy=crossprod(Oz,Ox) ; 
     32  double cos_alpha; 
     33 
     34  for(int n=0; n<na;n++) 
     35  { 
     36    cos_alpha=scalarprod(OC,dstPolygon[n]) ; 
     37    a_gno[n].x=scalarprod(dstPolygon[n],Ox)/cos_alpha ; 
     38    a_gno[n].y=scalarprod(dstPolygon[n],Oy)/cos_alpha ; 
     39    a_gno[n].z=scalarprod(dstPolygon[n],Oz)/cos_alpha ; // must be equal to 1 
     40 
     41    vect_points.push_back( array<double, 2>() ); 
     42    vect_points[n][0] = a_gno[n].x; 
     43    vect_points[n][1] = a_gno[n].y; 
     44 
     45  } 
     46 
     47  polyline.push_back(vect_points); 
     48  vector<N> indices_a_gno = mapbox::earcut<N>(polyline); 
     49   
     50  double area_a_gno=0 ; 
     51  for(int i=0;i<indices_a_gno.size()/3;++i) 
     52    { 
     53      Coord x0 = Ox * polyline[0][indices_a_gno[3*i]][0] + Oy* polyline[0][indices_a_gno[3*i]][1] + Oz ; 
     54      Coord x1 = Ox * polyline[0][indices_a_gno[3*i+1]][0] + Oy* polyline[0][indices_a_gno[3*i+1]][1] + Oz ; 
     55      Coord x2 = Ox * polyline[0][indices_a_gno[3*i+2]][0] + Oy* polyline[0][indices_a_gno[3*i+2]][1] + Oz ; 
     56      area_a_gno+=triarea(x0 * (1./norm(x0)),x1* (1./norm(x1)), x2* (1./norm(x2))) ; 
     57    } 
     58 
     59  vect_points.clear(); 
     60  polyline.clear(); 
     61  indices_a_gno.clear(); 
     62  return area_a_gno ; 
     63} 
     64 
    865 
    966void cptEltGeom(Elt& elt, const Coord &pole) 
     
    1471  elt.area = airbar(elt.n, elt.vertex, elt.edge, elt.d, pole, gg); 
    1572  elt.x = gg; 
    16 } 
     73// overwrite area computation  
     74 
     75  elt.area =  computePolygoneArea(elt, pole) ; 
     76} 
     77 
    1778 
    1879void cptAllEltsGeom(Elt *elt, int N, const Coord &pole) 
  • XIOS/dev/dev_trunk_omp/extern/remap/src/polyg.cpp

    r1579 r1619  
    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 
  • XIOS/dev/dev_trunk_omp/src/config/domain_attribute.conf

    r1553 r1619  
    5858 
    5959DECLARE_ARRAY(double, 2, area) 
     60DECLARE_ATTRIBUTE(double, radius) 
    6061 
    6162DECLARE_ENUM4(type,rectilinear,curvilinear,unstructured, gaussian) 
  • XIOS/dev/dev_trunk_omp/src/config/interpolate_domain_attribute.conf

    r1269 r1619  
    1010DECLARE_ATTRIBUTE(bool, write_weight) 
    1111DECLARE_ENUM2(read_write_convention, c, fortran) 
     12DECLARE_ATTRIBUTE(bool, use_area) 
  • XIOS/dev/dev_trunk_omp/src/context_client.cpp

    r1601 r1619  
    9696    { 
    9797      list<int> ranks = event.getRanks(); 
    98  
     98      info(100)<<"Event "<<timeLine<<" of context "<<context->getId()<<endl ; 
    9999      if (CXios::checkEventSync) 
    100100      { 
     
    124124        { 
    125125          event.send(timeLine, sizes, buffList); 
     126          info(100)<<"Event "<<timeLine<<" of context "<<context->getId()<<"  sent"<<endl ; 
    126127 
    127128          checkBuffers(ranks); 
     
    142143          info(100)<<"DEBUG : temporaly event created : timeline "<<timeLine<<endl ; 
    143144          event.send(timeLine, tmpBufferedEvent.sizes, tmpBufferedEvent.buffers); 
     145          info(100)<<"Event "<<timeLine<<" of context "<<context->getId()<<"  sent"<<endl ; 
    144146        } 
    145147      } 
  • XIOS/dev/dev_trunk_omp/src/io/onetcdf4_impl.hpp

    r1442 r1619  
    7373    } 
    7474    char *PtrArrayStr ; 
    75     PtrArrayStr=new char[stringArrayLen] ; 
     75    PtrArrayStr=new char[stringArrayLen*data.numElements()] ; 
     76    memset (PtrArrayStr,' ',stringArrayLen*data.numElements()); 
     77    size_t offset=0 ; 
    7678    Array<StdString,1>::const_iterator it, itb=data.begin(), ite=data.end() ; 
    77     int lineNb = 0; 
    78     for(it=itb;it!=ite;++it) 
     79    for(it=itb;it!=ite;++it, offset+=stringArrayLen) 
    7980    { 
    80       it->copy(PtrArrayStr,it->size()) ; 
    81       PtrArrayStr[it->size()]='\0' ; 
    82       sstart[0] = lineNb; 
    83       sstart[dimArrayLen] = 0; 
    84       scount[0] = 1; 
    85       scount[dimArrayLen] = it->size() + 1; 
    86       CTimer::get("CONetCDF4::writeData writeData_").resume(); 
    87       this->writeData_(grpid, varid, sstart, scount, PtrArrayStr); 
    88       CTimer::get("CONetCDF4::writeData writeData_").suspend(); 
    89       ++lineNb; 
     81      it->copy(PtrArrayStr+offset,it->size()) ; 
     82      PtrArrayStr[offset+it->size()]='\0' ; 
    9083    } 
     84 
     85     CTimer::get("CONetCDF4::writeData writeData_").resume(); 
     86     this->writeData_(grpid, varid, sstart, scount, PtrArrayStr); 
     87     CTimer::get("CONetCDF4::writeData writeData_").suspend(); 
     88 
    9189    delete []  PtrArrayStr; 
    9290  } 
  • XIOS/dev/dev_trunk_omp/src/node/interpolate_domain.cpp

    r1269 r1619  
    6565    if (this->read_write_convention.isEmpty()) this->read_write_convention.setValue(read_write_convention_attr::fortran); 
    6666 
     67 
    6768  } 
    6869 
  • XIOS/dev/dev_trunk_omp/src/transformation/domain_algorithm_interpolate.cpp

    r1601 r1619  
    305305  CArray<double,2> boundsLonSrcUnmasked(nVertexSrc,nSrcLocalUnmasked); 
    306306  CArray<double,2> boundsLatSrcUnmasked(nVertexSrc,nSrcLocalUnmasked); 
     307  CArray<double,1> areaSrcUnmasked(nSrcLocalUnmasked); 
     308   
    307309  long int * globalSrcUnmasked = new long int [nSrcLocalUnmasked]; 
    308310 
    309311  nSrcLocalUnmasked=0 ; 
     312  bool hasSrcArea=domainSrc_->hasArea && !domainSrc_->radius.isEmpty() && !interpDomain_->use_area.isEmpty() && interpDomain_->use_area==true  ; 
     313  double srcAreaFactor ; 
     314  if (hasSrcArea) srcAreaFactor=1./(domainSrc_->radius*domainSrc_->radius) ; 
     315   
    310316  for (int idx=0 ; idx < nSrcLocal; idx++) 
    311317  { 
     
    317323        boundsLatSrcUnmasked(n,nSrcLocalUnmasked) = boundsLatSrc(n,idx) ; 
    318324      } 
     325      if (hasSrcArea) areaSrcUnmasked(nSrcLocalUnmasked) = domainSrc_->areavalue(idx)*srcAreaFactor ; 
    319326      globalSrcUnmasked[nSrcLocalUnmasked]=globalSrc[idx] ; 
    320327      ++nSrcLocalUnmasked ; 
    321328    } 
    322329  } 
    323  
     330  
    324331 
    325332  int nDstLocalUnmasked = 0 ; 
     
    328335  CArray<double,2> boundsLonDestUnmasked(nVertexDest,nDstLocalUnmasked); 
    329336  CArray<double,2> boundsLatDestUnmasked(nVertexDest,nDstLocalUnmasked); 
     337  CArray<double,1>   areaDstUnmasked(nDstLocalUnmasked); 
     338 
    330339  long int * globalDstUnmasked = new long int [nDstLocalUnmasked]; 
    331340 
    332341  nDstLocalUnmasked=0 ; 
     342  bool hasDstArea=domainDest_->hasArea && !domainDest_->radius.isEmpty() && !interpDomain_->use_area.isEmpty() && interpDomain_->use_area==true ; 
     343  double dstAreaFactor ; 
     344  if (hasDstArea) dstAreaFactor=1./(domainDest_->radius*domainDest_->radius) ; 
    333345  for (int idx=0 ; idx < nDstLocal; idx++) 
    334346  { 
     
    340352        boundsLatDestUnmasked(n,nDstLocalUnmasked) = boundsLatDest(n,idx) ; 
    341353      } 
     354      if (hasDstArea) areaDstUnmasked(nDstLocalUnmasked) = domainDest_->areavalue(idx)*dstAreaFactor ; 
    342355      globalDstUnmasked[nDstLocalUnmasked]=globalDst[idx] ; 
    343356      ++nDstLocalUnmasked ; 
     
    345358  } 
    346359 
    347   mapper.setSourceMesh(boundsLonSrcUnmasked.dataFirst(), boundsLatSrcUnmasked.dataFirst(), nVertexSrc, nSrcLocalUnmasked, &srcPole[0], globalSrcUnmasked); 
    348   mapper.setTargetMesh(boundsLonDestUnmasked.dataFirst(), boundsLatDestUnmasked.dataFirst(), nVertexDest, nDstLocalUnmasked, &dstPole[0], globalDstUnmasked); 
     360  double* ptAreaSrcUnmasked = NULL ; 
     361  if (hasSrcArea) ptAreaSrcUnmasked=areaSrcUnmasked.dataFirst() ; 
     362 
     363  double* ptAreaDstUnmasked = NULL ; 
     364  if (hasDstArea) ptAreaDstUnmasked=areaDstUnmasked.dataFirst() ; 
     365 
     366  mapper.setSourceMesh(boundsLonSrcUnmasked.dataFirst(), boundsLatSrcUnmasked.dataFirst(), ptAreaSrcUnmasked, nVertexSrc, nSrcLocalUnmasked, &srcPole[0], globalSrcUnmasked); 
     367  mapper.setTargetMesh(boundsLonDestUnmasked.dataFirst(), boundsLatDestUnmasked.dataFirst(), ptAreaDstUnmasked, nVertexDest, nDstLocalUnmasked, &dstPole[0], globalDstUnmasked); 
    349368 
    350369  std::vector<double> timings = mapper.computeWeights(orderInterp,renormalize,quantity); 
Note: See TracChangeset for help on using the changeset viewer.