src/ParaMEDMEM/OverlapElementLocator.cxx

   1 // Copyright (C) 2007-2015  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 // Author : Anthony Geay (CEA/DEN)
  20
  21 #include "OverlapElementLocator.hxx"
  22
  23 #include "CommInterface.hxx"
  24 #include "Topology.hxx"
  25 #include "BlockTopology.hxx"
  26 #include "ParaFIELD.hxx"
  27 #include "ParaMESH.hxx"
  28 #include "ProcessorGroup.hxx"
  29 #include "MPIProcessorGroup.hxx"
  30 #include "OverlapInterpolationMatrix.hxx"
  31 #include "MEDCouplingFieldDouble.hxx"
  32 #include "MEDCouplingFieldDiscretization.hxx"
  33 #include "DirectedBoundingBox.hxx"
  34 #include "InterpKernelAutoPtr.hxx"
  35
  36 #include <limits>
  37
  38 using namespace std;
  39
  40 namespace ParaMEDMEM
  41 {
  42   OverlapElementLocator::OverlapElementLocator(const ParaFIELD *sourceField, const ParaFIELD *targetField, const ProcessorGroup& group)
  43     : _local_source_field(sourceField),
  44       _local_target_field(targetField),
  45       _local_source_mesh(0),
  46       _local_target_mesh(0),
  47       _domain_bounding_boxes(0),
  48       _group(group)
  49   {
  50     if(_local_source_field)
  51       _local_source_mesh=_local_source_field->getSupport()->getCellMesh();
  52     if(_local_target_field)
  53       _local_target_mesh=_local_target_field->getSupport()->getCellMesh();
  54     _comm=getCommunicator();
  55     computeBoundingBoxes();
  56   }
  57
  58   OverlapElementLocator::~OverlapElementLocator()
  59   {
  60     delete [] _domain_bounding_boxes;
  61   }
  62
  63   const MPI_Comm *OverlapElementLocator::getCommunicator() const
  64   {
  65     const MPIProcessorGroup* group=static_cast<const MPIProcessorGroup*>(&_group);
  66     return group->getComm();
  67   }
  68
  69   void OverlapElementLocator::computeBoundingBoxes()
  70   {
  71     CommInterface comm_interface=_group.getCommInterface();
  72     const MPIProcessorGroup* group=static_cast<const MPIProcessorGroup*> (&_group);
  73     _local_space_dim=0;
  74     if(_local_source_mesh)
  75       _local_space_dim=_local_source_mesh->getSpaceDimension();
  76     else
  77       _local_space_dim=_local_target_mesh->getSpaceDimension();
  78     //
  79     const MPI_Comm* comm = group->getComm();
  80     int bbSize=2*2*_local_space_dim;//2 (for source/target) 2 (min/max)
  81     _domain_bounding_boxes=new double[bbSize*_group.size()];
  82     INTERP_KERNEL::AutoPtr<double> minmax=new double[bbSize];
  83     //Format minmax : Xmin_src,Xmax_src,Ymin_src,Ymax_src,Zmin_src,Zmax_src,Xmin_trg,Xmax_trg,Ymin_trg,Ymax_trg,Zmin_trg,Zmax_trg
  84     if(_local_source_mesh)
  85       _local_source_mesh->getBoundingBox(minmax);
  86     else
  87       {
  88         for(int i=0;i<_local_space_dim;i++)
  89           {
  90             minmax[i*2]=std::numeric_limits<double>::max();
  91             minmax[i*2+1]=-std::numeric_limits<double>::max();
  92           }
  93       }
  94     if(_local_target_mesh)
  95       _local_target_mesh->getBoundingBox(minmax+2*_local_space_dim);
  96     else
  97       {
  98         for(int i=0;i<_local_space_dim;i++)
  99           {
 100             minmax[i*2+2*_local_space_dim]=std::numeric_limits<double>::max();
 101             minmax[i*2+1+2*_local_space_dim]=-std::numeric_limits<double>::max();
 102           }
 103       }
 104     comm_interface.allGather(minmax, bbSize, MPI_DOUBLE,
 105                              _domain_bounding_boxes,bbSize, MPI_DOUBLE,
 106                              *comm);
 107
 108     // Computation of all pairs needing an interpolation pairs are duplicated now !
 109
 110     _proc_pairs.clear();//first is source second is target
 111     _proc_pairs.resize(_group.size());
 112     for(int i=0;i<_group.size();i++)
 113       for(int j=0;j<_group.size();j++)
 114         {
 115           if(intersectsBoundingBox(i,j))
 116             _proc_pairs[i].push_back(j);
 117         }
 118
 119     // OK now let's assigning as balanced as possible, job to each proc of group
 120     std::vector< std::vector< std::pair<int,int> > > pairsToBeDonePerProc(_group.size());
 121     int i=0;
 122     for(std::vector< std::vector< int > >::const_iterator it1=_proc_pairs.begin();it1!=_proc_pairs.end();it1++,i++)
 123       for(std::vector< int >::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
 124         {
 125           if(pairsToBeDonePerProc[i].size()<=pairsToBeDonePerProc[*it2].size())//it includes the fact that i==*it2
 126             pairsToBeDonePerProc[i].push_back(std::pair<int,int>(i,*it2));
 127           else
 128             pairsToBeDonePerProc[*it2].push_back(std::pair<int,int>(i,*it2));
 129         }
 130     //Keeping todo list of current proc. _to_do_list contains a set of pair where at least _group.myRank() appears once.
 131     //This proc will be in charge to perform interpolation of any of element of '_to_do_list'
 132     //If _group.myRank()==myPair.first, current proc should fetch target mesh of myPair.second (if different from _group.myRank()).
 133     //If _group.myRank()==myPair.second, current proc should fetch source mesh of myPair.second.
 134
 135     int myProcId=_group.myRank();
 136     _to_do_list=pairsToBeDonePerProc[myProcId];
 137
 138     //Feeding now '_procs_to_send'. A same id can appears twice. The second parameter in pair means what to send true=source, false=target
 139     _procs_to_send.clear();
 140     for(int i=_group.size()-1;i>=0;i--)
 141       if(i!=myProcId)
 142         {
 143           const std::vector< std::pair<int,int> >& anRemoteProcToDoList=pairsToBeDonePerProc[i];
 144           for(std::vector< std::pair<int,int> >::const_iterator it=anRemoteProcToDoList.begin();it!=anRemoteProcToDoList.end();it++)
 145             {
 146               if((*it).first==myProcId)
 147                 _procs_to_send.push_back(std::pair<int,bool>(i,true));
 148               if((*it).second==myProcId)
 149                 _procs_to_send.push_back(std::pair<int,bool>(i,false));
 150             }
 151         }
 152   }
 153
 154   /*!
 155    * The aim of this method is to perform the communication to get data corresponding to '_to_do_list' attribute.
 156    * The principle is the following : if proc n1 and n2 need to perform a cross sending with n1<n2, then n1 will send first and receive then.
 157    */
 158   void OverlapElementLocator::exchangeMeshes(OverlapInterpolationMatrix& matrix)
 159   {
 160     int myProcId=_group.myRank();
 161     //starting to receive every procs whose id is lower than myProcId.
 162     std::vector< std::pair<int,int> > toDoListForFetchRemaining;
 163     for(std::vector< std::pair<int,int> >::const_iterator it=_to_do_list.begin();it!=_to_do_list.end();it++)
 164       {
 165         if((*it).first!=(*it).second)
 166           {
 167             if((*it).first==myProcId)
 168               {
 169                 if((*it).second<myProcId)
 170                   receiveRemoteMesh((*it).second,false);
 171                 else
 172                   toDoListForFetchRemaining.push_back(std::pair<int,int>((*it).first,(*it).second));
 173               }
 174             else
 175               {//(*it).second==myProcId
 176                 if((*it).first<myProcId)
 177                   receiveRemoteMesh((*it).first,true);
 178                 else
 179                   toDoListForFetchRemaining.push_back(std::pair<int,int>((*it).first,(*it).second));
 180               }
 181           }
 182       }
 183     //sending source or target mesh to remote procs
 184     for(std::vector< std::pair<int,bool> >::const_iterator it2=_procs_to_send.begin();it2!=_procs_to_send.end();it2++)
 185       sendLocalMeshTo((*it2).first,(*it2).second,matrix);
 186     //fetching remaining meshes
 187     for(std::vector< std::pair<int,int> >::const_iterator it=toDoListForFetchRemaining.begin();it!=toDoListForFetchRemaining.end();it++)
 188       {
 189         if((*it).first!=(*it).second)
 190           {
 191             if((*it).first==myProcId)
 192               receiveRemoteMesh((*it).second,false);
 193             else//(*it).second==myProcId
 194               receiveRemoteMesh((*it).first,true);
 195           }
 196       }
 197   }
 198
 199   std::string OverlapElementLocator::getSourceMethod() const
 200   {
 201     return _local_source_field->getField()->getDiscretization()->getStringRepr();
 202   }
 203
 204   std::string OverlapElementLocator::getTargetMethod() const
 205   {
 206     return _local_target_field->getField()->getDiscretization()->getStringRepr();
 207   }
 208
 209   const MEDCouplingPointSet *OverlapElementLocator::getSourceMesh(int procId) const
 210   {
 211     int myProcId=_group.myRank();
 212     if(myProcId==procId)
 213       return _local_source_mesh;
 214     std::pair<int,bool> p(procId,true);
 215     std::map<std::pair<int,bool>, MEDCouplingAutoRefCountObjectPtr< MEDCouplingPointSet > >::const_iterator it=_remote_meshes.find(p);
 216     return (*it).second;
 217   }
 218
 219   const DataArrayInt *OverlapElementLocator::getSourceIds(int procId) const
 220   {
 221     int myProcId=_group.myRank();
 222     if(myProcId==procId)
 223       return 0;
 224     std::pair<int,bool> p(procId,true);
 225     std::map<std::pair<int,bool>, MEDCouplingAutoRefCountObjectPtr< DataArrayInt > >::const_iterator it=_remote_elems.find(p);
 226     return (*it).second;
 227   }
 228
 229   const MEDCouplingPointSet *OverlapElementLocator::getTargetMesh(int procId) const
 230   {
 231     int myProcId=_group.myRank();
 232     if(myProcId==procId)
 233       return _local_target_mesh;
 234     std::pair<int,bool> p(procId,false);
 235     std::map<std::pair<int,bool>, MEDCouplingAutoRefCountObjectPtr< MEDCouplingPointSet > >::const_iterator it=_remote_meshes.find(p);
 236     return (*it).second;
 237   }
 238
 239   const DataArrayInt *OverlapElementLocator::getTargetIds(int procId) const
 240   {
 241     int myProcId=_group.myRank();
 242     if(myProcId==procId)
 243       return 0;
 244     std::pair<int,bool> p(procId,false);
 245     std::map<std::pair<int,bool>, MEDCouplingAutoRefCountObjectPtr< DataArrayInt > >::const_iterator it=_remote_elems.find(p);
 246     return (*it).second;
 247   }
 248
 249   bool OverlapElementLocator::intersectsBoundingBox(int isource, int itarget) const
 250   {
 251     const double *source_bb=_domain_bounding_boxes+isource*2*2*_local_space_dim;
 252     const double *target_bb=_domain_bounding_boxes+itarget*2*2*_local_space_dim+2*_local_space_dim;
 253
 254     for (int idim=0; idim < _local_space_dim; idim++)
 255       {
 256         const double eps = -1e-12;//tony to change
 257         bool intersects = (target_bb[idim*2]<source_bb[idim*2+1]+eps)
 258           && (source_bb[idim*2]<target_bb[idim*2+1]+eps);
 259         if (!intersects)
 260           return false;
 261       }
 262     return true;
 263   }
 264
 265   /*!
 266    * This methods sends local source if 'sourceOrTarget'==True to proc 'procId'.
 267    * This methods sends local target if 'sourceOrTarget'==False to proc 'procId'.
 268    *
 269    * This method prepares the matrix too, for matrix assembling and future matrix-vector computation.
 270    */
 271   void OverlapElementLocator::sendLocalMeshTo(int procId, bool sourceOrTarget, OverlapInterpolationMatrix& matrix) const
 272   {
 273    //int myProcId=_group.myRank();
 274    const double *distant_bb=0;
 275    MEDCouplingPointSet *local_mesh=0;
 276    const ParaFIELD *field=0;
 277    if(sourceOrTarget)//source for local but target for distant
 278      {
 279        distant_bb=_domain_bounding_boxes+procId*2*2*_local_space_dim+2*_local_space_dim;
 280        local_mesh=_local_source_mesh;
 281        field=_local_source_field;
 282      }
 283    else//target for local but source for distant
 284      {
 285        distant_bb=_domain_bounding_boxes+procId*2*2*_local_space_dim;
 286        local_mesh=_local_target_mesh;
 287        field=_local_target_field;
 288      }
 289    MEDCouplingAutoRefCountObjectPtr<DataArrayInt> elems=local_mesh->getCellsInBoundingBox(distant_bb,getBoundingBoxAdjustment());
 290    DataArrayInt *idsToSend;
 291    MEDCouplingPointSet *send_mesh=static_cast<MEDCouplingPointSet *>(field->getField()->buildSubMeshData(elems->begin(),elems->end(),idsToSend));
 292    if(sourceOrTarget)
 293      matrix.keepTracksOfSourceIds(procId,idsToSend);//Case#1 in Step2 of main algorithm.
 294    else
 295      matrix.keepTracksOfTargetIds(procId,idsToSend);//Case#0 in Step2 of main algorithm.
 296    sendMesh(procId,send_mesh,idsToSend);
 297    send_mesh->decrRef();
 298    idsToSend->decrRef();
 299   }
 300
 301   /*!
 302    * This method recieves source remote mesh on proc 'procId' if sourceOrTarget==True
 303    * This method recieves target remote mesh on proc 'procId' if sourceOrTarget==False
 304    */
 305   void OverlapElementLocator::receiveRemoteMesh(int procId, bool sourceOrTarget)
 306   {
 307     DataArrayInt *da=0;
 308     MEDCouplingPointSet *m=0;
 309     receiveMesh(procId,m,da);
 310     std::pair<int,bool> p(procId,sourceOrTarget);
 311     _remote_meshes[p]=m;
 312     _remote_elems[p]=da;
 313   }
 314
 315   void OverlapElementLocator::sendMesh(int procId, const MEDCouplingPointSet *mesh, const DataArrayInt *idsToSend) const
 316   {
 317     CommInterface comInterface=_group.getCommInterface();
 318     // First stage : exchanging sizes
 319     vector<double> tinyInfoLocalD;//tinyInfoLocalD not used for the moment
 320     vector<int> tinyInfoLocal;
 321     vector<string> tinyInfoLocalS;
 322     mesh->getTinySerializationInformation(tinyInfoLocalD,tinyInfoLocal,tinyInfoLocalS);
 323     const MPI_Comm *comm=getCommunicator();
 324     //
 325     int lgth[2];
 326     lgth[0]=tinyInfoLocal.size();
 327     lgth[1]=idsToSend->getNbOfElems();
 328     comInterface.send(&lgth,2,MPI_INT,procId,1140,*_comm);
 329     comInterface.send(&tinyInfoLocal[0],tinyInfoLocal.size(),MPI_INT,procId,1141,*comm);
 330     //
 331     DataArrayInt *v1Local=0;
 332     DataArrayDouble *v2Local=0;
 333     mesh->serialize(v1Local,v2Local);
 334     comInterface.send(v1Local->getPointer(),v1Local->getNbOfElems(),MPI_INT,procId,1142,*comm);
 335     comInterface.send(v2Local->getPointer(),v2Local->getNbOfElems(),MPI_DOUBLE,procId,1143,*comm);
 336     //finished for mesh, ids now
 337     comInterface.send(const_cast<int *>(idsToSend->getConstPointer()),lgth[1],MPI_INT,procId,1144,*comm);
 338     //
 339     v1Local->decrRef();
 340     v2Local->decrRef();
 341   }
 342
 343   void OverlapElementLocator::receiveMesh(int procId, MEDCouplingPointSet* &mesh, DataArrayInt *&ids) const
 344   {
 345     int lgth[2];
 346     MPI_Status status;
 347     const MPI_Comm *comm=getCommunicator();
 348     CommInterface comInterface=_group.getCommInterface();
 349     comInterface.recv(lgth,2,MPI_INT,procId,1140,*_comm,&status);
 350     std::vector<int> tinyInfoDistant(lgth[0]);
 351     ids=DataArrayInt::New();
 352     ids->alloc(lgth[1],1);
 353     comInterface.recv(&tinyInfoDistant[0],lgth[0],MPI_INT,procId,1141,*comm,&status);
 354     mesh=MEDCouplingPointSet::BuildInstanceFromMeshType((MEDCouplingMeshType)tinyInfoDistant[0]);
 355     std::vector<std::string> unusedTinyDistantSts;
 356     vector<double> tinyInfoDistantD(1);//tinyInfoDistantD not used for the moment
 357     DataArrayInt *v1Distant=DataArrayInt::New();
 358     DataArrayDouble *v2Distant=DataArrayDouble::New();
 359     mesh->resizeForUnserialization(tinyInfoDistant,v1Distant,v2Distant,unusedTinyDistantSts);
 360     comInterface.recv(v1Distant->getPointer(),v1Distant->getNbOfElems(),MPI_INT,procId,1142,*comm,&status);
 361     comInterface.recv(v2Distant->getPointer(),v2Distant->getNbOfElems(),MPI_DOUBLE,procId,1143,*comm,&status);
 362     mesh->unserialization(tinyInfoDistantD,tinyInfoDistant,v1Distant,v2Distant,unusedTinyDistantSts);
 363     //finished for mesh, ids now
 364     comInterface.recv(ids->getPointer(),lgth[1],MPI_INT,procId,1144,*comm,&status);
 365     //
 366     v1Distant->decrRef();
 367     v2Distant->decrRef();
 368   }
 369 }