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(MEDCouplingUMesh *mesh, DataArrayIdType *globalCellIds, DataArrayIdType *globalNodeIds)
  41 {
  42   _mesh.takeRef(mesh);
  43   _cell_global.takeRef(globalCellIds);
  44   _node_global.takeRef(globalNodeIds);
  45   _mesh.checkNotNull();
  46   _cell_global.checkNotNull();
  47   _node_global.checkNotNull();
  48   _mesh->checkConsistencyLight();
  49   if(_mesh->getNumberOfNodes() != _node_global->getNumberOfTuples())
  50     throw INTERP_KERNEL::Exception("ParaUMesh constructor : mismatch between # nodes and len of global # nodes.");
  51   if(_mesh->getNumberOfCells() != _cell_global->getNumberOfTuples())
  52     throw INTERP_KERNEL::Exception("ParaUMesh constructor : mismatch between # cells and len of global # cells.");
  53 }
  54
  55 /*!
  56 * This method computes the cells part of distributed mesh lying on \a globalNodeIds nodes.
  57 * The input \a globalNodeIds are not supposed to reside on the current process.
  58 */
  59 MCAuto<DataArrayIdType> ParaUMesh::getCellIdsLyingOnNodes(const DataArrayIdType *globalNodeIds, bool fullyIn) const
  60 {
  61   if(fullyIn)
  62     throw INTERP_KERNEL::Exception("ParaUMesh::getCellIdsLyingOnNodes : not implemented yet for fullyIn == True !");
  63   MPI_Comm comm(MPI_COMM_WORLD);
  64   CommInterface ci;
  65   int size;
  66   ci.commSize(comm,&size);
  67   std::unique_ptr<mcIdType[]> nbOfElems(new mcIdType[size]),nbOfElems2(new mcIdType[size]),nbOfElems3(new mcIdType[size]);
  68   mcIdType nbOfNodeIdsLoc(globalNodeIds->getNumberOfTuples());
  69   ci.allGather(&nbOfNodeIdsLoc,1,MPI_ID_TYPE,nbOfElems.get(),1,MPI_ID_TYPE,comm);
  70   std::vector< MCAuto<DataArrayIdType> > tabs(size);
  71   // loop to avoid to all procs to have all the nodes per proc
  72   for(int subDiv = 0 ; subDiv < size ; ++subDiv)
  73   {
  74     std::unique_ptr<mcIdType[]> nbOfElemsSp(CommInterface::SplitArrayOfLength(nbOfElems,size,subDiv,size));
  75     mcIdType nbOfNodeIdsSum(std::accumulate(nbOfElemsSp.get(),nbOfElemsSp.get()+size,0));
  76     std::unique_ptr<mcIdType[]> allGlobalNodeIds(new mcIdType[nbOfNodeIdsSum]);
  77     std::unique_ptr<int[]> nbOfElemsInt( CommInterface::ToIntArray<mcIdType>(nbOfElemsSp,size) );
  78     std::unique_ptr<int[]> offsetsIn( CommInterface::ComputeOffset(nbOfElemsInt,size) );
  79     mcIdType startGlobalNodeIds,endGlobalNodeIds;
  80     DataArray::GetSlice(0,globalNodeIds->getNumberOfTuples(),1,subDiv,size,startGlobalNodeIds,endGlobalNodeIds);
  81     ci.allGatherV(globalNodeIds->begin()+startGlobalNodeIds,endGlobalNodeIds-startGlobalNodeIds,MPI_ID_TYPE,allGlobalNodeIds.get(),nbOfElemsInt.get(),offsetsIn.get(),MPI_ID_TYPE,comm);
  82     mcIdType offset(0);
  83     for(int curRk = 0 ; curRk < size ; ++curRk)
  84     {
  85       MCAuto<DataArrayIdType> globalNodeIdsOfCurProc(DataArrayIdType::New());
  86       globalNodeIdsOfCurProc->useArray(allGlobalNodeIds.get()+offset,false,DeallocType::CPP_DEALLOC,nbOfElemsSp[curRk],1);
  87       offset += nbOfElemsSp[curRk];
  88       MCAuto<DataArrayIdType> globalNodeIdsCaptured(_node_global->buildIntersection(globalNodeIdsOfCurProc));
  89       MCAuto<DataArrayIdType> localNodeIdsToLocate(_node_global->findIdForEach(globalNodeIdsCaptured->begin(),globalNodeIdsCaptured->end()));
  90       MCAuto<DataArrayIdType> localCellCaptured(_mesh->getCellIdsLyingOnNodes(localNodeIdsToLocate->begin(),localNodeIdsToLocate->end(),false));
  91       MCAuto<DataArrayIdType> localCellCapturedGlob(_cell_global->selectByTupleIdSafe(localCellCaptured->begin(),localCellCaptured->end()));
  92       if(tabs[curRk].isNull())
  93         tabs[curRk] = localCellCapturedGlob;
  94       else
  95         tabs[curRk]->insertAtTheEnd(localCellCapturedGlob->begin(),localCellCapturedGlob->end());
  96     }
  97   }
  98   for(int curRk = 0 ; curRk < size ; ++curRk)
  99   {
 100     tabs[curRk] = tabs[curRk]->buildUniqueNotSorted();
 101     nbOfElems3[curRk] = tabs[curRk]->getNumberOfTuples();
 102   }
 103   std::vector<const DataArrayIdType *> tabss(tabs.begin(),tabs.end());
 104   MCAuto<DataArrayIdType> cells(DataArrayIdType::Aggregate(tabss));
 105   ci.allToAll(nbOfElems3.get(),1,MPI_ID_TYPE,nbOfElems2.get(),1,MPI_ID_TYPE,comm);
 106   mcIdType nbOfCellIdsSum(std::accumulate(nbOfElems2.get(),nbOfElems2.get()+size,0));
 107   MCAuto<DataArrayIdType> cellIdsFromProcs(DataArrayIdType::New());
 108   cellIdsFromProcs->alloc(nbOfCellIdsSum,1);
 109   {
 110     std::unique_ptr<int[]> nbOfElemsInt( CommInterface::ToIntArray<mcIdType>(nbOfElems3,size) ),nbOfElemsOutInt( CommInterface::ToIntArray<mcIdType>(nbOfElems2,size) );
 111     std::unique_ptr<int[]> offsetsIn( CommInterface::ComputeOffset(nbOfElemsInt,size) ), offsetsOut( CommInterface::ComputeOffset(nbOfElemsOutInt,size) );
 112     ci.allToAllV(cells->begin(),nbOfElemsInt.get(),offsetsIn.get(),MPI_ID_TYPE,
 113                  cellIdsFromProcs->getPointer(),nbOfElemsOutInt.get(),offsetsOut.get(),MPI_ID_TYPE,comm);
 114   }
 115   cellIdsFromProcs->sort();
 116   return cellIdsFromProcs;
 117 }
 118
 119 /*!
 120  */
 121 MCAuto<ParaUMesh> ParaUMesh::redistributeCells(const DataArrayIdType *globalCellIds) const
 122 {
 123   MPI_Comm comm(MPI_COMM_WORLD);
 124   CommInterface ci;
 125   int size;
 126   ci.commSize(comm,&size);
 127   std::unique_ptr<mcIdType[]> nbOfElems(new mcIdType[size]);
 128   mcIdType nbOfCellsRequested(globalCellIds->getNumberOfTuples());
 129   ci.allGather(&nbOfCellsRequested,1,MPI_ID_TYPE,nbOfElems.get(),1,MPI_ID_TYPE,comm);
 130   mcIdType nbOfCellIdsSum(std::accumulate(nbOfElems.get(),nbOfElems.get()+size,0));
 131   std::unique_ptr<mcIdType[]> allGlobalCellIds(new mcIdType[nbOfCellIdsSum]);
 132   std::unique_ptr<int[]> nbOfElemsInt( CommInterface::ToIntArray<mcIdType>(nbOfElems,size) );
 133   std::unique_ptr<int[]> offsetsIn( CommInterface::ComputeOffset(nbOfElemsInt,size) );
 134   ci.allGatherV(globalCellIds->begin(),nbOfCellsRequested,MPI_ID_TYPE,allGlobalCellIds.get(),nbOfElemsInt.get(),offsetsIn.get(),MPI_ID_TYPE,comm);
 135   mcIdType offset(0);
 136   // Prepare ParaUMesh parts to be sent : compute for each proc the contribution of current rank.
 137   std::vector< MCAuto<DataArrayIdType> > globalCellIdsToBeSent(size),globalNodeIdsToBeSent(size);
 138   std::vector< MCAuto<MEDCouplingUMesh> > meshPartsToBeSent(size);
 139   for(int curRk = 0 ; curRk < size ; ++curRk)
 140   {
 141     MCAuto<DataArrayIdType> globalCellIdsOfCurProc(DataArrayIdType::New());
 142     globalCellIdsOfCurProc->useArray(allGlobalCellIds.get()+offset,false,DeallocType::CPP_DEALLOC,nbOfElems[curRk],1);
 143     offset += nbOfElems[curRk];
 144     // the key call is here : compute for rank curRk the cells to be sent
 145     MCAuto<DataArrayIdType> globalCellIdsCaptured(_cell_global->buildIntersection(globalCellIdsOfCurProc));// OK for the global cellIds
 146     MCAuto<DataArrayIdType> localCellIdsCaptured(_node_global->findIdForEach(globalCellIdsCaptured->begin(),globalCellIdsCaptured->end()));
 147     MCAuto<MEDCouplingUMesh> meshPart(_mesh->buildPartOfMySelf(localCellIdsCaptured->begin(),localCellIdsCaptured->end(),true));
 148     MCAuto<DataArrayIdType> o2n(meshPart->zipCoordsTraducer());// OK for the mesh
 149     MCAuto<DataArrayIdType> n2o(o2n->invertArrayO2N2N2O(meshPart->getNumberOfNodes()));
 150     MCAuto<DataArrayIdType> globalNodeIdsPart(_node_global->selectByTupleIdSafe(n2o->begin(),n2o->end())); // OK for the global nodeIds
 151     meshPartsToBeSent[curRk] = meshPart;
 152     globalCellIdsToBeSent[curRk] = globalCellIdsCaptured;
 153     globalNodeIdsToBeSent[curRk] = globalNodeIdsPart;
 154   }
 155   // Receive
 156 }