Back to index

salome-med  6.5.0
OverlapElementLocator.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D
00002 //
00003 // This library is free software; you can redistribute it and/or
00004 // modify it under the terms of the GNU Lesser General Public
00005 // License as published by the Free Software Foundation; either
00006 // version 2.1 of the License.
00007 //
00008 // This library is distributed in the hope that it will be useful,
00009 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00010 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00011 // Lesser General Public License for more details.
00012 //
00013 // You should have received a copy of the GNU Lesser General Public
00014 // License along with this library; if not, write to the Free Software
00015 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00016 //
00017 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00018 //
00019 
00020 #include "OverlapElementLocator.hxx"
00021 
00022 #include "CommInterface.hxx"
00023 #include "Topology.hxx"
00024 #include "BlockTopology.hxx"
00025 #include "ParaFIELD.hxx"
00026 #include "ParaMESH.hxx"
00027 #include "ProcessorGroup.hxx"
00028 #include "MPIProcessorGroup.hxx"
00029 #include "OverlapInterpolationMatrix.hxx"
00030 #include "MEDCouplingFieldDouble.hxx"
00031 #include "MEDCouplingFieldDiscretization.hxx"
00032 #include "DirectedBoundingBox.hxx"
00033 #include "InterpKernelAutoPtr.hxx"
00034 
00035 #include <limits>
00036 
00037 using namespace std;
00038 
00039 namespace ParaMEDMEM 
00040 { 
00041   OverlapElementLocator::OverlapElementLocator(const ParaFIELD *sourceField, const ParaFIELD *targetField, const ProcessorGroup& group)
00042     : _local_source_field(sourceField),
00043       _local_target_field(targetField),
00044       _local_source_mesh(0),
00045       _local_target_mesh(0),
00046       _domain_bounding_boxes(0),
00047       _group(group)
00048   { 
00049     if(_local_source_field)
00050       _local_source_mesh=_local_source_field->getSupport()->getCellMesh();
00051     if(_local_target_field)
00052       _local_target_mesh=_local_target_field->getSupport()->getCellMesh();
00053     _comm=getCommunicator();
00054     computeBoundingBoxes();
00055   }
00056 
00057   OverlapElementLocator::~OverlapElementLocator()
00058   {
00059     delete [] _domain_bounding_boxes;
00060   }
00061 
00062   const MPI_Comm *OverlapElementLocator::getCommunicator() const
00063   {
00064     const MPIProcessorGroup* group=static_cast<const MPIProcessorGroup*>(&_group);
00065     return group->getComm();
00066   }
00067 
00068   void OverlapElementLocator::computeBoundingBoxes()
00069   {
00070     CommInterface comm_interface=_group.getCommInterface();
00071     const MPIProcessorGroup* group=static_cast<const MPIProcessorGroup*> (&_group);
00072     _local_space_dim=0;
00073     if(_local_source_mesh)
00074       _local_space_dim=_local_source_mesh->getSpaceDimension();
00075     else
00076       _local_space_dim=_local_target_mesh->getSpaceDimension();
00077     //
00078     const MPI_Comm* comm = group->getComm();
00079     int bbSize=2*2*_local_space_dim;//2 (for source/target) 2 (min/max)
00080     _domain_bounding_boxes=new double[bbSize*_group.size()];
00081     INTERP_KERNEL::AutoPtr<double> minmax=new double[bbSize];
00082     //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
00083     if(_local_source_mesh)
00084       _local_source_mesh->getBoundingBox(minmax);
00085     else
00086       {
00087         for(int i=0;i<_local_space_dim;i++)
00088           {
00089             minmax[i*2]=std::numeric_limits<double>::max();
00090             minmax[i*2+1]=-std::numeric_limits<double>::max();
00091           }
00092       }
00093     if(_local_target_mesh)
00094       _local_target_mesh->getBoundingBox(minmax+2*_local_space_dim);
00095     else
00096       {
00097         for(int i=0;i<_local_space_dim;i++)
00098           {
00099             minmax[i*2+2*_local_space_dim]=std::numeric_limits<double>::max();
00100             minmax[i*2+1+2*_local_space_dim]=-std::numeric_limits<double>::max();
00101           }
00102       }
00103     comm_interface.allGather(minmax, bbSize, MPI_DOUBLE,
00104                              _domain_bounding_boxes,bbSize, MPI_DOUBLE, 
00105                              *comm);
00106   
00107     // Computation of all pairs needing an interpolation pairs are duplicated now !
00108     
00109     _proc_pairs.clear();//first is source second is target
00110     _proc_pairs.resize(_group.size());
00111     for(int i=0;i<_group.size();i++)
00112       for(int j=0;j<_group.size();j++)
00113         {
00114           if(intersectsBoundingBox(i,j))
00115             _proc_pairs[i].push_back(j);
00116         }
00117 
00118     // OK now let's assigning as balanced as possible, job to each proc of group
00119     std::vector< std::vector< std::pair<int,int> > > pairsToBeDonePerProc(_group.size());
00120     int i=0;
00121     for(std::vector< std::vector< int > >::const_iterator it1=_proc_pairs.begin();it1!=_proc_pairs.end();it1++,i++)
00122       for(std::vector< int >::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
00123         {
00124           if(pairsToBeDonePerProc[i].size()<=pairsToBeDonePerProc[*it2].size())//it includes the fact that i==*it2
00125             pairsToBeDonePerProc[i].push_back(std::pair<int,int>(i,*it2));
00126           else
00127             pairsToBeDonePerProc[*it2].push_back(std::pair<int,int>(i,*it2));
00128         }
00129     //Keeping todo list of current proc. _to_do_list contains a set of pair where at least _group.myRank() appears once.
00130     //This proc will be in charge to perform interpolation of any of element of '_to_do_list'
00131     //If _group.myRank()==myPair.first, current proc should fetch target mesh of myPair.second (if different from _group.myRank()).
00132     //If _group.myRank()==myPair.second, current proc should fetch source mesh of myPair.second.
00133     
00134     int myProcId=_group.myRank();
00135     _to_do_list=pairsToBeDonePerProc[myProcId];
00136 
00137     //Feeding now '_procs_to_send'. A same id can appears twice. The second parameter in pair means what to send true=source, false=target
00138     _procs_to_send.clear();
00139     for(int i=_group.size()-1;i>=0;i--)
00140       if(i!=myProcId)
00141         {
00142           const std::vector< std::pair<int,int> >& anRemoteProcToDoList=pairsToBeDonePerProc[i];
00143           for(std::vector< std::pair<int,int> >::const_iterator it=anRemoteProcToDoList.begin();it!=anRemoteProcToDoList.end();it++)
00144             {
00145               if((*it).first==myProcId)
00146                 _procs_to_send.push_back(std::pair<int,bool>(i,true));
00147               if((*it).second==myProcId)
00148                 _procs_to_send.push_back(std::pair<int,bool>(i,false));
00149             }
00150         }
00151   }
00152 
00157   void OverlapElementLocator::exchangeMeshes(OverlapInterpolationMatrix& matrix)
00158   {
00159     int myProcId=_group.myRank();
00160     //starting to receive every procs whose id is lower than myProcId.
00161     std::vector< std::pair<int,int> > toDoListForFetchRemaining;
00162     for(std::vector< std::pair<int,int> >::const_iterator it=_to_do_list.begin();it!=_to_do_list.end();it++)
00163       {
00164         if((*it).first!=(*it).second)
00165           {
00166             if((*it).first==myProcId)
00167               {
00168                 if((*it).second<myProcId)
00169                   receiveRemoteMesh((*it).second,false);
00170                 else
00171                   toDoListForFetchRemaining.push_back(std::pair<int,int>((*it).first,(*it).second));
00172               }
00173             else
00174               {//(*it).second==myProcId
00175                 if((*it).first<myProcId)
00176                   receiveRemoteMesh((*it).first,true);
00177                 else
00178                   toDoListForFetchRemaining.push_back(std::pair<int,int>((*it).first,(*it).second));
00179               }
00180           }
00181       }
00182     //sending source or target mesh to remote procs
00183     for(std::vector< std::pair<int,bool> >::const_iterator it2=_procs_to_send.begin();it2!=_procs_to_send.end();it2++)
00184       sendLocalMeshTo((*it2).first,(*it2).second,matrix);
00185     //fetching remaining meshes
00186     for(std::vector< std::pair<int,int> >::const_iterator it=toDoListForFetchRemaining.begin();it!=toDoListForFetchRemaining.end();it++)
00187       {
00188         if((*it).first!=(*it).second)
00189           {
00190             if((*it).first==myProcId)
00191               receiveRemoteMesh((*it).second,false);
00192             else//(*it).second==myProcId
00193               receiveRemoteMesh((*it).first,true);
00194           }
00195       }
00196   }
00197   
00198   std::string OverlapElementLocator::getSourceMethod() const
00199   {
00200     return _local_source_field->getField()->getDiscretization()->getStringRepr();
00201   }
00202 
00203   std::string OverlapElementLocator::getTargetMethod() const
00204   {
00205     return _local_target_field->getField()->getDiscretization()->getStringRepr();
00206   }
00207 
00208   const MEDCouplingPointSet *OverlapElementLocator::getSourceMesh(int procId) const
00209   {
00210     int myProcId=_group.myRank();
00211     if(myProcId==procId)
00212       return _local_source_mesh;
00213     std::pair<int,bool> p(procId,true);
00214     std::map<std::pair<int,bool>, MEDCouplingAutoRefCountObjectPtr< MEDCouplingPointSet > >::const_iterator it=_remote_meshes.find(p);
00215     return (*it).second;
00216   }
00217 
00218   const DataArrayInt *OverlapElementLocator::getSourceIds(int procId) const
00219   {
00220     int myProcId=_group.myRank();
00221     if(myProcId==procId)
00222       return 0;
00223     std::pair<int,bool> p(procId,true);
00224     std::map<std::pair<int,bool>, MEDCouplingAutoRefCountObjectPtr< DataArrayInt > >::const_iterator it=_remote_elems.find(p);
00225     return (*it).second;
00226   }
00227 
00228   const MEDCouplingPointSet *OverlapElementLocator::getTargetMesh(int procId) const
00229   {
00230     int myProcId=_group.myRank();
00231     if(myProcId==procId)
00232       return _local_target_mesh;
00233     std::pair<int,bool> p(procId,false);
00234     std::map<std::pair<int,bool>, MEDCouplingAutoRefCountObjectPtr< MEDCouplingPointSet > >::const_iterator it=_remote_meshes.find(p);
00235     return (*it).second;
00236   }
00237 
00238   const DataArrayInt *OverlapElementLocator::getTargetIds(int procId) const
00239   {
00240     int myProcId=_group.myRank();
00241     if(myProcId==procId)
00242       return 0;
00243     std::pair<int,bool> p(procId,false);
00244     std::map<std::pair<int,bool>, MEDCouplingAutoRefCountObjectPtr< DataArrayInt > >::const_iterator it=_remote_elems.find(p);
00245     return (*it).second;
00246   }
00247 
00248   bool OverlapElementLocator::intersectsBoundingBox(int isource, int itarget) const
00249   {
00250     const double *source_bb=_domain_bounding_boxes+isource*2*2*_local_space_dim;
00251     const double *target_bb=_domain_bounding_boxes+itarget*2*2*_local_space_dim+2*_local_space_dim;
00252 
00253     for (int idim=0; idim < _local_space_dim; idim++)
00254       {
00255         const double eps = -1e-12;//tony to change
00256         bool intersects = (target_bb[idim*2]<source_bb[idim*2+1]+eps)
00257           && (source_bb[idim*2]<target_bb[idim*2+1]+eps);
00258         if (!intersects)
00259           return false; 
00260       }
00261     return true;
00262   }
00263 
00270   void OverlapElementLocator::sendLocalMeshTo(int procId, bool sourceOrTarget, OverlapInterpolationMatrix& matrix) const
00271   {
00272    vector<int> elems;
00273    //int myProcId=_group.myRank();
00274    const double *distant_bb=0;
00275    MEDCouplingPointSet *local_mesh=0;
00276    const ParaFIELD *field=0;
00277    if(sourceOrTarget)//source for local but target for distant
00278      {
00279        distant_bb=_domain_bounding_boxes+procId*2*2*_local_space_dim+2*_local_space_dim;
00280        local_mesh=_local_source_mesh;
00281        field=_local_source_field;
00282      }
00283    else//target for local but source for distant
00284      {
00285        distant_bb=_domain_bounding_boxes+procId*2*2*_local_space_dim;
00286        local_mesh=_local_target_mesh;
00287        field=_local_target_field;
00288      }
00289    local_mesh->getCellsInBoundingBox(distant_bb,getBoundingBoxAdjustment(),elems);
00290    DataArrayInt *idsToSend;
00291    MEDCouplingPointSet *send_mesh=static_cast<MEDCouplingPointSet *>(field->getField()->buildSubMeshData(&elems[0],&elems[elems.size()],idsToSend));
00292    if(sourceOrTarget)
00293      matrix.keepTracksOfSourceIds(procId,idsToSend);//Case#1 in Step2 of main algorithm.
00294    else
00295      matrix.keepTracksOfTargetIds(procId,idsToSend);//Case#0 in Step2 of main algorithm.
00296    sendMesh(procId,send_mesh,idsToSend);
00297    send_mesh->decrRef();
00298    idsToSend->decrRef();
00299   }
00300 
00305   void OverlapElementLocator::receiveRemoteMesh(int procId, bool sourceOrTarget)
00306   {
00307     DataArrayInt *da=0;
00308     MEDCouplingPointSet *m=0;
00309     receiveMesh(procId,m,da);
00310     std::pair<int,bool> p(procId,sourceOrTarget);
00311     _remote_meshes[p]=m;
00312     _remote_elems[p]=da;
00313   }
00314 
00315   void OverlapElementLocator::sendMesh(int procId, const MEDCouplingPointSet *mesh, const DataArrayInt *idsToSend) const
00316   {
00317     CommInterface comInterface=_group.getCommInterface();
00318     // First stage : exchanging sizes
00319     vector<double> tinyInfoLocalD;//tinyInfoLocalD not used for the moment
00320     vector<int> tinyInfoLocal;
00321     vector<string> tinyInfoLocalS;
00322     mesh->getTinySerializationInformation(tinyInfoLocalD,tinyInfoLocal,tinyInfoLocalS);
00323     const MPI_Comm *comm=getCommunicator();
00324     //
00325     int lgth[2];
00326     lgth[0]=tinyInfoLocal.size();
00327     lgth[1]=idsToSend->getNbOfElems();
00328     comInterface.send(&lgth,2,MPI_INT,procId,1140,*_comm);
00329     comInterface.send(&tinyInfoLocal[0],tinyInfoLocal.size(),MPI_INT,procId,1141,*comm);
00330     //
00331     DataArrayInt *v1Local=0;
00332     DataArrayDouble *v2Local=0;
00333     mesh->serialize(v1Local,v2Local);
00334     comInterface.send(v1Local->getPointer(),v1Local->getNbOfElems(),MPI_INT,procId,1142,*comm);
00335     comInterface.send(v2Local->getPointer(),v2Local->getNbOfElems(),MPI_DOUBLE,procId,1143,*comm);
00336     //finished for mesh, ids now
00337     comInterface.send((int *)idsToSend->getConstPointer(),lgth[1],MPI_INT,procId,1144,*comm);
00338     //
00339     v1Local->decrRef();
00340     v2Local->decrRef();
00341   }
00342 
00343   void OverlapElementLocator::receiveMesh(int procId, MEDCouplingPointSet* &mesh, DataArrayInt *&ids) const
00344   {
00345     int lgth[2];
00346     MPI_Status status;
00347     const MPI_Comm *comm=getCommunicator();
00348     CommInterface comInterface=_group.getCommInterface();
00349     comInterface.recv(lgth,2,MPI_INT,procId,1140,*_comm,&status);
00350     std::vector<int> tinyInfoDistant(lgth[0]);
00351     ids=DataArrayInt::New();
00352     ids->alloc(lgth[1],1);
00353     comInterface.recv(&tinyInfoDistant[0],lgth[0],MPI_INT,procId,1141,*comm,&status);
00354     mesh=MEDCouplingPointSet::BuildInstanceFromMeshType((MEDCouplingMeshType)tinyInfoDistant[0]);
00355     std::vector<std::string> unusedTinyDistantSts;
00356     vector<double> tinyInfoDistantD(1);//tinyInfoDistantD not used for the moment
00357     DataArrayInt *v1Distant=DataArrayInt::New();
00358     DataArrayDouble *v2Distant=DataArrayDouble::New();
00359     mesh->resizeForUnserialization(tinyInfoDistant,v1Distant,v2Distant,unusedTinyDistantSts);
00360     comInterface.recv(v1Distant->getPointer(),v1Distant->getNbOfElems(),MPI_INT,procId,1142,*comm,&status);
00361     comInterface.recv(v2Distant->getPointer(),v2Distant->getNbOfElems(),MPI_DOUBLE,procId,1143,*comm,&status);
00362     mesh->unserialization(tinyInfoDistantD,tinyInfoDistant,v1Distant,v2Distant,unusedTinyDistantSts);
00363     //finished for mesh, ids now
00364     comInterface.recv(ids->getPointer(),lgth[1],MPI_INT,procId,1144,*comm,&status);
00365     //
00366     v1Distant->decrRef();
00367     v2Distant->decrRef();
00368   }
00369 }