1 // Copyright (C) 2007-2015 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, or (at your option) any later version.
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 const int OverlapElementLocator::START_TAG_MESH_XCH = 1140;
44 OverlapElementLocator::OverlapElementLocator(const ParaFIELD *sourceField, const ParaFIELD *targetField,
45 const ProcessorGroup& group, double epsAbs, int workSharingAlgo)
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),
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();
60 computeBoundingBoxesAndInteractionList();
61 if (workSharingAlgo == 0)
62 computeTodoList_original();
64 if(workSharingAlgo == 1)
65 computeTodoList_new();
67 throw INTERP_KERNEL::Exception("OverlapElementLocator::OverlapElementLocator(): invalid algorithm selected!");
72 OverlapElementLocator::~OverlapElementLocator()
74 delete [] _domain_bounding_boxes;
77 const MPI_Comm *OverlapElementLocator::getCommunicator() const
79 const MPIProcessorGroup* group=static_cast<const MPIProcessorGroup*>(&_group);
80 return group->getComm();
83 void OverlapElementLocator::computeBoundingBoxesAndInteractionList()
85 CommInterface comm_interface=_group.getCommInterface();
86 const MPIProcessorGroup* group=static_cast<const MPIProcessorGroup*> (&_group);
88 if(_local_source_mesh)
89 _local_space_dim=_local_source_mesh->getSpaceDimension();
91 _local_space_dim=_local_target_mesh->getSpaceDimension();
93 const MPI_Comm* comm = group->getComm();
94 int bbSize=2*2*_local_space_dim;//2 (for source/target) 2 (min/max)
95 _domain_bounding_boxes=new double[bbSize*_group.size()];
96 INTERP_KERNEL::AutoPtr<double> minmax=new double[bbSize];
97 //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
98 if(_local_source_mesh)
99 _local_source_mesh->getBoundingBox(minmax);
102 for(int i=0;i<_local_space_dim;i++)
104 minmax[i*2]=std::numeric_limits<double>::max();
105 minmax[i*2+1]=-std::numeric_limits<double>::max();
108 if(_local_target_mesh)
109 _local_target_mesh->getBoundingBox(minmax+2*_local_space_dim);
112 for(int i=0;i<_local_space_dim;i++)
114 minmax[i*2+2*_local_space_dim]=std::numeric_limits<double>::max();
115 minmax[i*2+1+2*_local_space_dim]=-std::numeric_limits<double>::max();
118 comm_interface.allGather(minmax, bbSize, MPI_DOUBLE,
119 _domain_bounding_boxes,bbSize, MPI_DOUBLE,
122 // Computation of all pairs needing an interpolation pairs are duplicated now !
124 _proc_pairs.clear();//first is source second is target
125 _proc_pairs.resize(_group.size());
126 for(int i=0;i<_group.size();i++)
127 for(int j=0;j<_group.size();j++)
128 if(intersectsBoundingBox(i,j))
129 _proc_pairs[i].push_back(j);
132 void OverlapElementLocator::computeTodoList_original()
134 // OK now let's assigning as balanced as possible, job to each proc of group
135 _all_todo_lists.resize(_group.size());
137 for(std::vector< std::vector< int > >::const_iterator it1=_proc_pairs.begin();it1!=_proc_pairs.end();it1++,i++)
138 for(std::vector< int >::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
140 if(_all_todo_lists[i].size()<=_all_todo_lists[*it2].size())//it includes the fact that i==*it2
141 _all_todo_lists[i].push_back(ProcCouple(i,*it2));
143 _all_todo_lists[*it2].push_back(ProcCouple(i,*it2));
145 //Keeping todo list of current proc. _to_do_list contains a set of pair where at least _group.myRank() appears once.
146 //This proc will be in charge to perform interpolation of any of element of '_to_do_list'
147 //If _group.myRank()==myPair.first, current proc should fetch target mesh of myPair.second (if different from _group.myRank()).
148 //If _group.myRank()==myPair.second, current proc should fetch source mesh of myPair.second.
150 int myProcId=_group.myRank();
151 _to_do_list=_all_todo_lists[myProcId];
154 std::stringstream scout;
155 scout << "(" << myProcId << ") my TODO list is: ";
156 for (std::vector< ProcCouple >::const_iterator itdbg=_to_do_list.begin(); itdbg!=_to_do_list.end(); itdbg++)
157 scout << "(" << (*itdbg).first << "," << (*itdbg).second << ")";
158 std::cout << scout.str() << "\n";
162 /* More efficient (?) work sharing algorithm: a job (i,j) is initially assigned twice: to proc#i and to proc#j.
163 * Then try to reduce as much as possible the variance of the num of jobs per proc:
164 * - take the most loaded proc i,
165 * + select the job (i,j) for which proc#j is the less loaded
166 * + remove this job from proc#i
167 * - repeat until no more duplicates are found
169 void OverlapElementLocator::computeTodoList_new()
172 int infinity = std::numeric_limits<int>::max();
174 int grp_size = _group.size();
175 vector < map<ProcCouple, int> > full_set(grp_size );
177 for(vector< vector< int > >::const_iterator it = _proc_pairs.begin(); it != _proc_pairs.end(); it++, srcProcID++)
178 for (vector< int >::const_iterator it2=(*it).begin(); it2 != (*it).end(); it2++)
180 // Here a pair of the form (i,i) is added only once!
181 int tgtProcID = *it2;
182 ProcCouple cpl = make_pair(srcProcID, tgtProcID);
183 full_set[srcProcID][cpl] = -1;
184 full_set[tgtProcID][cpl] = -1;
187 vector < map<ProcCouple, int> > ::iterator itVector;
188 map<ProcCouple, int>::iterator itMap;
189 for(itVector = full_set.begin(); itVector != full_set.end(); itVector++, procID++)
190 for (itMap=(*itVector).begin(); itMap != (*itVector).end(); itMap++)
192 const ProcCouple & cpl = (*itMap).first;
193 if (cpl.first == cpl.second)
194 // special case - this couple can not be removed in the future
195 (*itMap).second = infinity;
198 if(cpl.first == procID)
199 (*itMap).second = full_set[cpl.second].size();
200 else // cpl.second == srcProcID
201 (*itMap).second = full_set[cpl.first].size();
204 INTERP_KERNEL::AutoPtr<bool> proc_valid = new bool[grp_size];
205 fill((bool *)proc_valid, proc_valid+grp_size, true);
208 while (find((bool *)proc_valid, proc_valid+grp_size, true) != proc_valid+grp_size)
211 int max_sz = -1, max_id = -1;
212 for(itVector = full_set.begin(), procID=0; itVector != full_set.end(); itVector++, procID++)
214 int sz = (*itVector).size();
215 if (proc_valid[procID] && sz > max_sz)
217 max_sz = (*itVector).size();
222 // Nothing more to do:
225 // For this proc, job with less loaded second proc:
226 int min_sz = infinity;
227 map<ProcCouple, int> & max_map = full_set[max_id];
228 ProcCouple hit_cpl = make_pair(-1,-1);
229 for(itMap=max_map.begin(); itMap != max_map.end(); itMap++)
230 if ((*itMap).second < min_sz)
231 hit_cpl = (*itMap).first;
232 if (hit_cpl.first == -1)
234 // Plouf. Current proc 'max_id' can not be reduced. Invalid it:
235 proc_valid[max_id] = false;
238 // Remove item from proc 'max_id'
239 full_set[max_id].erase(hit_cpl);
240 // And mark it as not removable on the other side:
241 if (hit_cpl.first == max_id)
242 full_set[hit_cpl.second][hit_cpl] = infinity;
243 else // hit_cpl.second == max_id
244 full_set[hit_cpl.first][hit_cpl] = infinity;
246 // Now update all counts of valid maps:
248 for(itVector = full_set.begin(); itVector != full_set.end(); itVector++, procID++)
249 if(proc_valid[procID] && procID != max_id)
250 for(itMap = (*itVector).begin(); itMap != (*itVector).end(); itMap++)
252 const ProcCouple & cpl = (*itMap).first;
254 if ((*itMap).second == infinity)
256 if (cpl.first == max_id || cpl.second == max_id)
260 // Final formatting - extract remaining keys in each map:
261 int myProcId=_group.myRank();
262 _all_todo_lists.resize(grp_size);
264 for(itVector = full_set.begin(); itVector != full_set.end(); itVector++, procID++)
265 for(itMap = (*itVector).begin(); itMap != (*itVector).end(); itMap++)
266 _all_todo_lists[procID].push_back((*itMap).first);
267 _to_do_list=_all_todo_lists[myProcId];
270 std::stringstream scout;
271 scout << "(" << myProcId << ") my TODO list is: ";
272 for (std::vector< ProcCouple >::const_iterator itdbg=_to_do_list.begin(); itdbg!=_to_do_list.end(); itdbg++)
273 scout << "(" << (*itdbg).first << "," << (*itdbg).second << ")";
274 std::cout << scout.str() << "\n";
278 void OverlapElementLocator::fillProcToSend()
280 // Feeding now '_procs_to_send*'. A same id can appears twice. The second parameter in pair means what
281 // to send true=source, false=target
282 int myProcId=_group.myRank();
283 _procs_to_send_mesh.clear();
284 _procs_to_send_field.clear();
285 for(int i=_group.size()-1;i>=0;i--)
287 const std::vector< ProcCouple >& anRemoteProcToDoList=_all_todo_lists[i];
288 for(std::vector< ProcCouple >::const_iterator it=anRemoteProcToDoList.begin();it!=anRemoteProcToDoList.end();it++)
290 if((*it).first==myProcId)
293 _procs_to_send_mesh.push_back(Proc_SrcOrTgt(i,true));
294 _procs_to_send_field.push_back((*it).second);
296 if((*it).second==myProcId)
298 _procs_to_send_mesh.push_back(Proc_SrcOrTgt(i,false));
305 * The aim of this method is to perform the communication to get data corresponding to '_to_do_list' attribute.
306 * 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.
308 void OverlapElementLocator::exchangeMeshes(OverlapInterpolationMatrix& matrix)
310 int myProcId=_group.myRank();
311 //starting to receive every procs whose id is lower than myProcId.
312 std::vector< ProcCouple > toDoListForFetchRemaining;
313 for(std::vector< ProcCouple >::const_iterator it=_to_do_list.begin();it!=_to_do_list.end();it++)
315 int first = (*it).first, second = (*it).second;
321 receiveRemoteMeshFrom(second,false);
323 toDoListForFetchRemaining.push_back(ProcCouple(first,second));
326 {//(*it).second==myProcId
328 receiveRemoteMeshFrom(first,true);
330 toDoListForFetchRemaining.push_back(ProcCouple(first,second));
334 //sending source or target mesh to remote procs
335 for(std::vector< Proc_SrcOrTgt >::const_iterator it2=_procs_to_send_mesh.begin();it2!=_procs_to_send_mesh.end();it2++)
336 sendLocalMeshTo((*it2).first,(*it2).second,matrix);
337 //fetching remaining meshes
338 for(std::vector< ProcCouple >::const_iterator it=toDoListForFetchRemaining.begin();it!=toDoListForFetchRemaining.end();it++)
340 if((*it).first!=(*it).second)
342 if((*it).first==myProcId)
343 receiveRemoteMeshFrom((*it).second,false);
344 else//(*it).second==myProcId
345 receiveRemoteMeshFrom((*it).first,true);
350 std::string OverlapElementLocator::getSourceMethod() const
352 return _local_source_field->getField()->getDiscretization()->getStringRepr();
355 std::string OverlapElementLocator::getTargetMethod() const
357 return _local_target_field->getField()->getDiscretization()->getStringRepr();
360 const MEDCouplingPointSet *OverlapElementLocator::getSourceMesh(int procId) const
362 int myProcId=_group.myRank();
364 return _local_source_mesh;
365 Proc_SrcOrTgt p(procId,true);
366 std::map<Proc_SrcOrTgt, AutoMCPointSet >::const_iterator it=_remote_meshes.find(p);
370 const DataArrayInt *OverlapElementLocator::getSourceIds(int procId) const
372 int myProcId=_group.myRank();
375 Proc_SrcOrTgt p(procId,true);
376 std::map<Proc_SrcOrTgt, AutoDAInt >::const_iterator it=_remote_elems.find(p);
380 const MEDCouplingPointSet *OverlapElementLocator::getTargetMesh(int procId) const
382 int myProcId=_group.myRank();
384 return _local_target_mesh;
385 Proc_SrcOrTgt p(procId,false);
386 std::map<Proc_SrcOrTgt, AutoMCPointSet >::const_iterator it=_remote_meshes.find(p);
390 const DataArrayInt *OverlapElementLocator::getTargetIds(int procId) const
392 int myProcId=_group.myRank();
395 Proc_SrcOrTgt p(procId,false);
396 std::map<Proc_SrcOrTgt, AutoDAInt >::const_iterator it=_remote_elems.find(p);
400 bool OverlapElementLocator::isInMyTodoList(int i, int j) const
402 ProcCouple cpl = std::make_pair(i, j);
403 return std::find(_to_do_list.begin(), _to_do_list.end(), cpl)!=_to_do_list.end();
406 bool OverlapElementLocator::intersectsBoundingBox(int isource, int itarget) const
408 const double *source_bb=_domain_bounding_boxes+isource*2*2*_local_space_dim;
409 const double *target_bb=_domain_bounding_boxes+itarget*2*2*_local_space_dim+2*_local_space_dim;
411 for (int idim=0; idim < _local_space_dim; idim++)
413 bool intersects = (target_bb[idim*2]<source_bb[idim*2+1]+_epsAbs)
414 && (source_bb[idim*2]<target_bb[idim*2+1]+_epsAbs);
422 * This methods sends (part of) local source if 'sourceOrTarget'==True to proc 'procId'.
423 * This methods sends (part of) local target if 'sourceOrTarget'==False to proc 'procId'.
425 * This method prepares the matrix too, for matrix assembling and future matrix-vector computation.
427 void OverlapElementLocator::sendLocalMeshTo(int procId, bool sourceOrTarget, OverlapInterpolationMatrix& matrix) const
429 //int myProcId=_group.myRank();
430 const double *distant_bb=0;
431 MEDCouplingPointSet *local_mesh=0;
432 const ParaFIELD *field=0;
433 if(sourceOrTarget)//source for local mesh but target for distant mesh
435 distant_bb=_domain_bounding_boxes+procId*2*2*_local_space_dim+2*_local_space_dim;
436 local_mesh=_local_source_mesh;
437 field=_local_source_field;
439 else//target for local but source for distant
441 distant_bb=_domain_bounding_boxes+procId*2*2*_local_space_dim;
442 local_mesh=_local_target_mesh;
443 field=_local_target_field;
445 AutoDAInt elems=local_mesh->getCellsInBoundingBox(distant_bb,getBoundingBoxAdjustment());
446 DataArrayInt *old2new_map;
447 MEDCouplingPointSet *send_mesh=static_cast<MEDCouplingPointSet *>(field->getField()->buildSubMeshData(elems->begin(),elems->end(),old2new_map));
449 matrix.keepTracksOfSourceIds(procId,old2new_map);
451 matrix.keepTracksOfTargetIds(procId,old2new_map);
452 sendMesh(procId,send_mesh,old2new_map);
453 send_mesh->decrRef();
454 old2new_map->decrRef();
458 * This method recieves source remote mesh on proc 'procId' if sourceOrTarget==True
459 * This method recieves target remote mesh on proc 'procId' if sourceOrTarget==False
461 void OverlapElementLocator::receiveRemoteMeshFrom(int procId, bool sourceOrTarget)
463 DataArrayInt *old2new_map=0;
464 MEDCouplingPointSet *m=0;
465 receiveMesh(procId,m,old2new_map);
466 Proc_SrcOrTgt p(procId,sourceOrTarget);
468 _remote_elems[p]=old2new_map;
471 void OverlapElementLocator::sendMesh(int procId, const MEDCouplingPointSet *mesh, const DataArrayInt *idsToSend) const
473 CommInterface comInterface=_group.getCommInterface();
475 // First stage : exchanging sizes
476 vector<double> tinyInfoLocalD;//tinyInfoLocalD not used for the moment
477 vector<int> tinyInfoLocal;
478 vector<string> tinyInfoLocalS;
479 mesh->getTinySerializationInformation(tinyInfoLocalD,tinyInfoLocal,tinyInfoLocalS);
480 const MPI_Comm *comm=getCommunicator();
483 lgth[0]=tinyInfoLocal.size();
484 lgth[1]=idsToSend->getNbOfElems();
485 comInterface.send(&lgth,2,MPI_INT,procId,START_TAG_MESH_XCH,*_comm);
486 comInterface.send(&tinyInfoLocal[0],tinyInfoLocal.size(),MPI_INT,procId,START_TAG_MESH_XCH+1,*comm);
488 DataArrayInt *v1Local=0;
489 DataArrayDouble *v2Local=0;
490 mesh->serialize(v1Local,v2Local);
491 comInterface.send(v1Local->getPointer(),v1Local->getNbOfElems(),MPI_INT,procId,START_TAG_MESH_XCH+2,*comm);
492 comInterface.send(v2Local->getPointer(),v2Local->getNbOfElems(),MPI_DOUBLE,procId,START_TAG_MESH_XCH+3,*comm);
493 //finished for mesh, ids now
494 comInterface.send(const_cast<int *>(idsToSend->getConstPointer()),lgth[1],MPI_INT,procId,START_TAG_MESH_XCH+4,*comm);
500 void OverlapElementLocator::receiveMesh(int procId, MEDCouplingPointSet* &mesh, DataArrayInt *&ids) const
504 const MPI_Comm *comm=getCommunicator();
505 CommInterface comInterface=_group.getCommInterface();
506 comInterface.recv(lgth,2,MPI_INT,procId,START_TAG_MESH_XCH,*_comm,&status);
507 std::vector<int> tinyInfoDistant(lgth[0]);
508 ids=DataArrayInt::New();
509 ids->alloc(lgth[1],1);
510 comInterface.recv(&tinyInfoDistant[0],lgth[0],MPI_INT,procId,START_TAG_MESH_XCH+1,*comm,&status);
511 mesh=MEDCouplingPointSet::BuildInstanceFromMeshType((MEDCouplingMeshType)tinyInfoDistant[0]);
512 std::vector<std::string> unusedTinyDistantSts;
513 vector<double> tinyInfoDistantD(1);//tinyInfoDistantD not used for the moment
514 DataArrayInt *v1Distant=DataArrayInt::New();
515 DataArrayDouble *v2Distant=DataArrayDouble::New();
516 mesh->resizeForUnserialization(tinyInfoDistant,v1Distant,v2Distant,unusedTinyDistantSts);
517 comInterface.recv(v1Distant->getPointer(),v1Distant->getNbOfElems(),MPI_INT,procId,START_TAG_MESH_XCH+2,*comm,&status);
518 comInterface.recv(v2Distant->getPointer(),v2Distant->getNbOfElems(),MPI_DOUBLE,procId,START_TAG_MESH_XCH+3,*comm,&status);
519 mesh->unserialization(tinyInfoDistantD,tinyInfoDistant,v1Distant,v2Distant,unusedTinyDistantSts);
520 //finished for mesh, ids now
521 comInterface.recv(ids->getPointer(),lgth[1],MPI_INT,procId,1144,*comm,&status);
523 v1Distant->decrRef();
524 v2Distant->decrRef();