Changeset 1630


Ignore:
Timestamp:
12/21/18 09:19:12 (5 years ago)
Author:
yushan
Message:

working branch @1608 with bug fix @1628

Location:
XIOS/dev/dev_trunk_omp
Files:
22 edited

Legend:

Unmodified
Added
Removed

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

    r1619 r1630  
    4848        int n; /* number of vertices */ 
    4949        double area; 
    50   double given_area ; 
    5150        Coord x; /* barycentre */ 
    5251}; 
    … …  
    8180                n    = rhs.n; 
    8281                area = rhs.area; 
    83                 given_area = rhs.given_area; 
    8482                x    = rhs.x; 
    8583                val  = rhs.val; 

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

    r1619 r1630  
    4343        assert(n_cell_dst >= 4); 
    4444        assert(1 <= order && order <= 2); 
    45   double* src_area=NULL ; 
    46   double* dst_area=NULL ; 
     45 
    4746  mapper = new Mapper(MPI_COMM_WORLD); 
    4847  mapper->setVerbosity(PROGRESS) ; 
    49   mapper->setSourceMesh(src_bounds_lon, src_bounds_lat, src_area, n_vert_per_cell_src, n_cell_src, src_pole ) ; 
    50   mapper->setTargetMesh(dst_bounds_lon, dst_bounds_lat, dst_area, n_vert_per_cell_dst, n_cell_dst, dst_pole ) ; 
     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 ) ; 
    5150 
    5251/* 

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

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

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

    r1619 r1630  
    2323       void setVerbosity(verbosity v) {verbose=v ;} 
    2424 
    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) ; 
     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) ; 
    2727       void setSourceValue(const double* val) ; 
    2828       void getTargetValue(double* val) ; 

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

    r1619 r1630  
    22#include "elt.hpp" 
    33#include "polyg.hpp" 
    4 #include "intersection_ym.hpp" 
    5 #include "earcut.hpp" 
    6 #include <vector> 
    74 
    85namespace sphereRemap { 
    96 
    107using namespace std; 
    11  
    12 double 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  
    658 
    669void cptEltGeom(Elt& elt, const Coord &pole) 
    … …  
    7114  elt.area = airbar(elt.n, elt.vertex, elt.edge, elt.d, pole, gg); 
    7215  elt.x = gg; 
    73 // overwrite area computation  
    74  
    75   elt.area =  computePolygoneArea(elt, pole) ; 
    7616} 
    77  
    7817 
    7918void cptAllEltsGeom(Elt *elt, int N, const Coord &pole) 

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

    r1619 r1630  
    218218int packedPolygonSize(const Elt& e) 
    219219{ 
    220   return sizeof(e.id) + sizeof(e.src_id) + sizeof(e.x) + sizeof(e.val)  + sizeof(e.given_area)+ 
     220  return sizeof(e.id) + sizeof(e.src_id) + sizeof(e.x) + sizeof(e.val) + 
    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  
    240237  *((int *) & (buffer[pos])) = e.n; 
    241238  pos += sizeof(e.n); 
    … …  
    265262  pos += sizeof(double); 
    266263 
    267   e.given_area = *((double *) & (buffer[pos])); 
    268   pos += sizeof(double); 
    269  
    270264  e.n = *((int *) & (buffer[pos])); 
    271265  pos += sizeof(int); 
    … …  
    297291    *((double *) &(buffer[pos])) = e.area; 
    298292    pos += sizeof(double); 
    299  
    300293 
    301294    *((GloId *) &(buffer[pos])) = (*it)->id; 
    … …  
    329322    pos += sizeof(double); 
    330323 
    331  
    332324    Polyg *polygon = new Polyg; 
    333325 

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

    r1629 r1630  
    119119        this->timestep = timestep; 
    120120      } 
    121 /* 
     121 
    122122      int CCalendar::getStep(void) const 
    123123      { 
    124124        return step; 
    125125      } 
    126 */ 
     126 
    127127      const CDate& CCalendar::update(int step) 
    128128      { 
    … …  
    153153      //----------------------------------------------------------------- 
    154154 
    155       //const CDuration& CCalendar::getTimeStep(void) const { return this->timestep; } 
    156       //const CDate& CCalendar::getInitDate(void) const     { return this->initDate; } 
    157       //const CDate& CCalendar::getTimeOrigin(void) const   { return this->timeOrigin; } 
    158       //const CDate& CCalendar::getCurrentDate(void) const  { return this->currentDate; } 
     155      const CDuration& CCalendar::getTimeStep(void) const { return this->timestep; } 
     156      const CDate& CCalendar::getInitDate(void) const     { return this->initDate; } 
     157      const CDate& CCalendar::getTimeOrigin(void) const   { return this->timeOrigin; } 
     158      const CDate& CCalendar::getCurrentDate(void) const  { return this->currentDate; } 
    159159 
    160160      //----------------------------------------------------------------- 
    … …  
    169169      StdString CCalendar::getType(void) const { return StdString(this->getId()); } 
    170170 
    171       //int CCalendar::getYearTotalLength(const CDate& date) const { return (365 * 86400); } 
     171      int CCalendar::getYearTotalLength(const CDate& date) const { return (365 * 86400); } 
    172172 
    173173      //int CCalendar::getYearLength  (void) const { return 12; } 
    … …  
    177177      //int CCalendar::getDayLengthInSeconds(void) const { return getDayLength() * getHourLength() * getMinuteLength(); } 
    178178 
    179       //bool CCalendar::hasLeapYear() const { return false; } 
    180  
    181       /*StdString CCalendar::getMonthName(int monthId) const 
     179      bool CCalendar::hasLeapYear() const { return false; } 
     180 
     181      StdString CCalendar::getMonthName(int monthId) const 
    182182      { 
    183183        static const StdString MonthNames[] = 
    … …  
    185185            "july",    "august",   "september", "october", "november", "december" }; 
    186186        return MonthNames[monthId - 1]; 
    187       }*/ 
     187      } 
    188188 
    189189      const StdString CCalendar::getMonthShortName(int monthId) const 

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

    r1629 r1630  
    5252 
    5353            /// Mutateur /// 
    54             void setTimeStep(const CDuration& timestep) ; 
     54            void setTimeStep(const CDuration& timestep); 
    5555            void setInitDate(const CDate& initDate); 
    5656            void setTimeOrigin(const CDate& timeOrigin); 
    … …  
    6060 
    6161            /// Accesseurs /// 
    62             const CDuration& getTimeStep(void) const { return this->timestep; }; 
    63             const CDate& getInitDate(void) const { return this->initDate; }; 
    64             const CDate& getTimeOrigin(void) const { return this->timeOrigin; }; 
    65             const CDate& getCurrentDate(void) const { return this->currentDate; }; 
     62            const CDuration& getTimeStep(void) const; 
     63            const CDate& getInitDate(void) const; 
     64            const CDate& getTimeOrigin(void) const; 
     65            const CDate& getCurrentDate(void) const; 
    6666 
    6767         public : 
    … …  
    7070            virtual StdString getType(void) const; 
    7171 
    72             int getStep(void) const {return step;}; 
     72            int getStep(void) const; 
    7373 
    7474            inline int getMonthLength(const CDate& date) const 
    … …  
    7878            }; 
    7979 
    80             //virtual int getYearTotalLength(const CDate& date) const; // Retourne la durée d'une année en seconde. 
    81             inline virtual int getYearTotalLength(const CDate& date) const { return (365 * 86400); } ; // Retourne la durée d'une année en seconde. 
     80            virtual int getYearTotalLength(const CDate& date) const; // Retourne la durée d'une année en seconde. 
    8281 
    8382            //virtual int getYearLength  (void) const; // Retourne la durée d'une année en mois. 
    84             inline virtual int getYearLength (void) const { return 12; } ; 
    85             inline virtual int getDayLength   (void) const { return 24; } ; // Retourne la durée d'un jour en heures. 
    86             inline virtual int getHourLength  (void) const { return 60; } ; // Retourne la durée d'une heure en minute. 
    87             inline virtual int getMinuteLength(void) const {return 60; } ; // Retourne la durée d'une minute en secondes. 
     83            inline int getYearLength (void) const { return 12; } ; 
     84            inline int getDayLength   (void) const { return 24; } ; // Retourne la durée d'un jour en heures. 
     85            inline int getHourLength  (void) const { return 60; } ; // Retourne la durée d'une heure en minute. 
     86            inline int getMinuteLength(void) const {return 60; } ; // Retourne la durée d'une minute en secondes. 
    8887            /*! Returns the day length expressed in seconds. */ 
    89             inline virtual int getDayLengthInSeconds(void) const { return 86400; } ; 
     88            inline int getDayLengthInSeconds(void) const { return 86400; } ; 
    9089 
    91             inline virtual StdString getMonthName(int monthId) const 
    92             { 
    93               static const StdString MonthNames[] = 
    94                          { "january", "february", "march",     "april" ,  "may",      "june", 
    95                            "july",    "august",   "september", "october", "november", "december" }; 
    96               return MonthNames[monthId - 1]; 
    97             }; 
     90            virtual StdString getMonthName(int monthId) const; 
    9891            virtual const StdString getMonthShortName(int monthId) const; 
    9992 
    10093            /*! Test if the calendar can have leap year. */ 
    101             inline virtual bool hasLeapYear() const {return false;}; 
     94            virtual bool hasLeapYear() const; 
    10295 
    10396            void initializeDate(int yr, int mth, int d, int hr = 0, int min = 0, int sec = 0); 

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

    r1629 r1630  
    11#include "calendar_util.hpp" 
     2#include "calendar.hpp" 
    23 
    34namespace xios 

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

    r1619 r1630  
    5858 
    5959DECLARE_ARRAY(double, 2, area) 
    60 DECLARE_ATTRIBUTE(double, radius) 
    6160 
    6261DECLARE_ENUM4(type,rectilinear,curvilinear,unstructured, gaussian) 

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

    r1619 r1630  
    1010DECLARE_ATTRIBUTE(bool, write_weight) 
    1111DECLARE_ENUM2(read_write_convention, c, fortran) 
    12 DECLARE_ATTRIBUTE(bool, use_area) 

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

    r1619 r1630  
    9696    { 
    9797      list<int> ranks = event.getRanks(); 
    98       info(100)<<"Event "<<timeLine<<" of context "<<context->getId()<<endl ; 
     98 
    9999      if (CXios::checkEventSync) 
    100100      { 
    … …  
    124124        { 
    125125          event.send(timeLine, sizes, buffList); 
    126           info(100)<<"Event "<<timeLine<<" of context "<<context->getId()<<"  sent"<<endl ; 
    127126 
    128127          checkBuffers(ranks); 
    … …  
    143142          info(100)<<"DEBUG : temporaly event created : timeline "<<timeLine<<endl ; 
    144143          event.send(timeLine, tmpBufferedEvent.sizes, tmpBufferedEvent.buffers); 
    145           info(100)<<"Event "<<timeLine<<" of context "<<context->getId()<<"  sent"<<endl ; 
    146144        } 
    147145      } 

  • XIOS/dev/dev_trunk_omp/src/date/allleap.cpp

    r1629 r1630  
    3232 
    3333      ///-------------------------------------------------------------- 
    34 /* 
     34 
    3535      int CAllLeapCalendar::getYearTotalLength(const CDate & date) const 
    3636      { return (366 * 86400); } 
    37 */ 
     37 
    3838      int CAllLeapCalendar::getMonthLength(const CDate & date) const 
    3939      { 
    … …  
    4141         return (CCalendar::getMonthLength(date)); 
    4242      } 
    43 /* 
     43 
    4444      StdString CAllLeapCalendar::getType(void) const 
    4545      { return (StdString("all_leap")); } 
    46 */ 
     46 
    4747      ///-------------------------------------------------------------- 
    4848} // namespace xios 

  • XIOS/dev/dev_trunk_omp/src/date/allleap.hpp

    r1629 r1630  
    2626 
    2727            /// Accesseurs /// 
    28             inline virtual int getYearTotalLength(const CDate & date) const { return (366 * 86400); }; 
     28            virtual int getYearTotalLength(const CDate & date) const; 
    2929            virtual int getMonthLength(const CDate & date) const; 
    30             virtual StdString getType(void) const { return (StdString("all_leap")); }; 
     30            virtual StdString getType(void) const; 
    3131 
    3232            /// Destructeur /// 
    33             virtual ~CAllLeapCalendar(void) ; 
     33            virtual ~CAllLeapCalendar(void); 
    3434 
    3535      }; // class CAllLeapCalendar 

  • XIOS/dev/dev_trunk_omp/src/date/d360.cpp

    r1629 r1630  
    3333      ///-------------------------------------------------------------- 
    3434 
    35       //int CD360Calendar::getYearTotalLength(const CDate & date) const 
    36       //{ return (360 * 86400); } 
     35      int CD360Calendar::getYearTotalLength(const CDate & date) const 
     36      { return (360 * 86400); } 
    3737 
    38       //int CD360Calendar::getMonthLength(const CDate & date) const 
    39       //{ return (30); } 
     38      int CD360Calendar::getMonthLength(const CDate & date) const 
     39      { return (30); } 
    4040 
    41       //StdString CD360Calendar::getType(void) const 
    42       //{ return (StdString("360_day")); } 
     41      StdString CD360Calendar::getType(void) const 
     42      { return (StdString("360_day")); } 
    4343 
    4444      ///-------------------------------------------------------------- 

  • XIOS/dev/dev_trunk_omp/src/date/d360.hpp

    r1629 r1630  
    2626 
    2727            /// Accesseurs /// 
    28             inline virtual int getYearTotalLength(const CDate & date) const { return (360 * 86400); }; 
    29             inline virtual int getMonthLength(const CDate & date) const { return 30; }; 
    30             inline virtual StdString getType(void) const { return (StdString("360_day")); }; 
     28            virtual int getYearTotalLength(const CDate & date) const; 
     29            virtual int getMonthLength(const CDate & date) const; 
     30            virtual StdString getType(void) const; 
    3131 
    3232            /// Destructeur /// 

  • XIOS/dev/dev_trunk_omp/src/date/gregorian.cpp

    r1629 r1630  
    5555      } 
    5656 
    57       //StdString CGregorianCalendar::getType(void) const 
    58       //{ return (StdString("gregorian")); } 
     57      StdString CGregorianCalendar::getType(void) const 
     58      { return (StdString("gregorian")); } 
    5959 
    60       //bool CGregorianCalendar::hasLeapYear() const { return true; } 
     60      bool CGregorianCalendar::hasLeapYear() const { return true; } 
    6161 
    6262      ///-------------------------------------------------------------- 

  • XIOS/dev/dev_trunk_omp/src/date/gregorian.hpp

    r1629 r1630  
    2828            virtual int getYearTotalLength(const CDate & date) const; 
    2929            virtual int getMonthLength(const CDate & date) const; 
    30             inline virtual StdString getType(void) const { return (StdString("gregorian")); }; 
    31             inline virtual int getYearLength (void) const { return 12; } ; 
     30            virtual StdString getType(void) const; 
    3231 
    33             inline virtual bool hasLeapYear() const { return true; }; 
     32            virtual bool hasLeapYear() const; 
    3433 
    3534            /// Destructeur /// 

  • XIOS/dev/dev_trunk_omp/src/io/onetcdf4_impl.hpp

    r1619 r1630  
    7373    } 
    7474    char *PtrArrayStr ; 
    75     PtrArrayStr=new char[stringArrayLen*data.numElements()] ; 
    76     memset (PtrArrayStr,' ',stringArrayLen*data.numElements()); 
    77     size_t offset=0 ; 
     75    PtrArrayStr=new char[stringArrayLen] ; 
    7876    Array<StdString,1>::const_iterator it, itb=data.begin(), ite=data.end() ; 
    79     for(it=itb;it!=ite;++it, offset+=stringArrayLen) 
     77    int lineNb = 0; 
     78    for(it=itb;it!=ite;++it) 
    8079    { 
    81       it->copy(PtrArrayStr+offset,it->size()) ; 
    82       PtrArrayStr[offset+it->size()]='\0' ; 
     80      it->copy(PtrArrayStr,it->size()) ; 
     81      PtrArrayStr[it->size()]='\0' ; 
     82      sstart[0] = lineNb; 
     83      sstart[dimArrayLen] = 0; 
     84      scount[0] = 1; 
     85      scount[dimArrayLen] = it->size() + 1; 
     86      CTimer::get("CONetCDF4::writeData writeData_").resume(); 
     87      this->writeData_(grpid, varid, sstart, scount, PtrArrayStr); 
     88      CTimer::get("CONetCDF4::writeData writeData_").suspend(); 
     89      ++lineNb; 
    8390    } 
    84  
    85      CTimer::get("CONetCDF4::writeData writeData_").resume(); 
    86      this->writeData_(grpid, varid, sstart, scount, PtrArrayStr); 
    87      CTimer::get("CONetCDF4::writeData writeData_").suspend(); 
    88  
    8991    delete []  PtrArrayStr; 
    9092  } 

  • XIOS/dev/dev_trunk_omp/src/node/interpolate_domain.cpp

    r1619 r1630  
    6565    if (this->read_write_convention.isEmpty()) this->read_write_convention.setValue(read_write_convention_attr::fortran); 
    6666 
    67  
    6867  } 
    6968 

  • XIOS/dev/dev_trunk_omp/src/transformation/domain_algorithm_interpolate.cpp

    r1619 r1630  
    305305  CArray<double,2> boundsLonSrcUnmasked(nVertexSrc,nSrcLocalUnmasked); 
    306306  CArray<double,2> boundsLatSrcUnmasked(nVertexSrc,nSrcLocalUnmasked); 
    307   CArray<double,1> areaSrcUnmasked(nSrcLocalUnmasked); 
    308    
    309307  long int * globalSrcUnmasked = new long int [nSrcLocalUnmasked]; 
    310308 
    311309  nSrcLocalUnmasked=0 ; 
    312   bool hasSrcArea=domainSrc_->hasArea && !domainSrc_->radius.isEmpty() && !interpDomain_->use_area.isEmpty() && interpDomain_->use_area==true  ; 
    313   double srcAreaFactor ; 
    314   if (hasSrcArea) srcAreaFactor=1./(domainSrc_->radius*domainSrc_->radius) ; 
    315    
    316310  for (int idx=0 ; idx < nSrcLocal; idx++) 
    317311  { 
    … …  
    323317        boundsLatSrcUnmasked(n,nSrcLocalUnmasked) = boundsLatSrc(n,idx) ; 
    324318      } 
    325       if (hasSrcArea) areaSrcUnmasked(nSrcLocalUnmasked) = domainSrc_->areavalue(idx)*srcAreaFactor ; 
    326319      globalSrcUnmasked[nSrcLocalUnmasked]=globalSrc[idx] ; 
    327320      ++nSrcLocalUnmasked ; 
    328321    } 
    329322  } 
    330   
     323 
    331324 
    332325  int nDstLocalUnmasked = 0 ; 
    … …  
    335328  CArray<double,2> boundsLonDestUnmasked(nVertexDest,nDstLocalUnmasked); 
    336329  CArray<double,2> boundsLatDestUnmasked(nVertexDest,nDstLocalUnmasked); 
    337   CArray<double,1>   areaDstUnmasked(nDstLocalUnmasked); 
    338  
    339330  long int * globalDstUnmasked = new long int [nDstLocalUnmasked]; 
    340331 
    341332  nDstLocalUnmasked=0 ; 
    342   bool hasDstArea=domainDest_->hasArea && !domainDest_->radius.isEmpty() && !interpDomain_->use_area.isEmpty() && interpDomain_->use_area==true ; 
    343   double dstAreaFactor ; 
    344   if (hasDstArea) dstAreaFactor=1./(domainDest_->radius*domainDest_->radius) ; 
    345333  for (int idx=0 ; idx < nDstLocal; idx++) 
    346334  { 
    … …  
    352340        boundsLatDestUnmasked(n,nDstLocalUnmasked) = boundsLatDest(n,idx) ; 
    353341      } 
    354       if (hasDstArea) areaDstUnmasked(nDstLocalUnmasked) = domainDest_->areavalue(idx)*dstAreaFactor ; 
    355342      globalDstUnmasked[nDstLocalUnmasked]=globalDst[idx] ; 
    356343      ++nDstLocalUnmasked ; 
    … …  
    358345  } 
    359346 
    360   double* ptAreaSrcUnmasked = NULL ; 
    361   if (hasSrcArea) ptAreaSrcUnmasked=areaSrcUnmasked.dataFirst() ; 
    362  
    363   double* ptAreaDstUnmasked = NULL ; 
    364   if (hasDstArea) ptAreaDstUnmasked=areaDstUnmasked.dataFirst() ; 
    365  
    366   mapper.setSourceMesh(boundsLonSrcUnmasked.dataFirst(), boundsLatSrcUnmasked.dataFirst(), ptAreaSrcUnmasked, nVertexSrc, nSrcLocalUnmasked, &srcPole[0], globalSrcUnmasked); 
    367   mapper.setTargetMesh(boundsLonDestUnmasked.dataFirst(), boundsLatDestUnmasked.dataFirst(), ptAreaDstUnmasked, nVertexDest, nDstLocalUnmasked, &dstPole[0], globalDstUnmasked); 
     347  mapper.setSourceMesh(boundsLonSrcUnmasked.dataFirst(), boundsLatSrcUnmasked.dataFirst(), nVertexSrc, nSrcLocalUnmasked, &srcPole[0], globalSrcUnmasked); 
     348  mapper.setTargetMesh(boundsLonDestUnmasked.dataFirst(), boundsLatDestUnmasked.dataFirst(), nVertexDest, nDstLocalUnmasked, &dstPole[0], globalDstUnmasked); 
    368349 
    369350  std::vector<double> timings = mapper.computeWeights(orderInterp,renormalize,quantity); 

  • XIOS/dev/dev_trunk_omp/src/transformation/domain_algorithm_zoom.cpp

    r1611 r1630  
    9393  domainDest_->ni.setValue(niDest); 
    9494  domainDest_->nj.setValue(njDest); 
    95   if ( (niDest==0) || (njDest==0)) 
    96   { 
    97     domainDest_->ibegin.setValue(0); 
    98     domainDest_->jbegin.setValue(0); 
    99   } 
    100   else 
    101   { 
    102     domainDest_->ibegin.setValue(ibeginDest); 
    103     domainDest_->jbegin.setValue(jbeginDest); 
    104   } 
     95  domainDest_->ibegin.setValue(ibeginDest); 
     96  domainDest_->jbegin.setValue(jbeginDest); 
    10597  domainDest_->i_index.resize(niDest*njDest); 
    10698  domainDest_->j_index.resize(niDest*njDest); 
Note: See TracChangeset for help on using the changeset viewer.