Changeset 1646


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

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

Location:
XIOS/dev/dev_trunk_omp
Files:
6 added
247 edited

Legend:

Unmodified
Added
Removed

  • XIOS/dev/dev_trunk_omp/bld.cfg

    r1628 r1646  
    5656#bld::target test_client.exe 
    5757bld::target test_omp.exe 
    58 bld::target test_omp2.exe 
    59 bld::target test_send.exe 
    60 bld::target test_send2.exe 
     58#bld::target test_omp2.exe 
     59#bld::target test_send.exe 
     60#bld::target test_send2.exe 
    6161#bld::target test_unstruct_complete.exe 
    6262#bld::target test_unstructured.exe 
    … …  
    7171bld::tool::ld        %LINKER 
    7272bld::tool::ldflags   %LD_FLAGS  
    73 bld::tool::cflags    %CFLAGS %CBASE_INC -I${PWD}/extern/src_netcdf -I${PWD}/extern/boost/include -I${PWD}/extern/rapidxml/include -I${PWD}/extern/blitz/include  
     73bld::tool::cflags    %CFLAGS %CBASE_INC -I${PWD}/extern/src_netcdf -I${PWD}/extern/boost/include -I${PWD}/extern/rapidxml/include -I${PWD}/extern/blitz/include -I${PWD}/extern/src_ep_dev  
    7474bld::tool::fflags    %FFLAGS %FBASE_INC  
    7575bld::tool::cppkeys   %CPP_KEY 

  • XIOS/dev/dev_trunk_omp/extern/ep_dev/ep_message.cpp

    r1604 r1646  
    133133  int Message_Check(MPI_Comm comm) 
    134134  { 
    135     if(comm->is_ep) return Message_Check_endpoint(comm); 
     135    if(comm->is_ep) return Message_Check_endpoint(comm);   
    136136  } 
    137137   

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

    r1602 r1646  
    11#include "mpi.hpp" 
     2#ifdef _usingEP 
    23using namespace ep_lib; 
     4#endif 
    35 
    46namespace sphereRemap { 

  • XIOS/dev/dev_trunk_omp/extern/remap/src/elt.hpp

    r1630 r1646  
    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

    r1630 r1646  
    1515#include "cputime.hpp" // cputime 
    1616 
     17#ifdef _usingEP 
    1718using namespace ep_lib; 
     19#endif 
    1820 
    1921using namespace sphereRemap ; 
    … …  
    4345        assert(n_cell_dst >= 4); 
    4446        assert(1 <= order && order <= 2); 
    45  
     47  double* src_area=NULL ; 
     48  double* dst_area=NULL ; 
    4649  mapper = new Mapper(MPI_COMM_WORLD); 
    4750  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 ) ; 
     51  mapper->setSourceMesh(src_bounds_lon, src_bounds_lat, src_area, n_vert_per_cell_src, n_cell_src, src_pole ) ; 
     52  mapper->setTargetMesh(dst_bounds_lon, dst_bounds_lat, dst_area, n_vert_per_cell_dst, n_cell_dst, dst_pole ) ; 
    5053 
    5154/* 

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

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

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

    r1630 r1646  
    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

    r1630 r1646  
    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/mpi_cascade.cpp

    r1602 r1646  
    11#include "mpi_cascade.hpp" 
    22#include <iostream> 
     3#ifdef _usingEP 
    34using namespace ep_lib; 
     5#endif 
    46 
    57namespace sphereRemap { 

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

    r1602 r1646  
    55#include "timerRemap.hpp" 
    66#include <iostream> 
     7#ifdef _usingEP 
    78using namespace ep_lib; 
     9#endif 
    810 
    911namespace sphereRemap { 

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

    r1602 r1646  
    1212 
    1313#include "parallel_tree.hpp" 
     14#ifdef _usingEP 
    1415using namespace ep_lib; 
     16#endif 
    1517 
    1618namespace sphereRemap { 

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

    r1630 r1646  
    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/extern/remap/src/timerRemap.cpp

    r1602 r1646  
    44#include <map> 
    55#include <iostream> 
     6#ifdef _usingEP 
    67using namespace ep_lib; 
     8#endif 
    79 
    810namespace sphereRemap { 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_accumulate.cpp

    r1603 r1646  
     1#ifdef _usingEP 
    12#include "ep_lib.hpp" 
    23#include <mpi.h> 
    … …  
    102103 
    103104} 
     105#endif 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_allgather.cpp

    r1603 r1646  
     1#ifdef _usingEP 
    12/*! 
    23   \file ep_gather.cpp 
    … …  
    135136 
    136137} 
     138#endif 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_allgatherv.cpp

    r1603 r1646  
     1#ifdef _usingEP 
    12/*! 
    23   \file ep_gather.cpp 
    … …  
    121122 
    122123} 
     124#endif 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_allocate.cpp

    r1603 r1646  
     1#ifdef _usingEP 
    12#include "ep_lib.hpp" 
    23#include <mpi.h> 
    … …  
    8586 
    8687} 
     88#endif 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_allreduce.cpp

    r1603 r1646  
     1#ifdef _usingEP 
    12/*! 
    23   \file ep_reduce.cpp 
    … …  
    7273 
    7374} 
    74  
     75#endif 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_alltoall.cpp

    r1603 r1646  
     1#ifdef _usingEP 
    12#include "ep_lib.hpp" 
    23#include <mpi.h> 
    … …  
    6465 
    6566 
     67#endif 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_barrier.cpp

    r1603 r1646  
     1#ifdef _usingEP 
    12#include "ep_lib.hpp" 
    23#include <mpi.h> 
    … …  
    5152} 
    5253 
    53  
     54#endif 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_barrier.hpp

    r1603 r1646  
    22#define EP_BARRIER_HPP_INCLUDED 
    33 
     4#ifdef _usingEP 
    45 
    56namespace ep_lib 
    … …  
    5960 
    6061} 
    61  
     62#endif 
    6263 
    6364 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_bcast.cpp

    r1603 r1646  
     1#ifdef _usingEP 
    12/*! 
    23   \file ep_bcast.cpp 
    … …  
    7879 
    7980} 
     81#endif 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_buffer.hpp

    r1603 r1646  
    11#ifndef EP_BUFFER_HPP_INCLUDED 
    22#define EP_BUFFER_HPP_INCLUDED 
    3  
     3#ifdef _usingEP 
    44 
    55namespace ep_lib 
    … …  
    1515 
    1616 
    17  
     17#endif 
    1818#endif // EP_BUFFER_HPP_INCLUDED 
    1919 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_comm.hpp

    r1603 r1646  
    11#ifndef EP_COMM_HPP_INCLUDED 
    22#define EP_COMM_HPP_INCLUDED 
    3  
     3#ifdef _usingEP 
    44#include "ep_message.hpp" 
    55#include "ep_barrier.hpp" 
    … …  
    6363 
    6464 
    65  
     65#endif 
    6666#endif // EP_COMM_HPP_INCLUDED 
    6767 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_compare.cpp

    r1603 r1646  
     1#ifdef _usingEP 
    12#include "ep_lib.hpp" 
    23#include <mpi.h> 
    … …  
    2223 
    2324} 
     25#endif 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_create.cpp

    r1605 r1646  
     1#ifdef _usingEP 
    12/*! 
    23   \file ep_create.cpp 
    … …  
    135136 
    136137} //namespace ep_lib 
     138#endif 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_declaration.cpp

    r1603 r1646  
     1#ifdef _usingEP 
    12#include "ep_lib.hpp" 
    23 
    … …  
    110111 
    111112 
     113#endif 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_declaration.hpp

    r1603 r1646  
    11#ifndef EP_DECLARATION_HPP_INCLUDED 
    22#define EP_DECLARATION_HPP_INCLUDED 
     3 
     4#ifdef _usingEP 
    35 
    46#undef MPI_INT 
    … …  
    4850extern ep_lib::MPI_Info MPI_INFO_NULL; 
    4951 
     52#endif 
    5053 
    5154#endif // EP_DECLARATION_HPP_INCLUDED 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_dup.cpp

    r1603 r1646  
     1#ifdef _usingEP 
    12#include "ep_lib.hpp" 
    23#include <mpi.h> 
    … …  
    7273 
    7374} 
     75#endif 
    7476 
    75  

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_exscan.cpp

    r1603 r1646  
     1#ifdef _usingEP 
    12/*! 
    23   \file ep_scan.cpp 
    … …  
    363364 
    364365} 
     366#endif 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_fetch.cpp

    r1603 r1646  
     1#ifdef _usingEP 
    12#include "ep_lib.hpp" 
    23#include <mpi.h> 
    … …  
    2223 
    2324} 
     25#endif 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_finalize.cpp

    r1603 r1646  
     1#ifdef _usingEP 
    12#include "ep_lib.hpp" 
    23#include <mpi.h> 
    … …  
    2324} 
    2425 
    25  
     26#endif 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_fortran.cpp

    r1605 r1646  
     1#ifdef _usingEP 
    12#include "ep_lib.hpp" 
    23#include "ep_lib_fortran.hpp" 
    … …  
    8687 
    8788//} 
     89#endif 
    8890 
    89  

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_free.cpp

    r1603 r1646  
     1#ifdef _usingEP 
    12#include "ep_lib.hpp" 
    23#include <mpi.h> 
    … …  
    108109 
    109110 
     111#endif 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_gather.cpp

    r1603 r1646  
     1#ifdef _usingEP 
    12/*! 
    23   \file ep_gather.cpp 
    … …  
    126127  } 
    127128} 
     129#endif 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_gatherv.cpp

    r1603 r1646  
     1#ifdef _usingEP 
    12/*! 
    23   \file ep_gather.cpp 
    … …  
    9394    else                          MPI_Gather_local(&sendcount, 1, MPI_INT, local_recvcounts.data(), 0, comm); 
    9495 
    95  
     96     
     97    //if(is_root) printf("rank %d =================local_recvcounts = %d %d\n", ep_rank, local_recvcounts[0], local_recvcounts[1]); 
     98    //if(!is_root && ep_rank_loc==0 && mpi_rank != root_mpi_rank) printf("rank %d =================local_recvcounts = %d %d\n", ep_rank, local_recvcounts[0], local_recvcounts[1]); 
    9699 
    97100    if(is_master) 
    … …  
    108111    else                          MPI_Gatherv_local(sendbuf, sendcount, sendtype, local_recvbuf, local_recvcounts.data(), local_displs.data(), 0, comm); 
    109112 
     113     
     114    //if(is_root) printf("rank %d =================local_recvbuf = %d %d\n", ep_rank, static_cast<int*>(local_recvbuf)[0], static_cast<int*>(local_recvbuf)[1]); 
     115    //if(!is_root && ep_rank_loc==0 && mpi_rank != root_mpi_rank && ep_rank!=6) printf("rank %d =================local_recvbuf = %d %d\n", ep_rank, static_cast<int*>(local_recvbuf)[0], static_cast<int*>(local_recvbuf)[1]); 
     116    //if(!is_root && ep_rank_loc==0 && mpi_rank != root_mpi_rank && ep_rank==6) printf("rank %d =================local_recvbuf = %d %d %d\n", ep_rank, static_cast<int*>(local_recvbuf)[0], static_cast<int*>(local_recvbuf)[1], static_cast<int*>(local_recvbuf)[2]); 
    110117 
    111118    void* tmp_recvbuf; 
    112119    int tmp_recvbuf_size = std::accumulate(recvcounts, recvcounts+ep_size, 0); 
     120 
     121 
    113122 
    114123    if(is_root) tmp_recvbuf = new void*[datasize * tmp_recvbuf_size]; 
    … …  
    121130    if(is_master) 
    122131    { 
     132      int sendcount_mpi = 0; 
     133      for(int i=0; i<num_ep; i++) 
     134      { 
     135        sendcount_mpi += local_recvcounts[i]; 
     136      } 
     137 
    123138      for(int i=0; i<ep_size; i++) 
    124139      { 
    … …  
    130145 
    131146 
    132       ::MPI_Gatherv(local_recvbuf, sendcount*num_ep, to_mpi_type(sendtype), tmp_recvbuf, mpi_recvcounts.data(), mpi_displs.data(), to_mpi_type(recvtype), root_mpi_rank, to_mpi_comm(comm->mpi_comm)); 
     147      ::MPI_Gatherv(local_recvbuf, sendcount_mpi, to_mpi_type(sendtype), tmp_recvbuf, mpi_recvcounts.data(), mpi_displs.data(), to_mpi_type(recvtype), root_mpi_rank, to_mpi_comm(comm->mpi_comm)); 
     148      //printf("****************** rank %d, sendcount*num_ep = %d, sendcount_mpi = %d\n", ep_rank, sendcount*num_ep, sendcount_mpi); 
     149 
     150      /*if(is_root) printf("rank %d =================tmp_recvbuf = %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d\n", ep_rank, static_cast<int*>(tmp_recvbuf)[0], static_cast<int*>(tmp_recvbuf)[1], static_cast<int*>(tmp_recvbuf)[2], 
     151                                                                                                                                 static_cast<int*>(tmp_recvbuf)[3], static_cast<int*>(tmp_recvbuf)[4], static_cast<int*>(tmp_recvbuf)[5], 
     152                                                                                                                                 static_cast<int*>(tmp_recvbuf)[6], static_cast<int*>(tmp_recvbuf)[7], static_cast<int*>(tmp_recvbuf)[8], 
     153                                                                                                                                 static_cast<int*>(tmp_recvbuf)[9], static_cast<int*>(tmp_recvbuf)[10], static_cast<int*>(tmp_recvbuf)[11], 
     154                                                                                                                                 static_cast<int*>(tmp_recvbuf)[12], static_cast<int*>(tmp_recvbuf)[13], static_cast<int*>(tmp_recvbuf)[14], 
     155                                                                                                                                 static_cast<int*>(tmp_recvbuf)[15], static_cast<int*>(tmp_recvbuf)[16]);*/ 
    133156    }    
    134157 
    … …  
    174197 
    175198} 
     199#endif 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_get.cpp

    r1603 r1646  
     1#ifdef _usingEP 
    12#include "ep_lib.hpp" 
    23#include <mpi.h> 
    … …  
    5758 
    5859} 
     60#endif 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_global.hpp

    r1603 r1646  
    11#ifndef EP_GLOBAL_HPP_INCLUDED 
    22#define EP_GLOBAL_HPP_INCLUDED 
     3#ifdef _usingEP 
    34 
    45extern int MPI_MODE_NOPRECEDE; 
    … …  
    910     
    1011 
    11  
     12#endif 
    1213#endif // EP_GLOBAL_HPP_INCLUDED 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_info.hpp

    r1603 r1646  
    22#define EP_INFO_HPP_INCLUDED 
    33 
     4#ifdef _usingEP 
    45 
    56namespace ep_lib 
    … …  
    2122} 
    2223 
    23  
     24#endif 
    2425 
    2526#endif // EP_INFO_HPP_INCLUDED 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_init.cpp

    r1603 r1646  
     1#ifdef _usingEP 
    12#include "ep_lib.hpp" 
    23#include <mpi.h> 
    … …  
    4344} 
    4445 
    45  
     46#endif 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_intercomm.cpp

    r1603 r1646  
     1#ifdef _usingEP 
    12#include "ep_lib.hpp" 
    23#include <mpi.h> 
    … …  
    354355      if(is_real_involved) 
    355356      { 
     357        #pragma omp critical (_mpi_call) 
    356358        ::MPI_Intercomm_create(extracted_comm, local_leader_rank_in_extracted_comm, to_mpi_comm(peer_comm->mpi_comm), remote_leader_rank_in_peer_mpi, tag, &mpi_inter_comm); 
    357359        ::MPI_Intercomm_merge(mpi_inter_comm, !priority, intracomm); 
    … …  
    550552 
    551553} 
     554#endif 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_lib.cpp

    r1603 r1646  
     1#ifdef _usingEP 
    12#include "ep_lib.hpp" 
    23#include <mpi.h> 
    … …  
    218219} 
    219220 
    220  
    221  
     221#endif 
     222 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_lib.hpp

    r1603 r1646  
    11#ifndef EP_LIB_HPP_INCLUDED 
    22#define EP_LIB_HPP_INCLUDED 
     3#ifdef _usingEP 
    34 
    45#include "ep_type.hpp" 
    … …  
    8788 
    8889} 
    89  
     90#endif 
    9091#endif // EP_LIB_HPP_INCLUDED 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_lib_collective.hpp

    r1603 r1646  
    11#ifndef EP_LIB_COLLECTIVE_HPP_INCLUDED 
    22#define EP_LIB_COLLECTIVE_HPP_INCLUDED 
     3#ifdef _usingEP 
    34 
    45namespace ep_lib 
    … …  
    6364 
    6465} 
    65  
     66#endif 
    6667#endif // EP_LIB_COLLECTIVE_HPP_INCLUDED 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_lib_endpoint.hpp

    r1603 r1646  
    11#ifndef EP_LIB_ENDPOINT_HPP_INCLUDED 
    22#define EP_LIB_ENDPOINT_HPP_INCLUDED 
     3#ifdef _usingEP 
    34 
    45namespace ep_lib 
    … …  
    2324 
    2425} 
    25  
     26#endif 
    2627#endif // EP_LIB_ENDPOINT_HPP_INCLUDED 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_lib_fortran.hpp

    r1605 r1646  
    11#ifndef EP_LIB_FORTRAN_HPP_INCLUDED 
    22#define EP_LIB_FORTRAN_HPP_INCLUDED 
     3#ifdef _usingEP 
    34 
    45#include "ep_type.hpp" 
    … …  
    1011  ep_lib::MPI_Comm EP_Comm_f2c(void* comm); 
    1112//} 
    12  
     13#endif 
    1314 
    1415#endif // EP_LIB_FORTRAN_HPP_INCLUDED 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_lib_intercomm.hpp

    r1603 r1646  
    11#ifndef EP_LIB_INTERCOMM_HPP_INCLUDED 
    22#define EP_LIB_INTERCOMM_HPP_INCLUDED 
     3#ifdef _usingEP 
    34 
    45namespace ep_lib 
    … …  
    3435} 
    3536 
    36  
     37#endif 
    3738#endif // EP_LIB_INTERCOMM_HPP_INCLUDED 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_lib_local.hpp

    r1603 r1646  
    11#ifndef EP_LIB_LOCAL_HPP_INCLUDED 
    22#define EP_LIB_LOCAL_HPP_INCLUDED 
     3#ifdef _usingEP 
    34 
    45namespace ep_lib 
    … …  
    3233 
    3334} 
    34  
     35#endif 
    3536#endif // EP_LIB_LOCAL_HPP_INCLUDED 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_lib_mpi.hpp

    r1603 r1646  
    11#ifndef EP_LIB_MPI_HPP_INCLUDED 
    22#define EP_LIB_MPI_HPP_INCLUDED 
     3#ifdef _usingEP 
    34 
    45namespace ep_lib 
    … …  
    3334 
    3435} 
     36#endif 
    3537 
    3638#endif // EP_LIB_MPI_HPP_INCLUDED 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_lib_win.hpp

    r1603 r1646  
    11#ifndef EP_LIB_WIN_HPP_INCLUDED 
    22#define EP_LIB_WIN_HPP_INCLUDED 
     3#ifdef _usingEP 
    34 
    45namespace ep_lib 
    … …  
    5152 
    5253} 
    53  
     54#endif 
    5455#endif // EP_LIB_COLLECTIVE_HPP_INCLUDED 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_merge.cpp

    r1603 r1646  
     1#ifdef _usingEP 
    12#include "ep_lib.hpp" 
    23#include <mpi.h> 
    … …  
    104105   
    105106} 
     107#endif 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_message.cpp

    r1628 r1646  
     1#ifdef _usingEP 
    12/*! 
    23   \file ep_message.cpp 
    … …  
    237238 
    238239} 
    239  
     240#endif 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_message.hpp

    r1603 r1646  
     1#ifdef _usingEP 
    12#ifndef EP_MESSAGE_HPP_INCLUDED 
    23#define EP_MESSAGE_HPP_INCLUDED 
    … …  
    2526 
    2627#endif // EP_MESSAGE_HPP_INCLUDED 
    27  
     28#endif 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_mpi.hpp

    r1603 r1646  
    22#define EP_MPI_HPP_INCLUDED 
    33 
     4#ifdef _usingEP 
    45#include "ep_type.hpp" 
    56 
    … …  
    1516MPI_Request* to_mpi_request_ptr(ep_lib::MPI_Request request); 
    1617MPI_Message* to_mpi_message_ptr(ep_lib::MPI_Message message); 
     18#endif 
    1719 
    1820#endif // EP_MPI_HPP_INCLUDED 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_probe.cpp

    r1628 r1646  
     1#ifdef _usingEP 
    12#include "ep_lib.hpp" 
    23#include <mpi.h> 
    … …  
    265266 
    266267 
     268#endif 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_put.cpp

    r1603 r1646  
     1#ifdef _usingEP 
    12#include "ep_lib.hpp" 
    23#include <mpi.h> 
    … …  
    5960 
    6061} 
     62#endif 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_rank.cpp

    r1603 r1646  
     1#ifdef _usingEP 
    12#include "ep_lib.hpp" 
    23#include <mpi.h> 
    … …  
    2930 
    3031 
     32#endif 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_recv.cpp

    r1603 r1646  
     1#ifdef _usingEP 
    12/*! 
    23   \file ep_recv.cpp 
    … …  
    176177 
    177178 
     179#endif 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_reduce.cpp

    r1603 r1646  
     1#ifdef _usingEP 
    12/*! 
    23   \file ep_reduce.cpp 
    … …  
    341342} 
    342343 
     344#endif 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_reduce_scatter.cpp

    r1603 r1646  
     1#ifdef _usingEP 
    12/*! 
    23   \file ep_reduce.cpp 
    … …  
    8889} 
    8990 
     91#endif 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_request.hpp

    r1603 r1646  
    11#ifndef EP_REQUEST_HPP_INCLUDED 
    22#define EP_REQUEST_HPP_INCLUDED 
     3#ifdef _usingEP 
    34 
    45#include "ep_comm.hpp" 
    … …  
    4243} 
    4344 
    44  
     45#endif 
    4546 
    4647#endif // EP_REQUEST_HPP_INCLUDED 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_scan.cpp

    r1603 r1646  
     1#ifdef _usingEP 
    12/*! 
    23   \file ep_scan.cpp 
    … …  
    528529 
    529530} 
     531#endif 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_scatter.cpp

    r1603 r1646  
     1#ifdef _usingEP 
    12/*! 
    23   \file ep_gather.cpp 
    … …  
    132133 
    133134} 
     135#endif 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_scatterv.cpp

    r1603 r1646  
     1#ifdef _usingEP 
    12/*! 
    23   \file ep_gather.cpp 
    … …  
    153154 
    154155} 
     156#endif 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_send.cpp

    r1603 r1646  
     1#ifdef _usingEP 
    12/*! 
    23   \file ep_send.hpp 
    … …  
    202203 
    203204} 
    204  
    205  
    206  
    207  
    208  
    209  
     205#endif 
     206 
     207 
     208 
     209 
     210 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_size.cpp

    r1603 r1646  
     1#ifdef _usingEP 
    12#include "ep_lib.hpp" 
    23#include <mpi.h> 
    … …  
    4243} 
    4344 
    44  
     45#endif 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_split.cpp

    r1605 r1646  
     1#ifdef _usingEP 
    12#include "ep_lib.hpp" 
    23#include <mpi.h> 
    … …  
    321322 
    322323} 
     324#endif 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_status.hpp

    r1603 r1646  
    22#define EP_STATUS_HPP_INCLUDED 
    33 
     4#ifdef _usingEP 
    45 
    56namespace ep_lib 
    … …  
    2425} 
    2526 
    26  
     27#endif 
    2728 
    2829#endif // EP_STATUS_HPP_INCLUDED 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_test.cpp

    r1603 r1646  
     1#ifdef _usingEP 
    12/*! 
    23   \file ep_test.cpp 
    … …  
    113114} 
    114115 
     116#endif 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_type.cpp

    r1603 r1646  
     1#ifdef _usingEP 
    12#include "ep_lib.hpp" 
    23#include <mpi.h> 
    … …  
    3839 
    3940 
     41#endif 
    4042 
    4143 
    4244 
    4345 
    44  

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_type.hpp

    r1603 r1646  
    22#define EP_TYPE_HPP_INCLUDED 
    33 
     4#ifdef _usingEP 
    45 
    56#include <iostream> 
    … …  
    9192 
    9293 
    93  
     94#endif 
    9495#endif // EP_TYPE_HPP_INCLUDED 
    9596 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_wait.cpp

    r1603 r1646  
     1#ifdef _usingEP 
    12/*! 
    23  \file ep_wait.cpp 
    … …  
    111112} 
    112113 
     114#endif 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_win.cpp

    r1603 r1646  
     1#ifdef _usingEP 
    12#include "ep_lib.hpp" 
    23#include <mpi.h> 
    … …  
    180181 
    181182 
     183#endif 

  • XIOS/dev/dev_trunk_omp/extern/src_ep_dev/ep_window.hpp

    r1603 r1646  
    22#define EP_WINDOW_HPP_INCLUDED 
    33 
     4#ifdef _usingEP 
    45 
    56namespace ep_lib 
    … …  
    2122} 
    2223 
    23  
     24#endif 
    2425 
    2526#endif // EP_WINDOW_HPP_INCLUDED 

  • XIOS/dev/dev_trunk_omp/inputs/COMPLETE/iodef.xml

    r1129 r1646  
    1919 
    2020        <variable_group id="parameters" > 
    21           <variable id="info_level" type="int">100</variable> 
     21          <variable id="info_level" type="int">50</variable> 
    2222          <variable id="print_file" type="bool">true</variable> 
    2323        </variable_group> 

  • XIOS/dev/dev_trunk_omp/inputs/REMAP/iodef.xml

    r1610 r1646  
    183183        <variable_group id="parameters" > 
    184184          <variable id="using_server" type="bool">true</variable> 
    185           <variable id="info_level" type="int">50</variable> 
     185          <variable id="info_level" type="int">200</variable> 
    186186          <variable id="print_file" type="bool">true</variable> 
    187187        </variable_group> 

  • XIOS/dev/dev_trunk_omp/inputs/iodef.xml

    r1628 r1646  
    2121        <field field_ref="field_A_zoom" name="field_B" /> 
    2222     </file> 
    23      <file id="output1" name="output1" enabled=".FALSE."> 
     23     <file id="output1" name="output1" enabled=".TRUE."> 
    2424        <!-- <field field_ref="field_Domain" name="field_A" /> --> 
    2525        <field field_ref="field_A" name="field_A" />         
    … …  
    7979        <variable_group id="parameters" > 
    8080          <variable id="using_server" type="bool">false</variable> 
    81           <variable id="info_level" type="int">50</variable> 
     81          <variable id="info_level" type="int">80</variable> 
    8282          <variable id="print_file" type="bool">true</variable> 
    8383        </variable_group> 

  • XIOS/dev/dev_trunk_omp/src/array_new.hpp

    r1290 r1646  
    532532      virtual void fromString(const string& str) { istringstream iss(str); iss >> *this; initialized = true; } 
    533533      virtual string toString(void) const { ostringstream oss; oss << *this; return oss.str(); } 
     534 
     535      virtual string dump(void) const 
     536      { 
     537        ostringstream oss; 
     538        oss << this->shape()<<" "; 
     539        if (this->shape().numElements() == 1 && this->shape().dataFirst()[0] == 1) 
     540          oss << this->dataFirst()[0]; 
     541        else 
     542          oss << this->dataFirst()[0] <<" ... "<< this->dataFirst()[this->numElements()-1]; 
     543        return oss.str(); 
     544      } 
     545 
    534546      virtual void reset(void) { this->free(); initialized = false; } 
    535547      virtual bool isEmpty(void) const { return !initialized; } 

  • XIOS/dev/dev_trunk_omp/src/attribute.hpp

    r1158 r1646  
    4141            virtual StdString toString(void) const = 0; 
    4242            virtual void fromString(const StdString & str) = 0; 
     43            virtual StdString dump(void) const = 0; 
    4344            virtual bool isEqual(const CAttribute& ) = 0; 
    4445 

  • XIOS/dev/dev_trunk_omp/src/attribute_array.hpp

    r1219 r1646  
    5555            virtual bool toBuffer  (CBufferOut& buffer) const { return _toBuffer(buffer);} 
    5656            virtual bool fromBuffer(CBufferIn& buffer) { return _fromBuffer(buffer); } 
     57            virtual string dump(void) const { return _dump();} 
    5758 
    5859            virtual void generateCInterface(ostream& oss,const string& className) ; 
    … …  
    6970          CArray<T_numtype, N_rank> inheritedValue ; 
    7071          StdString _toString(void) const; 
     72          StdString _dump(void) const; 
    7173          void _fromString(const StdString & str); 
    7274          bool _toBuffer  (CBufferOut& buffer) const; 

  • XIOS/dev/dev_trunk_omp/src/attribute_array_impl.hpp

    r1219 r1646  
    129129    } 
    130130 
     131    template <typename T_numtype, int N_rank> 
     132    StdString CAttributeArray<T_numtype,N_rank>::_dump(void) const 
     133    { 
     134      StdOStringStream oss; 
     135      if (! isEmpty() && this->hasId() && (this->numElements()!=0)) 
     136        oss << this->getName() << "=\"" << CArray<T_numtype, N_rank>::dump() << "\""; 
     137      return (oss.str()); 
     138    } 
     139 
     140 
    131141      template <typename T_numtype, int N_rank> 
    132142         void CAttributeArray<T_numtype, N_rank>::_fromString(const StdString & str) 

  • XIOS/dev/dev_trunk_omp/src/attribute_enum.hpp

    r1219 r1646  
    6262            virtual StdString toString(void) const { return _toString();} 
    6363            virtual void fromString(const StdString & str) { if (str==resetInheritanceStr) { reset(); _canInherite=false ;}  else _fromString(str);} 
     64            virtual StdString dump(void) const { return _toString();} 
    6465 
    6566            virtual bool toBuffer  (CBufferOut& buffer) const { return _toBuffer(buffer);}  

  • XIOS/dev/dev_trunk_omp/src/attribute_map.cpp

    r1158 r1646  
    2828            att.second->reset(); 
    2929         } 
     30      } 
     31 
     32      ///-------------------------------------------------------------- 
     33      /*! 
     34         Dump of all non-empty attributes of an object 
     35      */ 
     36      StdString CAttributeMap::dumpXiosAttributes(void) const 
     37      { 
     38        int maxNbChar = 250; 
     39        StdString str; 
     40        typedef std::pair<StdString, CAttribute*> StdStrAttPair; 
     41        auto it = SuperClassMap::begin(), end = SuperClassMap::end(); 
     42        for (; it != end; it++) 
     43        { 
     44          const StdStrAttPair& att = *it; 
     45          if (!att.second->isEmpty()) 
     46          { 
     47            if (str.length() < maxNbChar) 
     48            { 
     49              str.append(att.second->dump()); 
     50              str.append(" "); 
     51            } 
     52            else if (str.length() == maxNbChar) 
     53            { 
     54              str.append("..."); 
     55            } 
     56          } 
     57        } 
     58        return str; 
    3059      } 
    3160 

  • XIOS/dev/dev_trunk_omp/src/attribute_map.hpp

    r1601 r1646  
    3838            void duplicateAttributes(const CAttributeMap* const _parent); 
    3939            void clearAllAttributes(void); 
     40            StdString dumpXiosAttributes(void) const; 
    4041 
    4142            void clearAttribute(const StdString& key); 

  • XIOS/dev/dev_trunk_omp/src/attribute_template.hpp

    r1601 r1646  
    7373//            virtual void toBinary  (StdOStream & os) const; 
    7474//            virtual void fromBinary(StdIStream & is); 
     75            virtual StdString dump(void) const { return _dump();} 
    7576 
    7677            virtual bool toBuffer  (CBufferOut& buffer) const { return _toBuffer(buffer);} 
    … …  
    9798          bool isEqual_(const CAttributeTemplate& attr); 
    9899          StdString _toString(void) const; 
     100          StdString _dump(void) const; 
    99101          void _fromString(const StdString & str); 
    100102          bool _toBuffer  (CBufferOut& buffer) const; 

  • XIOS/dev/dev_trunk_omp/src/attribute_template_impl.hpp

    r1601 r1646  
    199199        CType<T>::fromString(str) ; 
    200200      } 
     201 
     202      //--------------------------------------------------------------- 
     203 
     204      template <class T> 
     205         StdString CAttributeTemplate<T>::_dump(void) const 
     206      { 
     207         StdOStringStream oss; 
     208         if (!CType<T>::isEmpty() && this->hasId()) 
     209            oss << this->getName() << "=\"" << CType<T>::dump() << "\""; 
     210         return (oss.str()); 
     211      } 
     212 
    201213 
    202214      //--------------------------------------------------------------- 

  • XIOS/dev/dev_trunk_omp/src/buffer_client.cpp

    r1601 r1646  
    99 
    1010 
     11#ifdef _usingEP 
    1112using namespace ep_lib; 
     13#endif 
    1214 
    1315namespace xios 
    … …  
    3133    retBuffer = new CBufferOut(buffer[current], bufferSize); 
    3234    #pragma omp critical (_output) 
    33     info(10) << "CClientBuffer: allocated 2 x " << bufferSize << " bytes for server " << serverRank << " with a maximum of " << maxBufferedEvents << " buffered events" << endl; 
     35    { 
     36      info(10) << "CClientBuffer: allocated 2 x " << bufferSize << " bytes for server " << serverRank << " with a maximum of " << maxBufferedEvents << " buffered events" << endl; 
     37    } 
    3438  } 
    3539 

  • XIOS/dev/dev_trunk_omp/src/client.cpp

    r1628 r1646  
    1212#include "buffer_client.hpp" 
    1313#include "string_tools.hpp" 
     14#ifdef _usingEP 
    1415using namespace ep_lib; 
     16#endif 
    1517 
    1618namespace xios 
    … …  
    226228        MPI_Intercomm_create(contextComm,0,CXios::globalComm,serverLeader,10+globalRank,&contextInterComm) ; 
    227229        #pragma omp critical (_output) 
    228         info(10)<<"Register new Context : "<<id<<endl ; 
     230        { 
     231          info(10)<<"Register new Context : "<<id<<endl ; 
     232        } 
    229233        MPI_Comm inter ; 
    230234        MPI_Intercomm_merge(contextInterComm,0,&inter) ; 
    … …  
    303307      } 
    304308      #pragma omp critical (_output) 
    305       info(20) << "Client side context is finalized"<<endl ; 
     309      { 
     310        info(20) << "Client side context is finalized"<<endl ; 
     311      } 
    306312 
    307313      #pragma omp critical (_output) 

  • XIOS/dev/dev_trunk_omp/src/client_server_mapping.cpp

    r1601 r1646  
    99#include "client_server_mapping.hpp" 
    1010 
     11#ifdef _usingEP 
    1112using namespace ep_lib; 
     13#endif 
    1214 
    1315namespace xios { 

  • XIOS/dev/dev_trunk_omp/src/client_server_mapping_distributed.cpp

    r1601 r1646  
    1515#include "context.hpp" 
    1616#include "context_client.hpp" 
     17#ifdef _usingEP 
    1718using namespace ep_lib; 
     19#endif 
    1820 
    1921namespace xios 

  • XIOS/dev/dev_trunk_omp/src/config/domain_attribute.conf

    r1630 r1646  
    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

    r1630 r1646  
    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

    r1630 r1646  
    1212#include "cxios.hpp" 
    1313#include "server.hpp" 
     14#ifdef _usingEP 
    1415using namespace ep_lib; 
     16#endif 
    1517 
    1618namespace xios 
    … …  
    9698    { 
    9799      list<int> ranks = event.getRanks(); 
    98  
     100      #pragma omp critical (_output) 
     101      { 
     102        info(100)<<"Event "<<timeLine<<" of context "<<context->getId()<<endl ; 
     103      } 
    99104      if (CXios::checkEventSync) 
    100105      { 
    … …  
    124129        { 
    125130          event.send(timeLine, sizes, buffList); 
     131          #pragma omp critical (_output) 
     132          { 
     133            info(100)<<"Event "<<timeLine<<" of context "<<context->getId()<<"  sent"<<endl ; 
     134          } 
    126135 
    127136          checkBuffers(ranks); 
    … …  
    140149          for (list<int>::const_iterator it = sizes.begin(); it != sizes.end(); it++) 
    141150            tmpBufferedEvent.buffers.push_back(new CBufferOut(*it)); 
    142           info(100)<<"DEBUG : temporaly event created : timeline "<<timeLine<<endl ; 
     151          #pragma omp critical (_output) 
     152          { 
     153            info(100)<<"DEBUG : temporaly event created : timeline "<<timeLine<<endl ; 
     154          } 
    143155          event.send(timeLine, tmpBufferedEvent.sizes, tmpBufferedEvent.buffers); 
     156          #pragma omp critical (_output) 
     157          { 
     158            info(100)<<"Event "<<timeLine<<" of context "<<context->getId()<<"  sent"<<endl ; 
     159          } 
    144160        } 
    145161      } 
    … …  
    167183            (*itBuffer)->put((char*)(*it)->start(), (*it)->count()); 
    168184 
    169           info(100)<<"DEBUG : temporaly event sent "<<endl ; 
     185          #pragma omp critical (_output) 
     186          { 
     187            info(100)<<"DEBUG : temporaly event sent "<<endl ; 
     188          } 
    170189          checkBuffers(tmpBufferedEvent.ranks); 
    171190 
    … …  
    341360       if (ratio < minBufferSizeEventSizeRatio) minBufferSizeEventSizeRatio = ratio; 
    342361     } 
     362      
     363     #ifdef _usingEP 
    343364     MPI_Allreduce(&minBufferSizeEventSizeRatio, &minBufferSizeEventSizeRatio, 1, MPI_DOUBLE, MPI_MIN, intraComm); 
     365     #elif _usingMPI 
     366     MPI_Allreduce(MPI_IN_PLACE, &minBufferSizeEventSizeRatio, 1, MPI_DOUBLE, MPI_MIN, intraComm); 
     367     #endif 
    344368 
    345369     if (minBufferSizeEventSizeRatio < 1.0) 
    … …  
    426450      { 
    427451        #pragma omp critical (_output) 
    428         info(100)<<"DEBUG : Sent context Finalize event to rank "<<*itRank<<endl ; 
     452        { 
     453          info(100)<<"DEBUG : Sent context Finalize event to rank "<<*itRank<<endl ; 
     454        } 
    429455        event.push(*itRank, 1, msg); 
    430456      } 
    … …  
    452478    { 
    453479      #pragma omp critical (_output) 
    454       report(10) << " Memory report : Context <" << context->getId() << "> : client side : memory used for buffer of each connection to server" << endl 
     480      { 
     481        report(10) << " Memory report : Context <" << context->getId() << "> : client side : memory used for buffer of each connection to server" << endl 
    455482                 << "  +) To server with rank " << itMap->first << " : " << itMap->second << " bytes " << endl; 
     483      } 
    456484      totalBuf += itMap->second; 
    457485    } 
    458486    #pragma omp critical (_output) 
    459     report(0) << " Memory report : Context <" << context->getId() << "> : client side : total memory used for buffer " << totalBuf << " bytes" << endl; 
     487    { 
     488      report(0) << " Memory report : Context <" << context->getId() << "> : client side : total memory used for buffer " << totalBuf << " bytes" << endl; 
     489    } 
    460490 
    461491    //releaseBuffers(); // moved to CContext::finalize() 

  • XIOS/dev/dev_trunk_omp/src/context_server.cpp

    r1601 r1646  
    1818#include <boost/functional/hash.hpp> 
    1919 
     20#ifdef _usingEP 
    2021using namespace ep_lib; 
     22#endif 
    2123 
    2224namespace xios 
    … …  
    265267      finished=true; 
    266268      #pragma omp critical (_output) 
    267       info(20)<<" CContextServer: Receive context <"<<context->getId()<<"> finalize."<<endl; 
     269      { 
     270        info(20)<<" CContextServer: Receive context <"<<context->getId()<<"> finalize."<<endl; 
     271      } 
    268272      context->finalize(); 
    269273      std::map<int, StdSize>::const_iterator itbMap = mapBufferSize_.begin(), 
    … …  
    273277        rank = itMap->first; 
    274278        #pragma omp critical (_output) 
    275         report(10)<< " Memory report : Context <"<<ctxId<<"> : server side : memory used for buffer of each connection to client" << endl 
     279        { 
     280          report(10)<< " Memory report : Context <"<<ctxId<<"> : server side : memory used for buffer of each connection to client" << endl 
    276281            << "  +) With client of rank " << rank << " : " << itMap->second << " bytes " << endl; 
     282        } 
    277283        totalBuf += itMap->second; 
    278284      } 
    279285      #pragma omp critical (_output) 
    280       report(0)<< " Memory report : Context <"<<ctxId<<"> : server side : total memory used for buffer "<<totalBuf<<" bytes"<<endl; 
     286      { 
     287        report(0)<< " Memory report : Context <"<<ctxId<<"> : server side : total memory used for buffer "<<totalBuf<<" bytes"<<endl; 
     288      } 
    281289    } 
    282290    else if (event.classId==CContext::GetType()) CContext::dispatchEvent(event); 

  • XIOS/dev/dev_trunk_omp/src/cxios.cpp

    r1601 r1646  
    1111#include "memtrack.hpp" 
    1212#include "registry.hpp" 
     13#include "timer.hpp" 
     14#ifdef _usingEP 
    1315using namespace ep_lib; 
     16#endif 
    1417 
    1518namespace xios 
    … …  
    2124  const string CXios::serverPrmFile="./xios_server1"; 
    2225  const string CXios::serverSndFile="./xios_server2"; 
     26 
     27  bool CXios::xiosStack = true; 
     28  bool CXios::systemStack = false; 
    2329 
    2430  bool CXios::isClient ; 
    … …  
    4854    #pragma omp critical 
    4955    { 
    50       std::cout<<"thread "<<tmp_rank<<"("<<omp_get_thread_num()<<")"<<" parsing rootfile"<<std::endl; 
     56      //std::cout<<"thread "<<tmp_rank<<"("<<omp_get_thread_num()<<")"<<" parsing rootfile"<<std::endl; 
    5157      parseFile(rootFile); 
    5258      std::cout<<"thread "<<tmp_rank<<"("<<omp_get_thread_num()<<")"<<" parsed rootfile"<<std::endl; 
    … …  
    7177    printLogs2Files=getin<bool>("print_file",false); 
    7278 
     79    xiosStack=getin<bool>("xios_stack",true) ; 
     80    systemStack=getin<bool>("system_stack",false) ; 
     81    if (xiosStack && systemStack) 
     82    { 
     83      xiosStack = false; 
     84    } 
     85 
    7386    StdString bufMemory("memory"); 
    7487    StdString bufPerformance("performance"); 
    … …  
    90103    checkEventSync = getin<bool>("check_event_sync", checkEventSync); 
    91104 
    92     //globalComm=MPI_COMM_WORLD ; 
     105    #ifdef _usingMPI 
     106    globalComm=MPI_COMM_WORLD ; 
     107    #elif _usingEP 
    93108    int num_ep; 
    94     if(isClient)   
    95     {  
    96       num_ep = omp_get_num_threads(); 
    97     } 
    98          
    99     if(isServer)  
    100     {  
    101       num_ep = 1; 
    102     } 
     109 
     110    if(isClient)  num_ep = omp_get_num_threads(); 
     111    if(isServer) num_ep = 1; 
    103112         
    104113    MPI_Info info; 
    … …  
    111120         
    112121    #pragma omp barrier 
    113      
    114            
    115122    CXios::globalComm = passage[omp_get_thread_num()]; 
     123    #endif 
    116124  } 
    117125 
    … …  
    123131  */ 
    124132  void CXios::initClientSide(const string& codeId, MPI_Comm& localComm, MPI_Comm& returnComm) 
     133  TRY 
    125134  { 
    126135    isClient = true; 
    … …  
    148157    } 
    149158  } 
     159  CATCH 
    150160 
    151161  void CXios::clientFinalize(void) 
    … …  
    155165     { 
    156166       #pragma omp critical (_output) 
    157        info(80)<<"Write data base Registry"<<endl<<globalRegistry->toString()<<endl ; 
     167       { 
     168        info(80)<<"Write data base Registry"<<endl<<globalRegistry->toString()<<endl ; 
     169       } 
    158170       globalRegistry->toFile("xios_registry.bin") ; 
    159171       delete globalRegistry ; 
    … …  
    232244      if (CServer::getRank()==0) 
    233245      { 
    234         info(80)<<"Write data base Registry"<<endl<<globalRegistry->toString()<<endl ; 
     246        #pragma omp critical (_output) 
     247        { 
     248          info(80)<<"Write data base Registry"<<endl<<globalRegistry->toString()<<endl ; 
     249        } 
    235250        globalRegistry->toFile("xios_registry.bin") ; 
    236251        delete globalRegistry ; 
    … …  
    289304          } 
    290305 
    291           info(80)<<"Write data base Registry"<<endl<<globalRegistrySndServers.toString()<<endl ; 
     306          #pragma omp critical (_output) 
     307          { 
     308            info(80)<<"Write data base Registry"<<endl<<globalRegistrySndServers.toString()<<endl ; 
     309          } 
    292310          globalRegistrySndServers.toFile("xios_registry.bin") ; 
    293311 
    … …  
    296314      delete globalRegistry; 
    297315    } 
     316    CTimer::get("XIOS").suspend() ; 
    298317    CServer::finalize(); 
    299318 

  • XIOS/dev/dev_trunk_omp/src/cxios.hpp

    r1601 r1646  
    3333     static const string serverPrmFile;  //!< Filename template for primary server in case of two server levels 
    3434     static const string serverSndFile;  //!< Filename template for secondary server in case of two server levels 
     35 
     36     static bool xiosStack;    //!< Exception handling 
     37     static bool systemStack;  //!< Exception handling 
     38     #pragma omp threadprivate(xiosStack, systemStack) 
    3539 
    3640     static bool isClient ; //!< Check if xios is client 

  • XIOS/dev/dev_trunk_omp/src/data_output.cpp

    r1542 r1646  
    4747 
    4848      void CDataOutput::writeGrid(CDomain* domain, CAxis* axis) 
     49      TRY 
    4950      { 
    5051         this->writeDomain_(domain); 
    5152         this->writeAxis_(axis); 
    5253      } 
     54      CATCH 
    5355 
    5456      void CDataOutput::writeGrid(std::vector<CDomain*> domains, std::vector<CAxis*> axis) 
     57      TRY 
    5558      { 
    5659        int domSize = domains.size(); 
    … …  
    5962        for (int i = 0; i < aSize; ++i) this->writeAxis_(axis[i]); 
    6063      } 
     64      CATCH 
    6165 
    6266      void CDataOutput::writeGrid(std::vector<CDomain*> domains, std::vector<CAxis*> axis, std::vector<CScalar*> scalars) 
     67      TRY 
    6368      { 
    6469        int domSize = domains.size(); 
    … …  
    6974        for (int i = 0; i < sSize; ++i) this->writeScalar_(scalars[i]); 
    7075      } 
     76      CATCH 
    7177 
    7278      //---------------------------------------------------------------- 
    7379 
    7480      void CDataOutput::writeGrid(CDomain* domain) 
     81      TRY 
    7582      { 
    7683         this->writeDomain_(domain); 
    7784      } 
     85      CATCH 
    7886 
    7987      void CDataOutput::writeTimeDimension(void) 
     88      TRY 
    8089      { 
    8190         this->writeTimeDimension_(); 
    8291      } 
     92      CATCH 
    8393 
    8494      //---------------------------------------------------------------- 
    8595 
    8696      void CDataOutput::writeFieldTimeAxis(CField* field) 
     97      TRY 
    8798      { 
    8899         CContext* context = CContext::getCurrent() ; 
    … …  
    91102         this->writeTimeAxis_(field, calendar); 
    92103      } 
    93        
     104      CATCH 
     105 
    94106      void CDataOutput::writeField(CField* field) 
     107      TRY 
    95108      { 
    96109         this->writeField_(field); 
    97110      } 
     111      CATCH 
    98112 
    99113      //---------------------------------------------------------------- 
    100114 
    101115      void CDataOutput::writeFieldGrid(CField* field) 
     116      TRY 
    102117      { 
    103118         this->writeGrid(field->getRelGrid(), 
    104119                         !field->indexed_output.isEmpty() && field->indexed_output); 
    105120      } 
    106  
     121      CATCH 
    107122      //---------------------------------------------------------------- 
    108123 
    109124      void CDataOutput::writeFieldData(CField* field) 
     125      TRY 
    110126      { 
    111127//         CGrid* grid = CGrid::get(field->grid_ref.getValue()); 
    … …  
    113129         this->writeFieldData_(field); 
    114130      } 
     131      CATCH 
    115132 
    116133      ///---------------------------------------------------------------- 

  • XIOS/dev/dev_trunk_omp/src/dht_auto_indexing.cpp

    r1601 r1646  
    88 */ 
    99#include "dht_auto_indexing.hpp" 
     10#ifdef _usingEP 
    1011using namespace ep_lib; 
     12#endif 
    1113 
    1214namespace xios 

  • XIOS/dev/dev_trunk_omp/src/distribution_client.cpp

    r1593 r1646  
    1515   , axisDomainOrder_() 
    1616   , nLocal_(), nGlob_(), nBeginLocal_(), nBeginGlobal_() 
    17    , dataNIndex_(), dataDims_(), dataBegin_(), dataIndex_(), domainMasks_(), axisMasks_() 
     17   , dataNIndex_(), dataDims_(), dataBegin_(), dataIndex_() 
    1818   , gridMask_(), indexMap_() 
    1919   , isDataDistributed_(true), axisNum_(0), domainNum_(0) 
    … …  
    3636  GlobalLocalMap void2 ; 
    3737  std::vector<int> void3 ; 
    38   std::vector<int> void4 ; 
     38  std::vector<bool> void4 ; 
    3939 
    4040  globalLocalDataSendToServerMap_.swap(void1) ; 
    … …  
    6161  // Then check mask of grid 
    6262  int gridDim = domList.size() * 2 + axisList.size(); 
    63   grid->checkMask(); 
    6463  switch (gridDim) { 
    6564    case 0: 
    … …  
    117116  axisDomainOrder_ = axisDomainOrder; 
    118117 
    119   // Each domain or axis has its mask, of course 
    120   domainMasks_.resize(domainNum_); 
    121   for (int i = 0; i < domainNum_;++i) 
    122   { 
    123     domainMasks_[i].resize(domList[i]->domainMask.numElements()); 
    124     domainMasks_[i] = domList[i]->domainMask; 
    125   } 
    126  
    127   axisMasks_.resize(axisNum_); 
    128   for (int i = 0; i < axisNum_; ++i) 
    129   { 
    130     axisMasks_[i].resize(axisList[i]->mask.numElements()); 
    131     axisMasks_[i] = axisList[i]->mask; 
    132   } 
    133  
    134118  // Because domain and axis can be in any order (axis1, domain1, axis2, axis3, ) 
    135119  // their position should be specified. In axisDomainOrder, domain == true, axis == false 
    … …  
    282266 
    283267        if ((iIdx >= nBeginLocal_[indexMap_[i]]) && (iIdx < nLocal_[indexMap_[i]]) && 
    284            (jIdx >= nBeginLocal_[indexMap_[i]+1]) && (jIdx < nLocal_[indexMap_[i]+1]) && 
    285            (domainMasks_[idxDomain](iIdx + jIdx*nLocal_[indexMap_[i]]))) 
     268           (jIdx >= nBeginLocal_[indexMap_[i]+1]) && (jIdx < nLocal_[indexMap_[i]+1])) 
    286269        { 
    287270          ++count; 
    … …  
    325308      elementIndexData_[i].resize(dataNIndex_[i]); 
    326309      elementIndexData_[i] = false; 
    327       int iIdx = 0, count = 0, localIndex = 0; 
     310      int iIdx = 0, count = 0; 
    328311      for (int j = 0; j < dataNIndex_[i]; ++j) 
    329312      { 
    330313        iIdx = getAxisIndex((dataIndex_[indexMap_[i]])(j), dataBegin_[indexMap_[i]], nLocal_[indexMap_[i]]); 
    331314        if ((iIdx >= nBeginLocal_[indexMap_[i]]) && 
    332            (iIdx < nLocal_[indexMap_[i]]) && (axisMasks_[idxAxis](iIdx))) 
     315           (iIdx < nLocal_[indexMap_[i]]) )//&& (axisMasks_[idxAxis](iIdx))) 
    333316        { 
    334317          ++count; 
    … …  
    413396 
    414397  for (int i = 0; i < numElement_; ++i) ssize *= eachElementSize[i]; 
    415   while (idx < ssize) 
    416   { 
    417     for (int i = 0; i < numElement_-1; ++i) 
    418     { 
    419       if (idxLoop[i] == eachElementSize[i]) 
    420       { 
    421         idxLoop[i] = 0; 
    422         ++idxLoop[i+1]; 
    423       } 
    424     } 
    425  
    426     // Find out outer index 
    427     // Depending the inner-most element is axis or domain, 
    428     // The outer loop index begins correspondingly at one (1) or zero (0) 
    429     for (int i = 1; i < numElement_; ++i) 
    430     { 
    431       currentIndex[i] = elementLocalIndex_[i](idxLoop[i]); 
    432     } 
    433  
    434     // Inner most index 
    435     for (int i = 0; i < innerLoopSize; ++i) 
    436     { 
    437       int gridMaskIndex = 0; 
    438       currentIndex[0] = elementLocalIndex_[0](i); 
    439  
    440       // If defined, iterate on grid mask 
    441       if (!gridMask_.isEmpty()) 
    442       { 
    443         for (int k = 0; k < this->numElement_; ++k) 
    444         { 
    445           gridMaskIndex += (currentIndex[k])*elementNLocal_[k]; 
    446         } 
    447         if (gridMask_(gridMaskIndex)) ++indexLocalDataOnClientCount; 
    448       } 
    449       // If grid mask is not defined, iterate on elements' mask 
    450       else 
    451       { 
    452         bool maskTmp = true; 
    453         int idxDomain = 0, idxAxis = 0; 
    454         for (int elem = 0; elem < numElement_; ++elem) 
    455         { 
    456           if (2 == axisDomainOrder_(elem)) 
    457           { 
    458             maskTmp = maskTmp && domainMasks_[idxDomain](currentIndex[elem]); 
    459             ++idxDomain; 
    460           } 
    461           else if (1 == axisDomainOrder_(elem)) 
    462           { 
    463             maskTmp = maskTmp && axisMasks_[idxAxis](currentIndex[elem]); 
    464             ++idxAxis; 
    465           } 
    466         } 
    467         if (maskTmp) ++indexLocalDataOnClientCount; 
    468       } 
    469  
    470     } 
    471     idxLoop[0] += innerLoopSize; 
    472     idx += innerLoopSize; 
    473   } 
    474  
    475   // Now allocate these arrays 
    476   localDataIndex_.resize(indexLocalDataOnClientCount); 
    477   localMaskIndex_.resize(indexLocalDataOnClientCount); 
    478   localMaskedDataIndex_.resize(indexLocalDataOnClientCount); 
    479   globalDataIndex_.rehash(std::ceil(indexLocalDataOnClientCount/globalDataIndex_.max_load_factor())); 
    480   globalLocalDataSendToServerMap_.rehash(std::ceil(indexLocalDataOnClientCount/globalLocalDataSendToServerMap_.max_load_factor())); 
     398 
     399  localDataIndex_.resize(ssize); 
     400  if (!gridMask_.isEmpty()) localMaskIndex_.resize(ssize); 
     401  localMaskedDataIndex_.resize(ssize); 
     402  globalDataIndex_.rehash(std::ceil(ssize/globalDataIndex_.max_load_factor())); 
     403  globalLocalDataSendToServerMap_.rehash(std::ceil(ssize/globalLocalDataSendToServerMap_.max_load_factor())); 
     404 
    481405 
    482406  // We need to loop with data index 
    … …  
    535459        if (isCurrentIndexDataCorrect) 
    536460        { 
    537           int gridMaskIndex = 0; 
    538           for (int k = 0; k < this->numElement_; ++k) 
     461          bool maskTmp = true; 
     462          bool maskGridTmp = true; 
     463          size_t globalIndex = 0; 
     464          for (int k = 0; k < numElement_; ++k) 
    539465          { 
    540             gridMaskIndex += (currentIndex[k])*elementNLocal_[k]; 
     466            globalIndex += (currentGlobalIndex[k])*elementNGlobal_[k]; 
    541467          } 
    542  
    543           bool maskTmp = true; 
    544           // If defined, apply grid mask 
    545          if (!gridMask_.isEmpty()) 
     468          globalDataIndex_[globalIndex] = indexLocalDataOnClientCount; 
     469          localDataIndex_[indexLocalDataOnClientCount] = countLocalData; 
     470          globalLocalDataSendToServerMap_[globalIndex] = indexLocalDataOnClientCount; 
     471          localMaskedDataIndex_[indexLocalDataOnClientCount] = indexLocalDataOnClientCount; 
     472 
     473          // Grid mask: unmasked values will be replaces by NaN and then all values will be sent 
     474          if (!gridMask_.isEmpty()) 
    546475          { 
    547             maskTmp =  gridMask_(gridMaskIndex); 
     476            int gridMaskIndex = 0; 
     477            for (int k = 0; k < this->numElement_; ++k) 
     478            { 
     479              gridMaskIndex += (currentIndex[k])*elementNLocal_[k]; 
     480            } 
     481            maskGridTmp =  gridMask_(gridMaskIndex); 
     482            if (maskGridTmp) 
     483              localMaskIndex_[indexLocalDataOnClientCount] = true; 
     484            else 
     485              localMaskIndex_[indexLocalDataOnClientCount] = false; 
    548486          } 
    549           // If grid mask is not defined, apply elements' mask 
    550           else 
    551           { 
    552             int idxDomain = 0, idxAxis = 0; 
    553             for (int elem = 0; elem < numElement_; ++elem) 
    554             { 
    555               if (2 == axisDomainOrder_(elem)) 
    556               { 
    557                 maskTmp = maskTmp && domainMasks_[idxDomain](currentIndex[elem]); 
    558                 ++idxDomain; 
    559               } 
    560               else if (1 == axisDomainOrder_(elem)) 
    561               { 
    562                 maskTmp = maskTmp && axisMasks_[idxAxis](currentIndex[elem]); 
    563                 ++idxAxis; 
    564               } 
    565             } 
    566           } 
    567  
    568           if (maskTmp) 
    569           { 
    570             size_t globalIndex = 0; 
    571             for (int k = 0; k < numElement_; ++k) 
    572             { 
    573               globalIndex += (currentGlobalIndex[k])*elementNGlobal_[k]; 
    574             } 
    575             globalDataIndex_[globalIndex] = indexLocalDataOnClientCount; 
    576             localDataIndex_[indexLocalDataOnClientCount] = countLocalData; 
    577             globalLocalDataSendToServerMap_[globalIndex] = indexLocalDataOnClientCount; 
    578             localMaskIndex_[indexLocalDataOnClientCount] = gridMaskIndex; 
    579             localMaskedDataIndex_[indexLocalDataOnClientCount] = indexLocalDataOnClientCount; 
    580             ++indexLocalDataOnClientCount; 
    581           } 
     487 
     488          ++indexLocalDataOnClientCount; 
     489 
    582490        } 
    583491        ++countLocalData; 
    … …  
    586494    } 
    587495    else countLocalData+=innerLoopSize ; 
    588      
     496 
    589497    idxLoop[0] += innerLoopSize; 
    590498    idx += innerLoopSize; 
    … …  
    614522                                        const int& dataDim, const int& ni, int& j) 
    615523{ 
     524  int i; 
    616525  int tempI = dataIIndex + dataIBegin, 
    617526      tempJ = (dataJIndex + dataJBegin); 
    618527  if (ni == 0) 
    619528  { 
    620     int i = 0; 
    621     j = 0; 
     529    i = -1; 
     530    j = -1; 
    622531    return i; 
    623532  } 
    624   int i = (dataDim == 1) ? (tempI) % ni 
    625                          : (tempI) ; 
    626   j = (dataDim == 1) ? (tempI) / ni 
    627                      : (tempJ) ; 
    628  
     533  if ((tempI < 0) || (tempJ < 0)) 
     534  { 
     535    i = -1; 
     536    j = -1; 
     537    return i; 
     538  } 
     539  else 
     540  { 
     541    i = (dataDim == 1) ? (tempI) % ni : (tempI) ; 
     542    j = (dataDim == 1) ? (tempI) / ni : (tempJ) ; 
     543  } 
    629544  return i; 
    630545} 
    … …  
    643558    return -1; 
    644559  } 
    645   int tempI = dataIndex + dataBegin; 
     560  int tempI = dataIndex; 
    646561  if ((tempI < 0) || (tempI > ni)) 
    647562    return -1; 
    … …  
    677592  Return local mask index of client 
    678593*/ 
    679 const std::vector<int>& CDistributionClient::getLocalMaskIndexOnClient() 
     594const std::vector<bool>& CDistributionClient::getLocalMaskIndexOnClient() 
    680595{ 
    681596  if (!isComputed_) createGlobalIndexSendToServer(); 

  • XIOS/dev/dev_trunk_omp/src/distribution_client.hpp

    r1562 r1646  
    3333 
    3434  public: 
    35     /** Default constructor */         
     35    /** Default constructor */ 
    3636    CDistributionClient(int rank, CGrid* grid); 
    3737 
    … …  
    4444    GlobalLocalDataMap& getGlobalLocalDataSendToServer(); 
    4545    GlobalLocalDataMap& getGlobalDataIndexOnClient(); 
    46     const std::vector<int>& getLocalMaskIndexOnClient(); 
     46    const std::vector<bool>& getLocalMaskIndexOnClient(); 
    4747    const std::vector<int>& getLocalMaskedDataIndexOnClient(); 
    4848 
    … …  
    8383    GlobalLocalDataMap globalLocalDataSendToServerMap_; 
    8484    GlobalLocalDataMap globalDataIndex_; 
     85 
     86    /*! Array holding masked data indexes. 
     87     * It includes: 
     88     *  masking on data (data_i/j_index or data_ni/nj and data_ibegin) 
     89     *  masking on grid elements (domains, axes or scalars) 
     90     * It DOES NOT include grid mask. 
     91     * The array size defines the data size entering the workflow. It is used by source filter of client or server1. 
     92    */ 
    8593    std::vector<int> localDataIndex_; 
    86     std::vector<int> localMaskIndex_; 
     94 
     95    /*! Array holding grid mask. If grid mask is not defined, its size is zero. 
     96     * It is used by source filter of client for replacing unmasked data by NaN. 
     97    */ 
     98    std::vector<bool> localMaskIndex_; 
     99 
    87100    std::vector<int> localMaskedDataIndex_; 
    88101 
    … …  
    104117    std::vector<CArray<int,1> > dataIndex_; //!< Data index 
    105118    std::vector<CArray<int,1> > infoIndex_; //!< i_index, j_index 
    106  
    107     std::vector<CArray<bool,1> > domainMasks_; //!< Domain mask 
    108     std::vector<CArray<bool,1> > axisMasks_; //!< Axis mask 
    109119 
    110120    std::vector<int> indexMap_; //!< Mapping element index to dimension index 

  • XIOS/dev/dev_trunk_omp/src/event_client.cpp

    r1601 r1646  
    5353     { 
    5454       #pragma omp critical(_output) 
    55        info(100)<<"Send event "<<timeLine<<" classId : "<<classId<<"  typeId : "<<typeId<<endl ; 
     55       { 
     56        info(100)<<"Send event "<<timeLine<<" classId : "<<classId<<"  typeId : "<<typeId<<endl ; 
     57       } 
    5658     } 
    5759     for (; itBuff != buffers.end(); ++itBuff, ++itSizes, ++itSenders, ++itMsg) 

  • XIOS/dev/dev_trunk_omp/src/event_scheduler.cpp

    r1601 r1646  
    44#include "tracer.hpp" 
    55 
     6#ifdef _usingEP 
    67using namespace ep_lib; 
     8#endif 
    79 
    810namespace xios 

  • XIOS/dev/dev_trunk_omp/src/exception.cpp

    r828 r1646  
    55#include "client.hpp" 
    66#include "server.hpp" 
    7 #include "cxios.hpp" 
    87#include "log.hpp" 
    98 
    … …  
    1211  /// ////////////////////// Définitions ////////////////////// /// 
    1312   CException::CException(void) 
    14       : CObject(), desc_rethrow(true) 
     13      : CObject(), desc_rethrow(true), stream() 
    1514   { /* Ne rien faire de plus */ } 
    1615 
    17    CException::CException(const StdString & id) 
    18       : CObject(id), desc_rethrow(true) 
     16   CException::CException(const std::string & id) 
     17      : CObject(id), desc_rethrow(true), stream() 
    1918   { /* Ne rien faire de plus */ } 
    2019 
    2120   CException::CException(const CException & exception) 
    22       : std::basic_ios<char>() 
    23       , CObject(exception.getId()) 
    24       , StdOStringStream() 
    25       , desc_rethrow(false) 
     21//      : std::basic_ios<char>() 
     22       : CObject(exception.getId()) 
     23//      , StdOStringStream() 
     24//      , desc_rethrow(false) 
    2625   { (*this) << exception.str(); } 
    2726 
    2827   CException::~CException(void) 
    2928   { 
    30       if (desc_rethrow) 
    31 #ifdef __XIOS_NOABORT 
    32       { 
    33         throw (*this); 
    34       } 
    35 #else 
    36      { 
    37       error << this->getMessage() << std::endl; 
    38       MPI_Abort(CXios::globalComm, -1); //abort(); 
    39       } 
    40 #endif 
     29//      if (desc_rethrow) 
     30//#ifdef __XIOS_NOABORT 
     31//      { 
     32//        throw (*this); 
     33//      } 
     34//#else 
     35//     { 
     36//      error << this->getMessage() << std::endl; 
     37//        throw 4; 
     38//      MPI_Abort(CXios::globalComm, -1); //abort(); 
     39//      } 
     40//#endif 
    4141   } 
    4242 
    4343   //--------------------------------------------------------------- 
    4444 
    45    StdString CException::getMessage(void) const 
     45   std::string CException::getMessage(void) const 
    4646   { 
    47       StdOStringStream oss; 
    48       oss << "> Error [" << this->getId() << "] : " << this->str(); 
    49       return (oss.str()); 
     47//     StdOStringStream oss; 
     48//     oss << "> Error [" << this->getId() << "] : " << this->str(); 
     49//     return (oss.str()); 
     50     return (stream.str()); 
    5051   } 
    5152 
    5253   StdOStringStream &  CException::getStream(void) 
    53    { return (*boost::polymorphic_cast<StdOStringStream*>(this)); } 
     54//   { return (*boost::polymorphic_cast<StdOStringStream*>(this)); } 
     55   { return stream; } 
    5456 
    55    StdString CException::toString(void) const 
    56    { return (StdString(this->getMessage())); } 
     57   std::string CException::toString(void) const 
     58//   { return (std::string(this->getMessage())); } 
     59   { return stream.str(); } 
    5760 
    58    void CException::fromString(const StdString & str) 
    59    { this->str(str); } 
     61   void CException::fromString(const std::string & str) 
     62   { } 
     63//   { this->str(str); } 
    6064 
    6165   //--------------------------------------------------------------- 

  • XIOS/dev/dev_trunk_omp/src/exception.hpp

    r792 r1646  
    55#include "xios_spl.hpp" 
    66#include "object.hpp" 
     7#include <iomanip> 
     8#include <stdexcept> 
    79 
    810namespace xios 
    911{ 
    1012   /// ////////////////////// Déclarations ////////////////////// /// 
     13 
    1114   class CException 
    1215      : private CObject, public StdOStringStream 
    … …  
    3336         virtual void fromString(const StdString & str); 
    3437 
     38         struct StackInfo 
     39         { 
     40           StdString file; 
     41           StdString function; 
     42           int line; 
     43           StdString info; 
     44         }; 
     45 
     46         std::list<StackInfo> stack; 
     47 
    3548      private : 
    3649 
    3750         /// Propriétés /// 
    3851         bool desc_rethrow; // throw destructor 
     52         StdOStringStream stream; 
    3953 
    4054   }; // CException 
    … …  
    4357/// //////////////////////////// Macros //////////////////////////// /// 
    4458 
     59#define FILE_NAME (std::strrchr("/" __FILE__, '/') + 1) 
     60 
     61#define FUNCTION_NAME (StdString(BOOST_CURRENT_FUNCTION).length() > 100 ? \ 
     62                       StdString(BOOST_CURRENT_FUNCTION).substr(0,100).append("...)") : BOOST_CURRENT_FUNCTION) 
     63 
    4564#define INFO(x) \ 
    46    "In file \'" __FILE__ "\', line " << __LINE__ << " -> " x << std::endl; 
     65   "In file \""<< FILE_NAME <<"\", function \"" << BOOST_CURRENT_FUNCTION <<"\", line " << __LINE__ << " -> " x << std::endl; 
    4766 
    4867#ifdef __XIOS_DEBUG 
    … …  
    5271#endif 
    5372 
    54 #define ERROR(id, x) xios::CException(id).getStream() << INFO(x) 
     73#define ERROR(id, x)  \ 
     74{                     \ 
     75       xios::CException exc(id);                \ 
     76       exc.getStream() << INFO(x);              \ 
     77       error << exc.getMessage() << std::endl;  \ 
     78       throw exc;                               \ 
     79} 
     80 
     81#ifdef __XIOS_EXCEPTION 
     82  #define TRY                          \ 
     83    {                                  \ 
     84      int funcFirstLine = __LINE__;    \ 
     85      try 
     86  #define CATCH                              \ 
     87      catch(CException& e)                   \ 
     88      {                                      \ 
     89        CException::StackInfo stk;           \ 
     90        stk.file = FILE_NAME;                \ 
     91        stk.function = FUNCTION_NAME;        \ 
     92        stk.line = funcFirstLine;            \ 
     93        e.stack.push_back(stk);              \ 
     94        if (CXios::xiosStack)                \ 
     95          throw;                             \ 
     96       else                                  \ 
     97         throw 0;                            \ 
     98      }                                      \ 
     99    } 
     100  #define CATCH_DUMP_ATTR                    \ 
     101      catch(CException& e)                   \ 
     102      {                                      \ 
     103        CException::StackInfo stk;           \ 
     104        stk.info.append(StdString("Object id=\"" + this->getId()) + "\" object type=\"" + this->getName() + "\"\n"); \ 
     105        stk.info.append("*** XIOS attributes as defined in XML file(s) or via Fortran interface:\n");                \ 
     106        stk.info.append("[" + this->dumpXiosAttributes() + "]\n");     \ 
     107        stk.info.append("*** Additional information:\n");              \ 
     108        stk.info.append("[" + this->dumpClassAttributes() + "]\n");    \ 
     109        stk.file = FILE_NAME;                \ 
     110        stk.function = FUNCTION_NAME;        \ 
     111        stk.line = funcFirstLine;            \ 
     112        e.stack.push_back(stk);              \ 
     113        if (CXios::xiosStack)                \ 
     114          throw;                             \ 
     115       else                                  \ 
     116         throw 0;                            \ 
     117      }                                      \ 
     118    } 
     119  #define CATCH_DUMP_STACK                                           \ 
     120      catch(CException& e)                                           \ 
     121      {                                                              \ 
     122        CException::StackInfo stk;                                   \ 
     123        int i = 1;                                                   \ 
     124        stk.file = FILE_NAME;                                        \ 
     125        stk.function = FUNCTION_NAME;                                \ 
     126        stk.line = funcFirstLine;                                    \ 
     127        e.stack.push_back(stk);                                      \ 
     128        for (auto itr = e.stack.crbegin(); itr!=e.stack.crend(); ++itr) \ 
     129        {                                                            \ 
     130          error << "("<< i <<") **************** ";                  \ 
     131          error << itr->function << std::endl;                       \ 
     132          error << itr->info << std::endl;                           \ 
     133          ++i;                                                       \ 
     134        }                                                            \ 
     135        error << left << "      ";                                   \ 
     136        error << left << std::setw(40) << "File " ;                  \ 
     137        error << left << std::setw(106) << "Function " ;             \ 
     138        error << std::setw(5)   << "Line " << std::endl;             \ 
     139        i = e.stack.size();                                          \ 
     140        for (auto itr = e.stack.begin(); itr!=e.stack.end(); itr++)  \ 
     141        {                                                            \ 
     142          StdOStringStream tmp;                                      \ 
     143          tmp <<  "("<< i <<")";                                     \ 
     144          error << left << std::setw(6) << tmp.str();                \ 
     145          error << left << std::setw(40)  << itr->file;              \ 
     146          error << left << std::setw(106) << itr->function;          \ 
     147          error << left << std::setw(5)   << itr->line << std::endl; \ 
     148          --i;                                                       \ 
     149        }                                                            \ 
     150        throw;                                                       \ 
     151      }                                                              \ 
     152    } 
     153#else 
     154  #define TRY 
     155  #define CATCH 
     156  #define CATCH_DUMP_ATTR 
     157  #define CATCH_DUMP_STACK 
     158#endif 
    55159 
    56160#endif // __XIOS_CException__ 
     161 

  • XIOS/dev/dev_trunk_omp/src/filter/file_writer_filter.cpp

    r1474 r1646  
    1717  void CFileWriterFilter::onInputReady(std::vector<CDataPacketPtr> data) 
    1818  { 
    19     const bool detectMissingValue = (!field->detect_missing_value.isEmpty() 
    20                                       && !field->default_value.isEmpty() 
    21                                       && field->detect_missing_value == true); 
     19    const bool detectMissingValue = ( !field->default_value.isEmpty() && 
     20                               ( (!field->detect_missing_value.isEmpty() || field->detect_missing_value == true) 
     21                                 || field->hasGridMask()) ); 
    2222 
    2323    CArray<double, 1> dataArray = (detectMissingValue) ? data[0]->data.copy() : data[0]->data; 

  • XIOS/dev/dev_trunk_omp/src/filter/source_filter.cpp

    r1250 r1646  
    77namespace xios 
    88{ 
    9   CSourceFilter::CSourceFilter(CGarbageCollector& gc, CGrid* grid, bool compression,  
     9  CSourceFilter::CSourceFilter(CGarbageCollector& gc, CGrid* grid, 
     10                               bool compression /*= true*/, bool mask /*= false*/, 
    1011                               const CDuration offset /*= NoneDu*/, bool manualTrigger /*= false*/, 
    1112                               bool hasMissingValue /*= false*/, 
    … …  
    1415    , grid(grid) 
    1516    , compression(compression) 
     17    , mask(mask) 
    1618    , offset(offset) 
    1719    , hasMissingValue(hasMissingValue), defaultValue(defaultValue) 
    … …  
    4042    } 
    4143    else 
    42       grid->inputField(data, packet->data); 
    43  
    44      
    45      
    46     // if (compression) grid->inputField(data, packet->data) ; 
    47     // else 
    48     // { 
    49     //   // just make a flat copy 
    50     //   CArray<double, N> data_tmp(data.copy()) ; // supress const attribute 
    51     //   CArray<double,1> dataTmp2(data_tmp.dataFirst(),shape(data.numElements()),neverDeleteData) ; 
    52     //   packet->data = dataTmp2 ; 
    53     // } 
     44    { 
     45      if (mask) 
     46        grid->maskField(data, packet->data); 
     47      else 
     48        grid->inputField(data, packet->data); 
     49    } 
    5450    // Convert missing values to NaN 
    5551    if (hasMissingValue) 

  • XIOS/dev/dev_trunk_omp/src/filter/source_filter.hpp

    r1241 r1646  
    2121       * \param gc the garbage collector associated with this filter 
    2222       * \param grid the grid to which the data is attached 
     23       * \param compression 
     24       * \param mask 
    2325       * \param offset the offset applied to the timestamp of all packets 
    2426       * \param manualTrigger whether the output should be triggered manually 
    … …  
    2729       */ 
    2830      CSourceFilter(CGarbageCollector& gc, CGrid* grid, 
    29                     bool compression=true, 
     31                    bool compression = true, 
     32                    bool mask = false, 
    3033                    const CDuration offset = NoneDu, bool manualTrigger = false, 
    3134                    bool hasMissingValue = false, 
    … …  
    6164 
    6265    private: 
    63       CGrid* grid; //!< The grid attached to the data the filter can accept 
    64       const CDuration offset; //!< The offset applied to the timestamp of all packets 
     66      CGrid* grid;             //!< The grid attached to the data the filter can accept 
     67      const CDuration offset;  //!< The offset applied to the timestamp of all packets 
    6568      const bool hasMissingValue; 
    6669      const double defaultValue; 
    67       const bool compression ; //!< indicate if the data need to be compressed : on client size : true, on server side : false 
     70      const bool compression ; //!< indicates if data need to be compressed : on client side : true, on server side : false 
     71      const bool mask ;        //!< indicates whether grid mask should be applied (true for clients, false for servers) 
    6872  }; // class CSourceFilter 
    6973} // namespace xios 

  • XIOS/dev/dev_trunk_omp/src/filter/spatial_transform_filter.cpp

    r1601 r1646  
    55#include "context_client.hpp" 
    66#include "timer.hpp" 
     7#ifdef _usingEP 
    78using namespace ep_lib; 
     9#endif 
    810 
    911namespace xios 
    … …  
    197199     
    198200    CContextClient* client = CContext::getCurrent()->client; 
     201    int rank; 
     202    MPI_Comm_rank (client->intraComm, &rank); 
    199203 
    200204    // Get default value for output data 
    … …  
    206210    const std::list<CGridTransformation::RecvIndexGridDestinationMap>& listLocalIndexToReceive = gridTransformation->getLocalIndexToReceiveOnGridDest(); 
    207211    const std::list<size_t>& listNbLocalIndexToReceive = gridTransformation->getNbLocalIndexToReceiveOnGridDest(); 
    208     const std::list<std::vector<bool> >& listLocalIndexMaskOnDest = gridTransformation->getLocalMaskIndexOnGridDest(); 
    209212    const std::vector<CGenericAlgorithmTransformation*>& listAlgos = gridTransformation->getAlgos(); 
    210213 
    … …  
    215218    std::list<CGridTransformation::RecvIndexGridDestinationMap>::const_iterator itListRecv = listLocalIndexToReceive.begin(); 
    216219    std::list<size_t>::const_iterator itNbListRecv = listNbLocalIndexToReceive.begin(); 
    217     std::list<std::vector<bool> >::const_iterator itLocalMaskIndexOnDest = listLocalIndexMaskOnDest.begin(); 
    218220    std::vector<CGenericAlgorithmTransformation*>::const_iterator itAlgo = listAlgos.begin(); 
    219221 
    220     for (; itListSend != iteListSend; ++itListSend, ++itListRecv, ++itNbListRecv, ++itLocalMaskIndexOnDest, ++itAlgo) 
     222    for (; itListSend != iteListSend; ++itListSend, ++itListRecv, ++itNbListRecv, ++itAlgo) 
    221223    { 
    222224      CArray<double,1> dataCurrentSrc(dataCurrentDest); 
    … …  
    228230      int idxSendBuff = 0; 
    229231      std::vector<double*> sendBuff(localIndexToSend.size()); 
     232      double* sendBuffRank; 
    230233      for (itSend = itbSend; itSend != iteSend; ++itSend, ++idxSendBuff) 
    231234      { 
     235        int destRank = itSend->first; 
    232236        if (0 != itSend->second.numElements()) 
    233           sendBuff[idxSendBuff] = new double[itSend->second.numElements()]; 
     237        { 
     238          if (rank != itSend->first) 
     239            sendBuff[idxSendBuff] = new double[itSend->second.numElements()]; 
     240          else 
     241            sendBuffRank = new double[itSend->second.numElements()]; 
     242        } 
    234243      } 
    235244 
    … …  
    242251        const CArray<int,1>& localIndex_p = itSend->second; 
    243252        int countSize = localIndex_p.numElements(); 
    244         for (int idx = 0; idx < countSize; ++idx) 
    245         { 
    246           sendBuff[idxSendBuff][idx] = dataCurrentSrc(localIndex_p(idx)); 
    247         } 
    248         MPI_Isend(sendBuff[idxSendBuff], countSize, MPI_DOUBLE, destRank, 12, client->intraComm, &sendRecvRequest[position++]); 
     253        if (destRank != rank) 
     254        { 
     255          for (int idx = 0; idx < countSize; ++idx) 
     256          { 
     257            sendBuff[idxSendBuff][idx] = dataCurrentSrc(localIndex_p(idx)); 
     258          }   
     259          MPI_Isend(sendBuff[idxSendBuff], countSize, MPI_DOUBLE, destRank, 12, client->intraComm, &sendRecvRequest[position++]); 
     260           
     261        } 
     262        else 
     263        { 
     264          for (int idx = 0; idx < countSize; ++idx) 
     265          { 
     266            sendBuffRank[idx] = dataCurrentSrc(localIndex_p(idx)); 
     267          } 
     268        } 
    249269      } 
    250270 
    … …  
    254274                                                                       iteRecv = localIndexToReceive.end(); 
    255275      int recvBuffSize = 0; 
    256       for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv) recvBuffSize += itRecv->second.size(); //(recvBuffSize < itRecv->second.size()) 
    257                                                                        //? itRecv->second.size() : recvBuffSize; 
     276      for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv) 
     277      { 
     278        if (itRecv->first != rank ) 
     279          recvBuffSize += itRecv->second.size(); 
     280      } 
     281      //(recvBuffSize < itRecv->second.size()) ? itRecv->second.size() : recvBuffSize; 
    258282      double* recvBuff; 
     283 
    259284      if (0 != recvBuffSize) recvBuff = new double[recvBuffSize]; 
    260285      int currentBuff = 0; 
    … …  
    262287      { 
    263288        int srcRank = itRecv->first; 
    264         int countSize = itRecv->second.size(); 
    265         MPI_Irecv(recvBuff + currentBuff, countSize, MPI_DOUBLE, srcRank, 12, client->intraComm, &sendRecvRequest[position++]); 
    266         currentBuff += countSize; 
     289        if (srcRank != rank) 
     290        { 
     291          int countSize = itRecv->second.size(); 
     292          MPI_Irecv(recvBuff + currentBuff, countSize, MPI_DOUBLE, srcRank, 12, client->intraComm, &sendRecvRequest[position++]); 
     293          currentBuff += countSize; 
     294        } 
    267295      } 
    268296      std::vector<MPI_Status> status(sendRecvRequest.size()); 
    269       MPI_Waitall(sendRecvRequest.size(), &sendRecvRequest[0], &status[0]); 
     297      MPI_Waitall(position, &sendRecvRequest[0], &status[0]); 
     298 
     299 
    270300 
    271301      dataCurrentDest.resize(*itNbListRecv); 
    272       const std::vector<bool>& localMaskDest = *itLocalMaskIndexOnDest; 
    273       for (int i = 0; i < localMaskDest.size(); ++i) 
    274         if (localMaskDest[i]) dataCurrentDest(i) = 0.0; 
    275         else dataCurrentDest(i) = defaultValue; 
     302      dataCurrentDest = 0.0; 
    276303 
    277304      std::vector<bool> localInitFlag(dataCurrentDest.numElements(), true); 
    … …  
    280307      for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv) 
    281308      { 
    282         int countSize = itRecv->second.size(); 
    283309        const std::vector<std::pair<int,double> >& localIndex_p = itRecv->second; 
    284         (*itAlgo)->apply(localIndex_p, 
    285                          recvBuff+currentBuff, 
    286                          dataCurrentDest, 
    287                          localInitFlag, 
    288                          ignoreMissingValue,firstPass); 
    289  
    290         currentBuff += countSize; 
     310        int srcRank = itRecv->first; 
     311        if (srcRank != rank) 
     312        { 
     313          int countSize = itRecv->second.size(); 
     314          (*itAlgo)->apply(localIndex_p, 
     315                           recvBuff+currentBuff, 
     316                           dataCurrentDest, 
     317                           localInitFlag, 
     318                           ignoreMissingValue,firstPass); 
     319          currentBuff += countSize; 
     320        } 
     321        else 
     322        { 
     323          (*itAlgo)->apply(localIndex_p, 
     324                           sendBuffRank, 
     325                           dataCurrentDest, 
     326                           localInitFlag, 
     327                           ignoreMissingValue,firstPass); 
     328        } 
     329 
    291330        firstPass=false ; 
    292331      } 
    … …  
    298337      { 
    299338        if (0 != itSend->second.numElements()) 
    300           delete [] sendBuff[idxSendBuff]; 
     339        { 
     340          if (rank != itSend->first) 
     341            delete [] sendBuff[idxSendBuff]; 
     342          else 
     343            delete [] sendBuffRank; 
     344        } 
    301345      } 
    302346      if (0 != recvBuffSize) delete [] recvBuff; 

  • XIOS/dev/dev_trunk_omp/src/filter/store_filter.cpp

    r1474 r1646  
    4848    { 
    4949      std::map<Time, CDataPacketPtr>::const_iterator it ; 
    50       info(0)<<"Impossible to get the packet with timestamp = " << timestamp<<std::endl<<"Available timestamp are : "<<std::endl ; 
    51       for(it=packets.begin();it!=packets.end();++it) info(0)<<it->first<<"  "; 
    52       info(0)<<std::endl ; 
     50      #pragma omp critical (_output) 
     51      { 
     52        info(0)<<"Impossible to get the packet with timestamp = " << timestamp<<std::endl<<"Available timestamp are : "<<std::endl ; 
     53      } 
     54      for(it=packets.begin();it!=packets.end();++it) 
     55      { 
     56        #pragma omp critical (_output) 
     57        { 
     58          info(0)<<it->first<<"  "; 
     59        } 
     60      } 
     61      #pragma omp critical (_output) 
     62      {  
     63        info(0)<<std::endl ; 
     64      } 
    5365      ERROR("CConstDataPacketPtr CStoreFilter::getPacket(Time timestamp) const", 
    5466            << "Impossible to get the packet with timestamp = " << timestamp); 

  • XIOS/dev/dev_trunk_omp/src/generate_fortran_interface.cpp

    r1558 r1646  
    317317  file.open((path+"ireorder_domain_attr.F90").c_str()); 
    318318  reorderDomain.generateFortranInterface(file); 
    319  
    320   file.open((path+"extract_domain_interface_attr.F90").c_str()); 
    321   extractDomain.generateFortran2003Interface(file); 
    322   file.close(); 
    323  
    324   file.open((path+"icextract_domain_attr.cpp").c_str()); 
    325   extractDomain.generateCInterface(file); 
    326   file.close(); 
    327  
    328   file.open((path+"iextract_domain_attr.F90").c_str()); 
    329   extractDomain.generateFortranInterface(file); 
    330  
    331319  file.close(); 
    332320   

  • XIOS/dev/dev_trunk_omp/src/generate_interface_impl.hpp

    r1158 r1646  
    215215    string fortranKindC=getStrFortranKindC<T>(); 
    216216 
    217     oss << "SUBROUTINE cxios_set_" << className << "_" << name << "(" << className << "_hdl, " << name << ") BIND(C)" << iendl; 
     217//    oss << "SUBROUTINE cxios_set_" << className << "_" << name << "(" << className << "_hdl, " << name << ") BIND(C)" << iendl; 
     218    int indent = oss.iword(iendl.index); 
     219    string str = "SUBROUTINE cxios_set_" + className + "_" + name + "(" + className + "_hdl, " + name + ") BIND(C)"; 
     220    if ((str.length() + indent) >132) 
     221    { 
     222      oss << str.substr(0,130-indent) ; 
     223      oss << "&" << endl; 
     224      oss << "&" << str.substr(130-indent,str.length()) ; 
     225    } 
     226    else 
     227    { 
     228      oss << str; 
     229    } 
     230    oss << iendl; 
    218231    oss << "  USE ISO_C_BINDING" << iendl; 
    219232    oss << "  INTEGER (kind = C_INTPTR_T), VALUE :: " << className << "_hdl" << iendl; 
    … …  
    231244  void CInterface::AttributeFortran2003Interface<string>(ostream& oss, const string& className, const string& name) 
    232245  { 
    233     oss << "SUBROUTINE cxios_set_" << className << "_" << name << "(" << className << "_hdl, " << name << ", " << name << "_size) BIND(C)" << iendl; 
     246//    oss << "SUBROUTINE cxios_set_" << className << "_" << name << "(" << className << "_hdl, " << name << ", " << name << "_size) BIND(C)" << iendl; 
     247    int indent = oss.iword(iendl.index); 
     248    string str ="SUBROUTINE cxios_set_" + className + "_" + name + "(" + className + "_hdl, " + name + ", " + name + "_size) BIND(C)"; 
     249    if ((str.length() + indent) >132) 
     250    { 
     251      oss << str.substr(0,130-indent) ; 
     252      oss << "&" << endl; 
     253      oss << "&" << str.substr(130-indent,str.length()) ; 
     254    } 
     255    else 
     256    { 
     257      oss << str; 
     258    } 
     259    oss << iendl; 
    234260    oss << "  USE ISO_C_BINDING" << iendl; 
    235261    oss << "  INTEGER (kind = C_INTPTR_T), VALUE :: " << className << "_hdl" << iendl; 
    … …  
    238264    oss << "END SUBROUTINE cxios_set_" << className << "_" << name << std::endl; 
    239265    oss << iendl; 
    240     oss << "SUBROUTINE cxios_get_" << className << "_" << name << "(" << className << "_hdl, " << name << ", " << name << "_size) BIND(C)" << iendl; 
     266//    oss << "SUBROUTINE cxios_get_" << className << "_" << name << "(" << className << "_hdl, " << name << ", " << name << "_size) BIND(C)" << iendl; 
     267    str = "SUBROUTINE cxios_get_" + className + "_" + name + "(" + className + "_hdl, " + name + ", " + name + "_size) BIND(C)"; 
     268    if ((str.length() + indent) >132) 
     269    { 
     270      oss << str.substr(0,130-indent) ; 
     271      oss << "&" << endl; 
     272      oss << "&" << str.substr(130-indent,str.length()) ; 
     273    } 
     274    else 
     275    { 
     276      oss << str; 
     277    } 
     278    oss << iendl; 
    241279    oss << "  USE ISO_C_BINDING" << iendl; 
    242280    oss << "  INTEGER (kind = C_INTPTR_T), VALUE :: " << className << "_hdl" << iendl; 
    … …  
    247285 
    248286  template <> 
    249   void CInterface::AttributeFortran2003Interface<CDate>(ostream& oss, const string& className, const string& name) 
    250   { 
    251     oss << "SUBROUTINE cxios_set_" << className << "_" << name << "(" << className << "_hdl, " << name << ") BIND(C)" << iendl; 
    252     oss << "  USE ISO_C_BINDING" << iendl; 
    253     oss << "  USE IDATE" << iendl; 
    254     oss << "  INTEGER (kind = C_INTPTR_T), VALUE :: " << className << "_hdl" << iendl; 
    255     oss << "  TYPE(txios(date)), VALUE :: " << name << iendl; 
    256     oss << "END SUBROUTINE cxios_set_" << className << "_" << name << std::endl; 
     287  void CInterface::AttributeFortran2003Interface<CDuration>(ostream& oss, const string& className, const string& name) 
     288  { 
     289//    oss << "SUBROUTINE cxios_set_" << className << "_" << name << "(" << className << "_hdl, " << name << ") BIND(C)" << iendl; 
     290    string str = "SUBROUTINE cxios_set_" + className + "_" + name + "(" + className + "_hdl, " + name + ") BIND(C)"; 
     291    int indent = oss.iword(iendl.index); 
     292    if ((str.length() + indent) >132) 
     293    { 
     294      oss << str.substr(0,130-indent) ; 
     295      oss << "&" << endl; 
     296      oss << "&" << str.substr(130-indent,str.length()) ; 
     297    } 
     298    else 
     299    { 
     300      oss << str; 
     301    } 
    257302    oss << iendl; 
    258     oss << "SUBROUTINE cxios_get_" << className << "_" << name << "(" << className << "_hdl, " << name << ") BIND(C)" << iendl; 
    259     oss << "  USE ISO_C_BINDING" << iendl; 
    260     oss << "  USE IDATE" << iendl; 
    261     oss << "  INTEGER (kind = C_INTPTR_T), VALUE :: " << className << "_hdl" << iendl; 
    262     oss << "  TYPE(txios(date)) :: " << name << iendl; 
    263     oss << "END SUBROUTINE cxios_get_" << className << "_" << name << std::endl; 
    264   } 
    265  
    266   template <> 
    267   void CInterface::AttributeFortran2003Interface<CDuration>(ostream& oss, const string& className, const string& name) 
    268   { 
    269     oss << "SUBROUTINE cxios_set_" << className << "_" << name << "(" << className << "_hdl, " << name << ") BIND(C)" << iendl; 
    270303    oss << "  USE ISO_C_BINDING" << iendl; 
    271304    oss << "  USE IDURATION" << iendl; 

  • XIOS/dev/dev_trunk_omp/src/interface/c/icaxis.cpp

    r1542 r1646  
    2828    
    2929   void cxios_axis_handle_create (XAxisPtr * _ret, const char * _id, int _id_len) 
     30   TRY 
    3031   { 
    3132      std::string id;  
    … …  
    3536      CTimer::get("XIOS").suspend() ; 
    3637   } 
     38   CATCH_DUMP_STACK 
    3739    
    3840   void cxios_axisgroup_handle_create (XAxisGroupPtr * _ret, const char * _id, int _id_len) 
     41   TRY 
    3942   { 
    4043      std::string id;  
    … …  
    4447      CTimer::get("XIOS").suspend() ; 
    4548    } 
     49   CATCH_DUMP_STACK 
    4650 
    4751   // -------------------- Vérification des identifiants ----------------------- 
    4852 
    4953   void cxios_axis_valid_id (bool * _ret, const char * _id, int _id_len) 
     54   TRY 
    5055   { 
    5156      std::string id; 
    … …  
    5661      CTimer::get("XIOS").suspend() ; 
    5762   } 
     63   CATCH_DUMP_STACK 
    5864 
    5965   void cxios_axisgroup_valid_id (bool * _ret, const char * _id, int _id_len) 
     66   TRY 
    6067   { 
    6168      std::string id; 
    … …  
    6774 
    6875   } 
     76   CATCH_DUMP_STACK 
    6977    
    7078} // extern "C" 

  • XIOS/dev/dev_trunk_omp/src/interface/c/iccalendar.cpp

    r1542 r1646  
    1010{ 
    1111  void cxios_update_calendar(int step) 
     12  TRY 
    1213  { 
    1314    CTimer::get("XIOS").resume(); 
    … …  
    1920    CTimer::get("XIOS").suspend(); 
    2021  } 
     22  CATCH_DUMP_STACK 
    2123 
    2224  void cxios_get_current_date(cxios_date* current_date_c) 
     25  TRY 
    2326  { 
    2427    CTimer::get("XIOS").resume(); 
    … …  
    3740    CTimer::get("XIOS").suspend(); 
    3841  } 
     42  CATCH_DUMP_STACK 
    3943 
    4044  int cxios_get_year_length_in_seconds(int year) 
     45  TRY 
    4146  { 
    4247    CTimer::get("XIOS").resume(); 
    … …  
    5055    return length; 
    5156  } 
     57  CATCH_DUMP_STACK 
    5258 
    5359  int cxios_get_day_length_in_seconds() 
     60  TRY 
    5461  { 
    5562    CTimer::get("XIOS").resume(); 
    … …  
    6370    return length; 
    6471  } 
     72  CATCH_DUMP_STACK 
    6573} 

  • XIOS/dev/dev_trunk_omp/src/interface/c/iccalendar_wrapper.cpp

    r1542 r1646  
    2929 
    3030  void cxios_calendar_wrapper_handle_create(XCalendarWrapperPtr* _ret, const char* _id, int _id_len) 
     31  TRY 
    3132  { 
    3233    std::string id; 
    … …  
    3637    CTimer::get("XIOS").suspend(); 
    3738  } 
     39  CATCH_DUMP_STACK 
    3840 
    3941  void cxios_get_current_calendar_wrapper(XCalendarWrapperPtr* _ret) 
     42  TRY 
    4043  { 
    4144    CTimer::get("XIOS").resume(); 
    … …  
    4346    CTimer::get("XIOS").suspend(); 
    4447  } 
     48  CATCH_DUMP_STACK 
    4549 
    4650  // -------------------- Vérification des identifiants ----------------------- 
    4751 
    4852  void cxios_calendar_wrapper_valid_id(bool* _ret, const char* _id, int _id_len) 
     53  TRY 
    4954  { 
    5055    std::string id; 
    … …  
    5459    CTimer::get("XIOS").suspend(); 
    5560  } 
     61  CATCH_DUMP_STACK 
    5662 
    5763  // ----------------------- Custom getters and setters ----------------------- 
    5864 
    5965  void cxios_set_calendar_wrapper_date_start_date(XCalendarWrapperPtr calendarWrapper_hdl, cxios_date start_date_c) 
     66  TRY 
    6067  { 
    6168    CTimer::get("XIOS").resume(); 
    … …  
    7077    CTimer::get("XIOS").suspend(); 
    7178  } 
     79  CATCH_DUMP_STACK 
    7280 
    7381  void cxios_get_calendar_wrapper_date_start_date(XCalendarWrapperPtr calendarWrapper_hdl, cxios_date* start_date_c) 
     82  TRY 
    7483  { 
    7584    CTimer::get("XIOS").resume(); 
    … …  
    8392    CTimer::get("XIOS").suspend(); 
    8493  } 
     94  CATCH_DUMP_STACK 
    8595 
    8696  void cxios_set_calendar_wrapper_date_time_origin(XCalendarWrapperPtr calendarWrapper_hdl, cxios_date time_origin_c) 
     97  TRY 
    8798  { 
    8899    CTimer::get("XIOS").resume(); 
    … …  
    97108    CTimer::get("XIOS").suspend(); 
    98109  } 
     110  CATCH_DUMP_STACK 
    99111 
    100112  void cxios_get_calendar_wrapper_date_time_origin(XCalendarWrapperPtr calendarWrapper_hdl, cxios_date* time_origin_c) 
     113  TRY 
    101114  { 
    102115    CTimer::get("XIOS").resume(); 
    … …  
    110123    CTimer::get("XIOS").suspend(); 
    111124  } 
     125  CATCH_DUMP_STACK 
    112126 
    113127  // ----------------------- Calendar creation and update ---------------------- 
    114128 
    115129  void cxios_create_calendar(XCalendarWrapperPtr calendarWrapper_hdl) 
     130  TRY 
    116131  { 
    117132    CTimer::get("XIOS").resume(); 
    … …  
    119134    CTimer::get("XIOS").suspend(); 
    120135  } 
     136  CATCH_DUMP_STACK 
    121137 
    122138  void cxios_update_calendar_timestep(XCalendarWrapperPtr calendarWrapper_hdl) 
     139  TRY 
    123140  { 
    124141    CTimer::get("XIOS").resume(); 
    … …  
    126143    CTimer::get("XIOS").suspend(); 
    127144  } 
     145  CATCH_DUMP_STACK 
    128146} // extern "C" 

  • XIOS/dev/dev_trunk_omp/src/interface/c/iccompute_connectivity_domain.cpp

    r1542 r1646  
    2525   // ------------------------ Création des handle ----------------------------- 
    2626   void cxios_compute_connectivity_domain_handle_create(XComConDomainPtr * _ret, const char * _id, int _id_len) 
     27   TRY 
    2728   { 
    2829      std::string id; 
    … …  
    3233      CTimer::get("XIOS").suspend() ; 
    3334   } 
     35   CATCH_DUMP_STACK 
    3436 
    3537   // -------------------- Vérification des identifiants ----------------------- 
    3638   void cxios_compute_connectivity_domain_valid_id(bool * _ret, const char * _id, int _id_len) 
     39   TRY 
    3740   { 
    3841      std::string id; 
    … …  
    4346      CTimer::get("XIOS").suspend() ; 
    4447   } 
     48   CATCH_DUMP_STACK 
    4549} // extern "C" 

  • XIOS/dev/dev_trunk_omp/src/interface/c/iccontext.cpp

    r1542 r1646  
    3232 
    3333   void cxios_context_handle_create (XContextPtr * _ret, const char * _id, int _id_len) 
     34   TRY 
    3435   { 
    3536      std::string id; 
    … …  
    5455      // Lever une exeception ici 
    5556   } 
     57   CATCH_DUMP_STACK 
    5658 
    5759   // ------------------------ Changements de contextes ------------------------ 
    5860</