src/ParaMEDMEM/ParaUMesh.cxx

   1 //
   2 // Copyright (C) 2020  CEA/DEN, EDF R&D
   3 //
   4 // This library is free software; you can redistribute it and/or
   5 // modify it under the terms of the GNU Lesser General Public
   6 // License as published by the Free Software Foundation; either
   7 // version 2.1 of the License, or (at your option) any later version.
   8 //
   9 // This library is distributed in the hope that it will be useful,
  10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12 // Lesser General Public License for more details.
  13 //
  14 // You should have received a copy of the GNU Lesser General Public
  15 // License along with this library; if not, write to the Free Software
  16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  17 //
  18 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  19 //
  20 // Author : Anthony Geay (EDF R&D)
  21
  22 #include "ParaUMesh.hxx"
  23 #include "ProcessorGroup.hxx"
  24 #include "MPIProcessorGroup.hxx"
  25 #include "Topology.hxx"
  26 #include "BlockTopology.hxx"
  27 #include "CommInterface.hxx"
  28 #include "MEDCouplingMemArray.hxx"
  29
  30 #include "mpi.h"
  31
  32 #include <fstream>
  33 #include <sstream>
  34 #include <numeric>
  35 #include <memory>
  36 #include <vector>
  37
  38 using namespace MEDCoupling;
  39
  40 ParaUMesh *ParaUMesh::New(MEDCouplingUMesh *mesh, DataArrayIdType *globalCellIds, DataArrayIdType *globalNodeIds)
  41 {
  42   return new ParaUMesh(mesh,globalCellIds,globalNodeIds);
  43 }
  44
  45 ParaUMesh::ParaUMesh(MEDCouplingUMesh *mesh, DataArrayIdType *globalCellIds, DataArrayIdType *globalNodeIds)
  46 {
  47   _mesh.takeRef(mesh);
  48   _cell_global.takeRef(globalCellIds);
  49   _node_global.takeRef(globalNodeIds);
  50   _mesh.checkNotNull();
  51   _cell_global.checkNotNull();
  52   _node_global.checkNotNull();
  53   _mesh->checkConsistencyLight();
  54   if(_mesh->getNumberOfNodes() != _node_global->getNumberOfTuples())
  55     throw INTERP_KERNEL::Exception("ParaUMesh constructor : mismatch between # nodes and len of global # nodes.");
  56   if(_mesh->getNumberOfCells() != _cell_global->getNumberOfTuples())
  57     throw INTERP_KERNEL::Exception("ParaUMesh constructor : mismatch between # cells and len of global # cells.");
  58 }
  59
  60 std::size_t ParaUMesh::getHeapMemorySizeWithoutChildren() const
  61 {
  62   return 0;
  63 }
  64
  65 std::vector<const BigMemoryObject *> ParaUMesh::getDirectChildrenWithNull() const
  66 {
  67   return {_mesh,_cell_global,_node_global};
  68 }
  69
  70 /*!
  71 * This method computes the cells part of distributed mesh lying on \a globalNodeIds nodes.
  72 * The input \a globalNodeIds are not supposed to reside on the current process.
  73 */
  74 MCAuto<DataArrayIdType> ParaUMesh::getCellIdsLyingOnNodes(const DataArrayIdType *globalNodeIds, bool fullyIn) const
  75 {
  76   if(fullyIn)
  77     throw INTERP_KERNEL::Exception("ParaUMesh::getCellIdsLyingOnNodes : not implemented yet for fullyIn == True !");
  78   MPI_Comm comm(MPI_COMM_WORLD);
  79   CommInterface ci;
  80   int size;
  81   ci.commSize(comm,&size);
  82   std::unique_ptr<mcIdType[]> nbOfElems(new mcIdType[size]),nbOfElems2(new mcIdType[size]),nbOfElems3(new mcIdType[size]);
  83   mcIdType nbOfNodeIdsLoc(globalNodeIds->getNumberOfTuples());
  84   ci.allGather(&nbOfNodeIdsLoc,1,MPI_ID_TYPE,nbOfElems.get(),1,MPI_ID_TYPE,comm);
  85   std::vector< MCAuto<DataArrayIdType> > tabs(size);
  86   // loop to avoid to all procs to have all the nodes per proc
  87   for(int subDiv = 0 ; subDiv < size ; ++subDiv)
  88   {
  89     std::unique_ptr<mcIdType[]> nbOfElemsSp(CommInterface::SplitArrayOfLength(nbOfElems,size,subDiv,size));
  90     mcIdType nbOfNodeIdsSum(std::accumulate(nbOfElemsSp.get(),nbOfElemsSp.get()+size,0));
  91     std::unique_ptr<mcIdType[]> allGlobalNodeIds(new mcIdType[nbOfNodeIdsSum]);
  92     std::unique_ptr<int[]> nbOfElemsInt( CommInterface::ToIntArray<mcIdType>(nbOfElemsSp,size) );
  93     std::unique_ptr<int[]> offsetsIn( CommInterface::ComputeOffset(nbOfElemsInt,size) );
  94     mcIdType startGlobalNodeIds,endGlobalNodeIds;
  95     DataArray::GetSlice(0,globalNodeIds->getNumberOfTuples(),1,subDiv,size,startGlobalNodeIds,endGlobalNodeIds);
  96     ci.allGatherV(globalNodeIds->begin()+startGlobalNodeIds,endGlobalNodeIds-startGlobalNodeIds,MPI_ID_TYPE,allGlobalNodeIds.get(),nbOfElemsInt.get(),offsetsIn.get(),MPI_ID_TYPE,comm);
  97     mcIdType offset(0);
  98     for(int curRk = 0 ; curRk < size ; ++curRk)
  99     {
 100       MCAuto<DataArrayIdType> globalNodeIdsOfCurProc(DataArrayIdType::New());
 101       globalNodeIdsOfCurProc->useArray(allGlobalNodeIds.get()+offset,false,DeallocType::CPP_DEALLOC,nbOfElemsSp[curRk],1);
 102       offset += nbOfElemsSp[curRk];
 103       MCAuto<DataArrayIdType> globalNodeIdsCaptured(_node_global->buildIntersection(globalNodeIdsOfCurProc));
 104       MCAuto<DataArrayIdType> localNodeIdsToLocate(_node_global->findIdForEach(globalNodeIdsCaptured->begin(),globalNodeIdsCaptured->end()));
 105       MCAuto<DataArrayIdType> localCellCaptured(_mesh->getCellIdsLyingOnNodes(localNodeIdsToLocate->begin(),localNodeIdsToLocate->end(),false));
 106       MCAuto<DataArrayIdType> localCellCapturedGlob(_cell_global->selectByTupleIdSafe(localCellCaptured->begin(),localCellCaptured->end()));
 107       if(tabs[curRk].isNull())
 108         tabs[curRk] = localCellCapturedGlob;
 109       else
 110         tabs[curRk]->insertAtTheEnd(localCellCapturedGlob->begin(),localCellCapturedGlob->end());
 111     }
 112   }
 113   for(int curRk = 0 ; curRk < size ; ++curRk)
 114   {
 115     tabs[curRk] = tabs[curRk]->buildUniqueNotSorted();
 116     nbOfElems3[curRk] = tabs[curRk]->getNumberOfTuples();
 117   }
 118   std::vector<const DataArrayIdType *> tabss(tabs.begin(),tabs.end());
 119   MCAuto<DataArrayIdType> cells(DataArrayIdType::Aggregate(tabss));
 120   ci.allToAll(nbOfElems3.get(),1,MPI_ID_TYPE,nbOfElems2.get(),1,MPI_ID_TYPE,comm);
 121   mcIdType nbOfCellIdsSum(std::accumulate(nbOfElems2.get(),nbOfElems2.get()+size,0));
 122   MCAuto<DataArrayIdType> cellIdsFromProcs(DataArrayIdType::New());
 123   cellIdsFromProcs->alloc(nbOfCellIdsSum,1);
 124   {
 125     std::unique_ptr<int[]> nbOfElemsInt( CommInterface::ToIntArray<mcIdType>(nbOfElems3,size) ),nbOfElemsOutInt( CommInterface::ToIntArray<mcIdType>(nbOfElems2,size) );
 126     std::unique_ptr<int[]> offsetsIn( CommInterface::ComputeOffset(nbOfElemsInt,size) ), offsetsOut( CommInterface::ComputeOffset(nbOfElemsOutInt,size) );
 127     ci.allToAllV(cells->begin(),nbOfElemsInt.get(),offsetsIn.get(),MPI_ID_TYPE,
 128                  cellIdsFromProcs->getPointer(),nbOfElemsOutInt.get(),offsetsOut.get(),MPI_ID_TYPE,comm);
 129   }
 130   cellIdsFromProcs->sort();
 131   return cellIdsFromProcs;
 132 }
 133
 134 DataArrayIdType *ParaUMesh::redistributeCellField(const DataArrayIdType *globalCellIds, const DataArrayIdType *fieldValueToRed) const
 135 {
 136   MPI_Comm comm(MPI_COMM_WORLD);
 137   CommInterface ci;
 138   std::unique_ptr<mcIdType[]> allGlobalCellIds,allGlobalCellIdsIndex;
 139   int size(ci.allGatherArrays(comm,globalCellIds,allGlobalCellIds,allGlobalCellIdsIndex));
 140   // Prepare ParaUMesh parts to be sent : compute for each proc the contribution of current rank.
 141   std::vector< MCAuto<DataArrayIdType> > globalCellIdsToBeSent(size);
 142   std::vector< MCAuto<DataArrayIdType> > fieldToBeSent(size);
 143   for(int curRk = 0 ; curRk < size ; ++curRk)
 144   {
 145     mcIdType offset(allGlobalCellIdsIndex[curRk]);
 146     MCAuto<DataArrayIdType> globalCellIdsOfCurProc(DataArrayIdType::New());
 147     globalCellIdsOfCurProc->useArray(allGlobalCellIds.get()+offset,false,DeallocType::CPP_DEALLOC,allGlobalCellIdsIndex[curRk+1]-offset,1);
 148     // the key call is here : compute for rank curRk the cells to be sent
 149     MCAuto<DataArrayIdType> globalCellIdsCaptured(_cell_global->buildIntersection(globalCellIdsOfCurProc));// OK for the global cellIds
 150     MCAuto<DataArrayIdType> localCellIdsCaptured(_cell_global->findIdForEach(globalCellIdsCaptured->begin(),globalCellIdsCaptured->end()));
 151     globalCellIdsToBeSent[curRk] = globalCellIdsCaptured;
 152     fieldToBeSent[curRk] = fieldValueToRed->selectByTupleIdSafe(localCellIdsCaptured->begin(),localCellIdsCaptured->end());
 153   }
 154   // Receive
 155   std::vector< MCAuto<DataArrayIdType> > globalCellIdsReceived;
 156   ci.allToAllArrays(comm,globalCellIdsToBeSent,globalCellIdsReceived);
 157   std::vector< MCAuto<DataArrayIdType> > fieldValueReceived;
 158   ci.allToAllArrays(comm,fieldToBeSent,fieldValueReceived);
 159   // use globalCellIdsReceived to reorganize everything
 160   MCAuto<DataArrayIdType> aggregatedCellIds( DataArrayIdType::Aggregate(FromVecAutoToVecOfConst<DataArrayIdType>(globalCellIdsReceived)) );
 161   MCAuto<DataArrayIdType> aggregatedCellIdsSorted(aggregatedCellIds->copySorted());
 162   MCAuto<DataArrayIdType> idsIntoAggregatedIds(DataArrayIdType::FindPermutationFromFirstToSecondDuplicate(aggregatedCellIdsSorted,aggregatedCellIds));
 163   MCAuto<DataArrayIdType> cellIdsOfSameNodeIds(aggregatedCellIdsSorted->indexOfSameConsecutiveValueGroups());
 164   MCAuto<DataArrayIdType> n2o_cells(idsIntoAggregatedIds->selectByTupleIdSafe(cellIdsOfSameNodeIds->begin(),cellIdsOfSameNodeIds->end()-1));//new == new ordering so that global cell ids are sorted . old == coarse ordering implied by the aggregation
 165   //
 166   MCAuto<DataArrayIdType> fieldAggregated(DataArrayIdType::Aggregate(FromVecAutoToVecOfConst<DataArrayIdType>(fieldValueReceived)));
 167   MCAuto<DataArrayIdType> ret(fieldAggregated->selectByTupleIdSafe(n2o_cells->begin(),n2o_cells->end()));
 168   return ret.retn();
 169 }
 170
 171 DataArrayIdType *ParaUMesh::redistributeNodeField(const DataArrayIdType *globalCellIds, const DataArrayIdType *fieldValueToRed) const
 172 {
 173   MPI_Comm comm(MPI_COMM_WORLD);
 174   CommInterface ci;
 175   std::unique_ptr<mcIdType[]> allGlobalCellIds,allGlobalCellIdsIndex;
 176   int size(ci.allGatherArrays(comm,globalCellIds,allGlobalCellIds,allGlobalCellIdsIndex));
 177   // Prepare ParaUMesh parts to be sent : compute for each proc the contribution of current rank.
 178   std::vector< MCAuto<DataArrayIdType> > globalNodeIdsToBeSent(size);
 179   std::vector< MCAuto<DataArrayIdType> > fieldToBeSent(size);
 180   for(int curRk = 0 ; curRk < size ; ++curRk)
 181   {
 182     mcIdType offset(allGlobalCellIdsIndex[curRk]);
 183     MCAuto<DataArrayIdType> globalCellIdsOfCurProc(DataArrayIdType::New());
 184     globalCellIdsOfCurProc->useArray(allGlobalCellIds.get()+offset,false,DeallocType::CPP_DEALLOC,allGlobalCellIdsIndex[curRk+1]-offset,1);
 185     // the key call is here : compute for rank curRk the cells to be sent
 186     MCAuto<DataArrayIdType> globalCellIdsCaptured(_cell_global->buildIntersection(globalCellIdsOfCurProc));// OK for the global cellIds
 187     MCAuto<DataArrayIdType> localCellIdsCaptured(_cell_global->findIdForEach(globalCellIdsCaptured->begin(),globalCellIdsCaptured->end()));
 188     MCAuto<MEDCouplingUMesh> meshPart(_mesh->buildPartOfMySelf(localCellIdsCaptured->begin(),localCellIdsCaptured->end(),true));
 189     MCAuto<DataArrayIdType> o2n(meshPart->zipCoordsTraducer());// OK for the mesh
 190     MCAuto<DataArrayIdType> n2o(o2n->invertArrayO2N2N2O(meshPart->getNumberOfNodes()));
 191     MCAuto<DataArrayIdType> globalNodeIdsPart(_node_global->selectByTupleIdSafe(n2o->begin(),n2o->end())); // OK for the global nodeIds
 192     globalNodeIdsToBeSent[curRk] = globalNodeIdsPart;
 193     fieldToBeSent[curRk] = fieldValueToRed->selectByTupleIdSafe(n2o->begin(),n2o->end());
 194   }
 195   // Receive
 196   std::vector< MCAuto<DataArrayIdType> > globalNodeIdsReceived;
 197   ci.allToAllArrays(comm,globalNodeIdsToBeSent,globalNodeIdsReceived);
 198   std::vector< MCAuto<DataArrayIdType> > fieldValueReceived;
 199   ci.allToAllArrays(comm,fieldToBeSent,fieldValueReceived);
 200   // firstly deal with nodes.
 201   MCAuto<DataArrayIdType> aggregatedNodeIds( DataArrayIdType::Aggregate(FromVecAutoToVecOfConst<DataArrayIdType>(globalNodeIdsReceived)) );
 202   MCAuto<DataArrayIdType> aggregatedNodeIdsSorted(aggregatedNodeIds->copySorted());
 203   MCAuto<DataArrayIdType> nodeIdsIntoAggregatedIds(DataArrayIdType::FindPermutationFromFirstToSecondDuplicate(aggregatedNodeIdsSorted,aggregatedNodeIds));
 204   MCAuto<DataArrayIdType> idxOfSameNodeIds(aggregatedNodeIdsSorted->indexOfSameConsecutiveValueGroups());
 205   MCAuto<DataArrayIdType> n2o_nodes(nodeIdsIntoAggregatedIds->selectByTupleIdSafe(idxOfSameNodeIds->begin(),idxOfSameNodeIds->end()-1));//new == new ordering so that global node ids are sorted . old == coarse ordering implied by the aggregation
 206   //
 207   MCAuto<DataArrayIdType> fieldAggregated(DataArrayIdType::Aggregate(FromVecAutoToVecOfConst<DataArrayIdType>(fieldValueReceived)));
 208   MCAuto<DataArrayIdType> ret(fieldAggregated->selectByTupleIdSafe(n2o_nodes->begin(),n2o_nodes->end()));
 209   //
 210   return ret.retn();
 211 }
 212
 213 /*!
 214  * Return part of \a this mesh split over COMM_WORLD. Part is defined by global cell ids array \a globaCellIds.
 215  */
 216 ParaUMesh *ParaUMesh::redistributeCells(const DataArrayIdType *globalCellIds) const
 217 {
 218   MPI_Comm comm(MPI_COMM_WORLD);
 219   CommInterface ci;
 220   std::unique_ptr<mcIdType[]> allGlobalCellIds,allGlobalCellIdsIndex;
 221   int size(ci.allGatherArrays(comm,globalCellIds,allGlobalCellIds,allGlobalCellIdsIndex));
 222   // Prepare ParaUMesh parts to be sent : compute for each proc the contribution of current rank.
 223   std::vector< MCAuto<DataArrayIdType> > globalCellIdsToBeSent(size),globalNodeIdsToBeSent(size);
 224   std::vector< MCAuto<MEDCouplingUMesh> > meshPartsToBeSent(size);
 225   for(int curRk = 0 ; curRk < size ; ++curRk)
 226   {
 227     mcIdType offset(allGlobalCellIdsIndex[curRk]);
 228     MCAuto<DataArrayIdType> globalCellIdsOfCurProc(DataArrayIdType::New());
 229     globalCellIdsOfCurProc->useArray(allGlobalCellIds.get()+offset,false,DeallocType::CPP_DEALLOC,allGlobalCellIdsIndex[curRk+1]-offset,1);
 230     // the key call is here : compute for rank curRk the cells to be sent
 231     MCAuto<DataArrayIdType> globalCellIdsCaptured(_cell_global->buildIntersection(globalCellIdsOfCurProc));// OK for the global cellIds
 232     MCAuto<DataArrayIdType> localCellIdsCaptured(_cell_global->findIdForEach(globalCellIdsCaptured->begin(),globalCellIdsCaptured->end()));
 233     MCAuto<MEDCouplingUMesh> meshPart(_mesh->buildPartOfMySelf(localCellIdsCaptured->begin(),localCellIdsCaptured->end(),true));
 234     MCAuto<DataArrayIdType> o2n(meshPart->zipCoordsTraducer());// OK for the mesh
 235     MCAuto<DataArrayIdType> n2o(o2n->invertArrayO2N2N2O(meshPart->getNumberOfNodes()));
 236     MCAuto<DataArrayIdType> globalNodeIdsPart(_node_global->selectByTupleIdSafe(n2o->begin(),n2o->end())); // OK for the global nodeIds
 237     meshPartsToBeSent[curRk] = meshPart;
 238     globalCellIdsToBeSent[curRk] = globalCellIdsCaptured;
 239     globalNodeIdsToBeSent[curRk] = globalNodeIdsPart;
 240   }
 241   // Receive
 242   std::vector< MCAuto<DataArrayIdType> > globalCellIdsReceived,globalNodeIdsReceived;
 243   ci.allToAllArrays(comm,globalCellIdsToBeSent,globalCellIdsReceived);
 244   ci.allToAllArrays(comm,globalNodeIdsToBeSent,globalNodeIdsReceived);
 245   //now exchange the 3 arrays for the umesh : connectivity, connectivityindex and coordinates
 246   std::vector<const MEDCouplingUMesh *> meshPartsToBeSent2(FromVecAutoToVecOfConst<MEDCouplingUMesh>(meshPartsToBeSent));
 247   //connectivityindex
 248   std::vector< MCAuto<DataArrayIdType> > connectivityIndexReceived,connectivityReceived;
 249   {
 250     std::vector<const DataArrayIdType *> connectivityIndexToBeSent(UMeshConnectivityIndexIterator(0,&meshPartsToBeSent2),UMeshConnectivityIndexIterator(meshPartsToBeSent2.size(),&meshPartsToBeSent2));
 251     ci.allToAllArrays(comm,FromVecConstToVecAuto<DataArrayIdType>(connectivityIndexToBeSent),connectivityIndexReceived);
 252   }
 253   //connectivity
 254   {
 255     std::vector<const DataArrayIdType *> connectivityToBeSent(UMeshConnectivityIterator(0,&meshPartsToBeSent2),UMeshConnectivityIterator(meshPartsToBeSent2.size(),&meshPartsToBeSent2));
 256     ci.allToAllArrays(comm,FromVecConstToVecAuto<DataArrayIdType>(connectivityToBeSent),connectivityReceived);
 257   }
 258   //coordinates
 259   MCAuto<DataArrayDouble> coords;
 260   {
 261     std::vector<const DataArrayDouble *> coordsToBeSent(UMeshCoordsIterator(0,&meshPartsToBeSent2),UMeshCoordsIterator(meshPartsToBeSent2.size(),&meshPartsToBeSent2));
 262     ci.allToAllArrays(comm,FromVecConstToVecAuto<DataArrayDouble>(coordsToBeSent),coords);
 263   }
 264   /////// Sort it all !
 265   // firstly deal with nodes.
 266   MCAuto<DataArrayIdType> aggregatedNodeIds( DataArrayIdType::Aggregate(FromVecAutoToVecOfConst<DataArrayIdType>(globalNodeIdsReceived)) );
 267   MCAuto<DataArrayIdType> aggregatedNodeIdsSorted(aggregatedNodeIds->copySorted());
 268   MCAuto<DataArrayIdType> nodeIdsIntoAggregatedIds(DataArrayIdType::FindPermutationFromFirstToSecondDuplicate(aggregatedNodeIdsSorted,aggregatedNodeIds));
 269   MCAuto<DataArrayIdType> idxOfSameNodeIds(aggregatedNodeIdsSorted->indexOfSameConsecutiveValueGroups());
 270   MCAuto<DataArrayIdType> n2o_nodes(nodeIdsIntoAggregatedIds->selectByTupleIdSafe(idxOfSameNodeIds->begin(),idxOfSameNodeIds->end()-1));//new == new ordering so that global node ids are sorted . old == coarse ordering implied by the aggregation
 271   MCAuto<DataArrayIdType> finalGlobalNodeIds(aggregatedNodeIdsSorted->selectByTupleIdSafe(idxOfSameNodeIds->begin(),idxOfSameNodeIds->end()-1));
 272   MCAuto<DataArrayDouble> finalCoords(coords->selectByTupleIdSafe(n2o_nodes->begin(),n2o_nodes->end()));
 273   finalCoords->copyStringInfoFrom(*_mesh->getCoords());
 274   // secondly renumbering of node ids in connectivityReceived
 275   for(int curRk = 0 ; curRk < size ; ++curRk)
 276   {
 277     auto current(globalNodeIdsReceived[curRk]);
 278     MCAuto<DataArrayIdType> aa(finalGlobalNodeIds->findIdForEach(current->begin(),current->end()));
 279     // work on connectivityReceived[curRk] with transformWithIndArr but do not forget type of cells that should be excluded !
 280     auto connectivityToModify(connectivityReceived[curRk]);
 281     auto connectivityIndex(connectivityIndexReceived[curRk]);
 282     MCAuto<DataArrayIdType> types(connectivityToModify->selectByTupleIdSafe(connectivityIndex->begin(),connectivityIndex->end()-1));
 283     connectivityToModify->setPartOfValuesSimple3(0,connectivityIndex->begin(),connectivityIndex->end()-1,0,1,1);
 284     connectivityToModify->transformWithIndArr(aa->begin(),aa->end());
 285     connectivityToModify->setPartOfValues3(types,connectivityIndex->begin(),connectivityIndex->end()-1,0,1,1,true);
 286   }
 287   // thirdly renumber cells
 288   MCAuto<DataArrayIdType> aggregatedCellIds( DataArrayIdType::Aggregate(FromVecAutoToVecOfConst<DataArrayIdType>(globalCellIdsReceived)) );
 289   MCAuto<DataArrayIdType> aggregatedCellIdsSorted(aggregatedCellIds->copySorted());
 290   MCAuto<DataArrayIdType> idsIntoAggregatedIds(DataArrayIdType::FindPermutationFromFirstToSecondDuplicate(aggregatedCellIdsSorted,aggregatedCellIds));
 291   MCAuto<DataArrayIdType> cellIdsOfSameNodeIds(aggregatedCellIdsSorted->indexOfSameConsecutiveValueGroups());
 292   MCAuto<DataArrayIdType> n2o_cells(idsIntoAggregatedIds->selectByTupleIdSafe(cellIdsOfSameNodeIds->begin(),cellIdsOfSameNodeIds->end()-1));//new == new ordering so that global cell ids are sorted . old == coarse ordering implied by the aggregation
 293   // TODO : check coordsReceived==globalCellIds
 294   MCAuto<DataArrayIdType> connSorted,indicesSorted;
 295   {
 296     MCAuto<DataArrayIdType> conn(DataArrayIdType::Aggregate(FromVecAutoToVecOfConst<DataArrayIdType>(connectivityReceived)));
 297     MCAuto<DataArrayIdType> connIndex(DataArrayIdType::AggregateIndexes(FromVecAutoToVecOfConst<DataArrayIdType>(connectivityIndexReceived)));
 298     {
 299       DataArrayIdType *indicesSortedTmp(nullptr),*valuesSortedTmp(nullptr);
 300       DataArrayIdType::ExtractFromIndexedArrays(n2o_cells->begin(),n2o_cells->end(),conn,connIndex,valuesSortedTmp,indicesSortedTmp);
 301       indicesSorted = indicesSortedTmp; connSorted=valuesSortedTmp;
 302     }
 303   }
 304   // finalize all
 305   MCAuto<MEDCouplingUMesh> mesh(MEDCouplingUMesh::New(_mesh->getName(),_mesh->getMeshDimension()));
 306   mesh->setConnectivity(connSorted,indicesSorted,true);
 307   mesh->setCoords(finalCoords);
 308   mesh->setDescription(_mesh->getDescription());
 309   MCAuto<ParaUMesh> ret(ParaUMesh::New(mesh,aggregatedCellIdsSorted,finalGlobalNodeIds));
 310   return ret.retn();
 311 }