src/MEDPartitioner/MEDPARTITIONER_JointFinder.cxx

   1 // Copyright (C) 2007-2021  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19
  20 #include "MEDPARTITIONER_JointFinder.hxx"
  21 #include "MEDPARTITIONER_MeshCollection.hxx"
  22 #include "MEDPARTITIONER_Topology.hxx"
  23 #include "MEDPARTITIONER_ParaDomainSelector.hxx"
  24 #include "MEDPARTITIONER_Utils.hxx"
  25
  26 #include "MEDCouplingUMesh.hxx"
  27 #include "BBTree.txx"
  28
  29 /*!
  30  * Method contributing to the distant cell graph
  31  */
  32 MEDPARTITIONER::JointFinder::JointFinder(const MeshCollection& mc):_mesh_collection(mc),_domain_selector(mc.getParaDomainSelector()),_topology(mc.getTopology())
  33 {
  34 }
  35
  36 MEDPARTITIONER::JointFinder::~JointFinder()
  37 {
  38 }
  39
  40 void MEDPARTITIONER::JointFinder::findCommonDistantNodes()
  41 {
  42   int nbdomain=_topology->nbDomain();
  43   _distant_node_cell.resize(nbdomain);
  44   _node_node.resize(nbdomain);
  45   for (int i=0; i<nbdomain; i++)
  46     {
  47       _distant_node_cell[i].resize(nbdomain);
  48       _node_node[i].resize(nbdomain);
  49     }
  50   int nbproc=_domain_selector->nbProcs();
  51   std::vector<BBTreeOfDim* > bbtree(nbdomain,(BBTreeOfDim*) 0);
  52   std::vector<double* > bbxi(nbdomain,(double*) 0);
  53   std::vector<MEDCoupling::DataArrayIdType*> rev(nbdomain,(MEDCoupling::DataArrayIdType*) 0);
  54   std::vector<MEDCoupling::DataArrayIdType*> revIndx(nbdomain,(MEDCoupling::DataArrayIdType*) 0);
  55   //int meshDim=-1;
  56   int spaceDim=-1;
  57
  58   //init rev and revIndx and bbtree for my domain (of me:proc n)
  59   for (int mydomain=0; mydomain<nbdomain; mydomain++)
  60     {
  61       if(!_domain_selector->isMyDomain(mydomain))
  62         continue;
  63       const MEDCoupling::MEDCouplingUMesh* myMesh=_mesh_collection.getMesh(mydomain);
  64       //meshDim = myMesh->getMeshDimension();
  65       spaceDim= myMesh->getSpaceDimension();
  66       rev[mydomain] = MEDCoupling::DataArrayIdType::New();
  67       revIndx[mydomain] = MEDCoupling::DataArrayIdType::New();
  68       myMesh->getReverseNodalConnectivity(rev[mydomain],revIndx[mydomain]);
  69       double* bbx=new double[2*spaceDim*myMesh->getNumberOfNodes()];
  70       for (int i=0; i<myMesh->getNumberOfNodes()*spaceDim; i++)
  71         {
  72           const double* coords=myMesh->getCoords()->getConstPointer();
  73           bbx[2*i]=(coords[i])-1e-12;
  74           bbx[2*i+1]=bbx[2*i]+2e-12;
  75         }
  76       bbtree[mydomain]=new BBTreeOfDim( spaceDim, bbx,0,0,myMesh->getNumberOfNodes(),-1e-12);
  77       //keep bbx because need it in getIntersectingElems
  78       //no delete [] bbx yet
  79       bbxi[mydomain]=bbx;
  80     }
  81
  82   //send my domains to other proc an receive other domains from other proc
  83   for (int isource=0; isource<nbdomain; isource++)
  84     {
  85       for (int itarget=0; itarget<nbdomain; itarget++)
  86         {
  87           const MEDCoupling::MEDCouplingUMesh* sourceMesh=_mesh_collection.getMesh(isource);
  88           if (_domain_selector->isMyDomain(isource)&&_domain_selector->isMyDomain(itarget))
  89             continue;
  90           if (_domain_selector->isMyDomain(isource))
  91             {
  92               //preparing data for treatment on target proc
  93               int targetProc = _domain_selector->getProcessorID(itarget);
  94
  95               std::vector<double> vec(spaceDim*sourceMesh->getNumberOfNodes());
  96               std::copy(sourceMesh->getCoords()->getConstPointer(),sourceMesh->getCoords()->getConstPointer()+sourceMesh->getNumberOfNodes()*spaceDim,&vec[0]);
  97               SendDoubleVec(vec,targetProc);
  98
  99               //retrieving target data for storage in commonDistantNodes array
 100               std::vector<mcIdType> localCorrespondency;
 101               RecvIntVec(localCorrespondency, targetProc);
 102               for (std::size_t i=0; i<localCorrespondency.size()/2; i++)
 103                 {
 104                   _distant_node_cell[isource][itarget].insert(std::make_pair(localCorrespondency[2*i],localCorrespondency[2*i+1]));
 105                 }
 106
 107             }
 108
 109           if (_domain_selector->isMyDomain(itarget))
 110             {
 111               //receiving data from source proc
 112               int sourceProc = isource%nbproc;
 113               std::vector<double> recvVec;
 114               RecvDoubleVec(recvVec,sourceProc);
 115               std::map<mcIdType,mcIdType> commonNodes; // (local nodes, distant nodes) list
 116               for (mcIdType inode=0; inode<ToIdType(recvVec.size()/spaceDim); inode++)
 117                 {
 118                   double* bbox=new double[2*spaceDim];
 119                   for (int i=0; i<spaceDim; i++)
 120                     {
 121                       bbox[2*i]=recvVec[inode*spaceDim+i]-1e-12;
 122                       bbox[2*i+1]=bbox[2*i]+2e-12;
 123                     }
 124                   std::vector<mcIdType> inodes;
 125                   bbtree[itarget]->getIntersectingElems(bbox,inodes);
 126                   delete [] bbox;
 127
 128                   if (inodes.size()>0)
 129                     {
 130                       commonNodes.insert(std::make_pair(inodes[0],inode));
 131                     }
 132
 133                 }
 134               std::vector<mcIdType> nodeCellCorrespondency;
 135               for (std::map<mcIdType,mcIdType>::iterator iter=commonNodes.begin(); iter!=commonNodes.end(); iter++)
 136                 {
 137                   _node_node[itarget][isource].push_back(std::make_pair(iter->first, iter->second));//storing node pairs in a vector
 138                   const mcIdType* revIndxPtr=revIndx[itarget]->getConstPointer();
 139                   const mcIdType* revPtr=rev[itarget]->getConstPointer();
 140                   for (mcIdType icell=revIndxPtr[iter->first]; icell<revIndxPtr[iter->first+1]; icell++)
 141                     {
 142                       nodeCellCorrespondency.push_back(iter->second); //
 143                       mcIdType globalCell=_topology->convertCellToGlobal(itarget,revPtr[icell]);
 144                       nodeCellCorrespondency.push_back(globalCell);
 145                     }
 146                 }
 147               SendIntVec(nodeCellCorrespondency, sourceProc); //itarget proc send to other (otherLocalNode-itargetGlobalCell)
 148             }
 149         }
 150     }
 151
 152   //free  rev(nbdomain) revIndx(nbdomain) bbtree(nbdomain) bbxi(nbdomain)
 153   for (int i=0; i<nbdomain; i++)
 154     {
 155       if (rev[i]!=0)
 156         rev[i]->decrRef();
 157       if (revIndx[i]!=0)
 158         revIndx[i]->decrRef();
 159       if (bbtree[i]!=0)
 160         delete bbtree[i];
 161       if (bbxi[i]!=0)
 162         delete [] bbxi[i];
 163     }
 164
 165   if (MyGlobals::_Verbose>100)
 166     std::cout << "proc " << _domain_selector->rank() << " : end JointFinder::findCommonDistantNodes" << std::endl;
 167 }
 168
 169 std::vector<std::vector<std::multimap<mcIdType,mcIdType> > >& MEDPARTITIONER::JointFinder::getDistantNodeCell()
 170 {
 171   return _distant_node_cell;
 172 }
 173
 174 std::vector<std::vector<std::vector<std::pair<mcIdType,mcIdType> > > >& MEDPARTITIONER::JointFinder::getNodeNode()
 175 {
 176   return _node_node;
 177 }
 178
 179 void MEDPARTITIONER::JointFinder::print()
 180 //it is for debug on small arrays under mpi 2,3 cpus
 181 {
 182   int nbdomain=_topology->nbDomain();
 183   //MPI_Barrier(MPI_COMM_WORLD);
 184   if (MyGlobals::_Is0verbose>0)
 185     std::cout << "\nJointFinder print node-node (nn)iproc|itarget|isource|i|inodefirst-inodesecond\n\n" <<
 186       "JointFinder print distantNode=cell (nc)iproc|itarget|isource|inode=icell\n\n";
 187   for (int isource=0; isource<nbdomain; isource++)
 188     {
 189       for (int itarget=0; itarget<nbdomain; itarget++)
 190         {
 191           for (std::size_t i=0; i<_node_node[itarget][isource].size(); i++)
 192             std::cout << " nn" << _domain_selector->rank() << itarget << "|" << isource << "|" << i << "|" <<
 193               _node_node[itarget][isource][i].first << "-" <<
 194               _node_node[itarget][isource][i].second;
 195         }
 196     }
 197   std::cout<<std::endl;
 198   //MPI_Barrier(MPI_COMM_WORLD);
 199   for (int isource=0; isource<nbdomain; isource++)
 200     {
 201       for (int itarget=0; itarget<nbdomain; itarget++)
 202         {
 203           std::multimap<mcIdType,mcIdType>::iterator it;
 204           for (it=_distant_node_cell[isource][itarget].begin() ; it!=_distant_node_cell[isource][itarget].end(); it++)
 205             {
 206               std::cout << " nc" << _domain_selector->rank() << "|" << itarget << "|" << isource << "|" << (*it).first << "=" << (*it).second;
 207             }
 208         }
 209     }
 210   std::cout << std::endl;
 211   //MPI_Barrier(MPI_COMM_WORLD);
 212 }