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