1 // Copyright (C) 2007-2013 CEA/DEN, EDF R&D
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.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 // Author : Anthony Geay (CEA/DEN)
21 #include "OverlapElementLocator.hxx"
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"
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),
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();
58 OverlapElementLocator::~OverlapElementLocator()
60 delete [] _domain_bounding_boxes;
63 const MPI_Comm *OverlapElementLocator::getCommunicator() const
65 const MPIProcessorGroup* group=static_cast<const MPIProcessorGroup*>(&_group);
66 return group->getComm();
69 void OverlapElementLocator::computeBoundingBoxes()
71 CommInterface comm_interface=_group.getCommInterface();
72 const MPIProcessorGroup* group=static_cast<const MPIProcessorGroup*> (&_group);
74 if(_local_source_mesh)
75 _local_space_dim=_local_source_mesh->getSpaceDimension();
77 _local_space_dim=_local_target_mesh->getSpaceDimension();
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);
88 for(int i=0;i<_local_space_dim;i++)
90 minmax[i*2]=std::numeric_limits<double>::max();
91 minmax[i*2+1]=-std::numeric_limits<double>::max();
94 if(_local_target_mesh)
95 _local_target_mesh->getBoundingBox(minmax+2*_local_space_dim);
98 for(int i=0;i<_local_space_dim;i++)
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();
104 comm_interface.allGather(minmax, bbSize, MPI_DOUBLE,
105 _domain_bounding_boxes,bbSize, MPI_DOUBLE,
108 // Computation of all pairs needing an interpolation pairs are duplicated now !
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++)
115 if(intersectsBoundingBox(i,j))
116 _proc_pairs[i].push_back(j);
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());
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++)
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));
128 pairsToBeDonePerProc[*it2].push_back(std::pair<int,int>(i,*it2));
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.
135 int myProcId=_group.myRank();
136 _to_do_list=pairsToBeDonePerProc[myProcId];
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--)
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++)
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));
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.
158 void OverlapElementLocator::exchangeMeshes(OverlapInterpolationMatrix& matrix)
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++)
165 if((*it).first!=(*it).second)
167 if((*it).first==myProcId)
169 if((*it).second<myProcId)
170 receiveRemoteMesh((*it).second,false);
172 toDoListForFetchRemaining.push_back(std::pair<int,int>((*it).first,(*it).second));
175 {//(*it).second==myProcId
176 if((*it).first<myProcId)
177 receiveRemoteMesh((*it).first,true);
179 toDoListForFetchRemaining.push_back(std::pair<int,int>((*it).first,(*it).second));
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++)
189 if((*it).first!=(*it).second)
191 if((*it).first==myProcId)
192 receiveRemoteMesh((*it).second,false);
193 else//(*it).second==myProcId
194 receiveRemoteMesh((*it).first,true);
199 std::string OverlapElementLocator::getSourceMethod() const
201 return _local_source_field->getField()->getDiscretization()->getStringRepr();
204 std::string OverlapElementLocator::getTargetMethod() const
206 return _local_target_field->getField()->getDiscretization()->getStringRepr();
209 const MEDCouplingPointSet *OverlapElementLocator::getSourceMesh(int procId) const
211 int myProcId=_group.myRank();
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);
219 const DataArrayInt *OverlapElementLocator::getSourceIds(int procId) const
221 int myProcId=_group.myRank();
224 std::pair<int,bool> p(procId,true);
225 std::map<std::pair<int,bool>, MEDCouplingAutoRefCountObjectPtr< DataArrayInt > >::const_iterator it=_remote_elems.find(p);
229 const MEDCouplingPointSet *OverlapElementLocator::getTargetMesh(int procId) const
231 int myProcId=_group.myRank();
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);
239 const DataArrayInt *OverlapElementLocator::getTargetIds(int procId) const
241 int myProcId=_group.myRank();
244 std::pair<int,bool> p(procId,false);
245 std::map<std::pair<int,bool>, MEDCouplingAutoRefCountObjectPtr< DataArrayInt > >::const_iterator it=_remote_elems.find(p);
249 bool OverlapElementLocator::intersectsBoundingBox(int isource, int itarget) const
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;
254 for (int idim=0; idim < _local_space_dim; idim++)
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);
266 * This methods sends local source if 'sourceOrTarget'==True to proc 'procId'.
267 * This methods sends local target if 'sourceOrTarget'==False to proc 'procId'.
269 * This method prepares the matrix too, for matrix assembling and future matrix-vector computation.
271 void OverlapElementLocator::sendLocalMeshTo(int procId, bool sourceOrTarget, OverlapInterpolationMatrix& matrix) const
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
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;
283 else//target for local but source for distant
285 distant_bb=_domain_bounding_boxes+procId*2*2*_local_space_dim;
286 local_mesh=_local_target_mesh;
287 field=_local_target_field;
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));
293 matrix.keepTracksOfSourceIds(procId,idsToSend);//Case#1 in Step2 of main algorithm.
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();
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
305 void OverlapElementLocator::receiveRemoteMesh(int procId, bool sourceOrTarget)
308 MEDCouplingPointSet *m=0;
309 receiveMesh(procId,m,da);
310 std::pair<int,bool> p(procId,sourceOrTarget);
315 void OverlapElementLocator::sendMesh(int procId, const MEDCouplingPointSet *mesh, const DataArrayInt *idsToSend) const
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();
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);
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);
343 void OverlapElementLocator::receiveMesh(int procId, MEDCouplingPointSet* &mesh, DataArrayInt *&ids) const
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);
366 v1Distant->decrRef();
367 v2Distant->decrRef();