1 // Copyright (C) 2007-2024 CEA, EDF
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 switch(workSharingAlgo)
64 computeTodoList_original(); break;
66 computeTodoList_new(false); break;
68 computeTodoList_new(true); break;
70 throw INTERP_KERNEL::Exception("OverlapElementLocator::OverlapElementLocator(): invalid algorithm selected!");
76 OverlapElementLocator::~OverlapElementLocator()
78 delete [] _domain_bounding_boxes;
81 const MPI_Comm *OverlapElementLocator::getCommunicator() const
83 const MPIProcessorGroup* group=static_cast<const MPIProcessorGroup*>(&_group);
84 return group->getComm();
87 void OverlapElementLocator::computeBoundingBoxesAndInteractionList()
89 CommInterface comm_interface=_group.getCommInterface();
90 const MPIProcessorGroup* group=static_cast<const MPIProcessorGroup*> (&_group);
92 if(_local_source_mesh)
93 _local_space_dim=_local_source_mesh->getSpaceDimension();
95 _local_space_dim=_local_target_mesh->getSpaceDimension();
97 const MPI_Comm* comm = group->getComm();
98 int bbSize=2*2*_local_space_dim;//2 (for source/target) 2 (min/max)
99 _domain_bounding_boxes=new double[bbSize*_group.size()];
100 INTERP_KERNEL::AutoPtr<double> minmax=new double[bbSize];
101 //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
102 if(_local_source_mesh)
103 _local_source_mesh->getBoundingBox(minmax);
106 for(int i=0;i<_local_space_dim;i++)
108 minmax[i*2]=std::numeric_limits<double>::max();
109 minmax[i*2+1]=-std::numeric_limits<double>::max();
112 if(_local_target_mesh)
113 _local_target_mesh->getBoundingBox(minmax+2*_local_space_dim);
116 for(int i=0;i<_local_space_dim;i++)
118 minmax[i*2+2*_local_space_dim]=std::numeric_limits<double>::max();
119 minmax[i*2+1+2*_local_space_dim]=-std::numeric_limits<double>::max();
122 comm_interface.allGather(minmax, bbSize, MPI_DOUBLE,
123 _domain_bounding_boxes,bbSize, MPI_DOUBLE,
126 // Computation of all pairs needing an interpolation - pairs are duplicated now !
127 _proc_pairs.clear();//first is source second is target
128 _proc_pairs.resize(_group.size());
129 for(int i=0;i<_group.size();i++)
130 for(int j=0;j<_group.size();j++)
131 if(intersectsBoundingBox(i,j))
132 _proc_pairs[i].push_back(j);
135 /*! See main OverlapDEC documentation for an explanation on this one. This is the original work sharing algorithm.
137 void OverlapElementLocator::computeTodoList_original()
139 // OK now let's assigning as balanced as possible, job to each proc of group
140 _all_todo_lists.resize(_group.size());
141 for(int i = 0; i < _group.size(); i++)
142 for(const int j: _proc_pairs[i])
144 if(_all_todo_lists[i].size()<=_all_todo_lists[j].size())//it includes the fact that i==j
145 _all_todo_lists[i].push_back(ProcCouple(i,j));
147 _all_todo_lists[j].push_back(ProcCouple(i,j));
149 //Keeping todo list of current proc. _to_do_list contains a set of pair where at least _group.myRank() appears once.
150 //This proc will be in charge to perform interpolation of any of element of '_to_do_list'
151 //If _group.myRank()==myPair.first, current proc should fetch target mesh of myPair.second (if different from _group.myRank()).
152 //If _group.myRank()==myPair.second, current proc should fetch source mesh of myPair.second.
154 int myProcId=_group.myRank();
155 _to_do_list=_all_todo_lists[myProcId];
158 std::stringstream scout;
159 scout << "(" << myProcId << ") my TODO list is: ";
160 for (const ProcCouple& pc: _to_do_list)
161 scout << "(" << pc.first << "," << pc.second << ")";
162 std::cout << scout.str() << "\n";
166 /*! More efficient (?) work sharing algorithm: a job (i,j) is initially assigned twice: to proc#i and to proc#j.
167 * Then try to reduce as much as possible the variance of the num of jobs per proc by selecting the right duplicate
169 * - take the most loaded proc i,
170 * + select the job (i,j) for which proc#j is the less loaded
171 * + remove this job from proc#i, and mark it as 'unremovable' from proc#j
172 * - repeat until no more duplicates are found
174 void OverlapElementLocator::computeTodoList_new(bool revertIter)
177 int infinity = std::numeric_limits<int>::max();
182 int grp_size = _group.size();
183 // For each proc, a map giving for a job (=an interaction (i,j)) its 'load', i.e. the amount of work found on proc #j
184 vector < map<ProcCouple, int> > full_set(grp_size );
185 for(int srcProcID = 0; srcProcID < _group.size(); srcProcID++)
186 for(const int tgtProcID: _proc_pairs[srcProcID])
188 // Here a pair of the form (i,i) is added only once!
189 ProcCouple cpl = make_pair(srcProcID, tgtProcID);
190 // The interaction (i,j) is initially given to both procs #i and #j - load is initialised at -1:
191 full_set[srcProcID][cpl] = -1;
192 full_set[tgtProcID][cpl] = -1;
196 for(auto& itVector : full_set)
198 for (auto &mapIt : itVector)
200 const ProcCouple & cpl = mapIt.first;
201 if (cpl.first == cpl.second) // interaction (i,i) : can not be removed:
202 // special case - this couple can not be removed in the future
203 mapIt.second = infinity;
206 if(cpl.first == procID)
207 mapIt.second = (int)full_set[cpl.second].size();
208 else // cpl.second == srcProcID
209 mapIt.second = (int)full_set[cpl.first].size();
214 INTERP_KERNEL::AutoPtr<bool> proc_valid = new bool[grp_size];
215 fill((bool *)proc_valid, proc_valid+grp_size, true);
220 while (find((bool *)proc_valid, proc_valid+grp_size, true) != proc_valid+grp_size) // as long as proc_valid is not full of 'false'
223 int max_sz = -1, max_id = -1;
225 for(const auto& a_set: full_set)
227 int sz = (int)a_set.size();
228 if (proc_valid[procID] && sz > max_sz)
236 // Nothing more to do:
237 if (max_sz == -1) break;
238 // For this proc, job with less loaded second proc:
239 int min_sz = infinity;
240 map<ProcCouple, int> & max_map = full_set[max_id];
241 ProcCouple hit_cpl = make_pair(-1,-1);
244 // Use a reverse iterator here increases our chances to hit a couple of the form (i, myProcId)
245 // meaning that the final matrix computed won't have to be sent: save some comm.
246 for(auto ritMap=max_map.rbegin(); ritMap != max_map.rend(); ritMap++)
247 if ((*ritMap).second < min_sz)
249 hit_cpl = (*ritMap).first;
250 min_sz = (*ritMap).second;
255 for(const auto& mapIt : max_map)
256 if (mapIt.second < min_sz)
258 hit_cpl = mapIt.first;
259 min_sz = mapIt.second;
262 if (hit_cpl.first == -1)
264 // Plouf. Current proc 'max_id' can not be reduced. Invalid it and move next:
265 proc_valid[max_id] = false;
268 // Remove item from proc 'max_id'
269 full_set[max_id].erase(hit_cpl);
270 // And mark it as not removable on the other side:
271 if (hit_cpl.first == max_id)
272 full_set[hit_cpl.second][hit_cpl] = infinity;
273 else // hit_cpl.second == max_id
274 full_set[hit_cpl.first][hit_cpl] = infinity;
276 // Now update all counts of valid maps:
278 for(auto& itVector: full_set)
280 if(proc_valid[procID] && procID != max_id)
281 for(auto& mapIt: itVector)
283 const ProcCouple & cpl = mapIt.first;
284 if (mapIt.second == infinity) // Unremovable item:
286 if (cpl.first == max_id || cpl.second == max_id)
294 // Final formatting - extract remaining keys in each map:
296 int myProcId = _group.myRank();
297 _all_todo_lists.resize(grp_size);
299 for(const auto& itVector: full_set)
301 for(const auto& mapIt: itVector)
302 _all_todo_lists[procID].push_back(mapIt.first);
305 _to_do_list=_all_todo_lists[myProcId];
308 std::stringstream scout;
309 scout << "(" << myProcId << ") my TODO list is: ";
310 for (const ProcCouple& pc: _to_do_list)
311 scout << "(" << pc.first << "," << pc.second << ")";
312 std::cout << scout.str() << "\n";
316 void OverlapElementLocator::debugPrintWorkSharing(std::ostream & ostr) const
318 std::vector< std::vector< ProcCouple > >::const_iterator it = _all_todo_lists.begin();
319 ostr << "TODO list lengths: ";
320 for(; it != _all_todo_lists.end(); ++it)
321 ostr << (*it).size() << " ";
325 void OverlapElementLocator::fillProcToSend()
327 // Feeding now '_procs_to_send*'. A same id can appears twice. The second parameter in pair means what
328 // to send true=source, false=target
329 int myProcId=_group.myRank();
330 _procs_to_send_mesh.clear();
331 _procs_to_send_field.clear();
332 for(int i=0;i<_group.size();i++)
334 for(const ProcCouple& pc: _all_todo_lists[i])
336 if(pc.first==myProcId)
339 _procs_to_send_mesh.push_back(Proc_SrcOrTgt(i,true));
340 _procs_to_send_field.push_back(pc.second);
342 if(pc.second==myProcId)
344 _procs_to_send_mesh.push_back(Proc_SrcOrTgt(i,false));
348 std::stringstream scout;
349 scout << "(" << _group.myRank() << ") PROC TO SEND list is: ";
350 for (const auto& pc: _procs_to_send_mesh)
351 scout << "(" << pc.first << "," << (pc.second ? "src":"tgt") << ") ";
352 std::cout << scout.str() << "\n";
358 * The aim of this method is to perform the communication to get data corresponding to '_to_do_list' attribute.
359 * 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.
361 void OverlapElementLocator::exchangeMeshes(OverlapInterpolationMatrix& matrix)
363 int myProcId=_group.myRank();
364 //starting to receive every procs whose id is lower than myProcId.
365 std::vector<ProcCouple> toDoListForFetchRemaining;
366 for (const ProcCouple& pc: _to_do_list)
368 if(pc.first == pc.second) continue; // no xchg needed
370 if(pc.first==myProcId)
372 if(pc.second<myProcId)
373 receiveRemoteMeshFrom(pc.second,false);
375 toDoListForFetchRemaining.push_back(ProcCouple(pc.first, pc.second));
378 {//pc.second==myProcId
379 if(pc.first<myProcId)
380 receiveRemoteMeshFrom(pc.first,true);
382 toDoListForFetchRemaining.push_back(ProcCouple(pc.first, pc.second));
385 //sending source or target mesh to remote procs
386 for (const Proc_SrcOrTgt& pst: _procs_to_send_mesh)
387 sendLocalMeshTo(pst.first, pst.second,matrix);
388 //fetching remaining meshes
389 for (const ProcCouple& pc: toDoListForFetchRemaining)
390 { // NB: here pc.first always != from pc.second
391 if(pc.first==myProcId)
392 receiveRemoteMeshFrom(pc.second,false);
393 else//pc.second==myProcId
394 receiveRemoteMeshFrom(pc.first,true);
398 std::string OverlapElementLocator::getSourceMethod() const
400 return _local_source_field->getField()->getDiscretization()->getStringRepr();
403 std::string OverlapElementLocator::getTargetMethod() const
405 return _local_target_field->getField()->getDiscretization()->getStringRepr();
408 const MEDCouplingPointSet *OverlapElementLocator::getSourceMesh(int procId) const
410 int myProcId=_group.myRank();
412 return _local_source_mesh;
413 Proc_SrcOrTgt p(procId,true);
414 std::map<Proc_SrcOrTgt, AutoMCPointSet >::const_iterator it=_remote_meshes.find(p);
418 const DataArrayIdType *OverlapElementLocator::getSourceIds(int procId) const
420 int myProcId=_group.myRank();
423 Proc_SrcOrTgt p(procId,true);
424 std::map<Proc_SrcOrTgt, AutoDAInt >::const_iterator it=_remote_elems.find(p);
428 const MEDCouplingPointSet *OverlapElementLocator::getTargetMesh(int procId) const
430 int myProcId=_group.myRank();
432 return _local_target_mesh;
433 Proc_SrcOrTgt p(procId,false);
434 std::map<Proc_SrcOrTgt, AutoMCPointSet >::const_iterator it=_remote_meshes.find(p);
438 const DataArrayIdType *OverlapElementLocator::getTargetIds(int procId) const
440 int myProcId=_group.myRank();
443 Proc_SrcOrTgt p(procId,false);
444 std::map<Proc_SrcOrTgt, AutoDAInt >::const_iterator it=_remote_elems.find(p);
448 bool OverlapElementLocator::isInMyTodoList(int i, int j) const
450 ProcCouple cpl = std::make_pair(i, j);
451 return std::find(_to_do_list.begin(), _to_do_list.end(), cpl)!=_to_do_list.end();
454 bool OverlapElementLocator::intersectsBoundingBox(int isource, int itarget) const
456 const double *source_bb=_domain_bounding_boxes+isource*2*2*_local_space_dim;
457 const double *target_bb=_domain_bounding_boxes+itarget*2*2*_local_space_dim+2*_local_space_dim;
459 for (int idim=0; idim < _local_space_dim; idim++)
461 bool intersects = (target_bb[idim*2]<source_bb[idim*2+1]+_epsAbs)
462 && (source_bb[idim*2]<target_bb[idim*2+1]+_epsAbs);
470 * This methods sends (part of) local source if 'sourceOrTarget'==True to proc 'procId'.
471 * This methods sends (part of) local target if 'sourceOrTarget'==False to proc 'procId'.
473 * This method prepares the matrix too, for matrix assembling and future matrix-vector computation.
475 void OverlapElementLocator::sendLocalMeshTo(int procId, bool sourceOrTarget, OverlapInterpolationMatrix& matrix) const
478 int rank = _group.myRank();
479 std::string st = sourceOrTarget ? "src" : "tgt";
480 std::stringstream scout;
481 scout << "(" << rank << ") SEND part of " << st << " TO: " << procId;
482 std::cout << scout.str() << "\n";
485 //int myProcId=_group.myRank();
486 const double *distant_bb=0;
487 MEDCouplingPointSet *local_mesh=0;
488 const ParaFIELD *field=0;
489 if(sourceOrTarget)//source for local mesh but target for distant mesh
491 distant_bb=_domain_bounding_boxes+procId*2*2*_local_space_dim+2*_local_space_dim;
492 local_mesh=_local_source_mesh;
493 field=_local_source_field;
495 else//target for local but source for distant
497 distant_bb=_domain_bounding_boxes+procId*2*2*_local_space_dim;
498 local_mesh=_local_target_mesh;
499 field=_local_target_field;
501 AutoDAInt elems=local_mesh->getCellsInBoundingBox(distant_bb,getBoundingBoxAdjustment());
502 DataArrayIdType *old2new_map;
503 MEDCouplingPointSet *send_mesh=static_cast<MEDCouplingPointSet *>(field->getField()->buildSubMeshData(elems->begin(),elems->end(),old2new_map));
505 matrix.keepTracksOfSourceIds(procId,old2new_map);
507 matrix.keepTracksOfTargetIds(procId,old2new_map);
508 sendMesh(procId,send_mesh,old2new_map);
509 send_mesh->decrRef();
510 old2new_map->decrRef();
514 * This method receives source remote mesh on proc 'procId' if sourceOrTarget==True
515 * This method receives target remote mesh on proc 'procId' if sourceOrTarget==False
517 void OverlapElementLocator::receiveRemoteMeshFrom(int procId, bool sourceOrTarget)
520 int rank = _group.myRank();
521 std::string st = sourceOrTarget ? "src" : "tgt";
522 std::stringstream scout;
523 scout << "(" << rank << ") RCV part of " << st << " FROM: " << procId;
524 std::cout << scout.str() << "\n";
526 DataArrayIdType *old2new_map=0;
527 MEDCouplingPointSet *m=0;
528 receiveMesh(procId,m,old2new_map);
529 Proc_SrcOrTgt p(procId,sourceOrTarget);
531 _remote_elems[p]=old2new_map;
534 void OverlapElementLocator::sendMesh(int procId, const MEDCouplingPointSet *mesh, const DataArrayIdType *idsToSend) const
536 CommInterface comInterface=_group.getCommInterface();
538 // First stage : exchanging sizes
539 vector<double> tinyInfoLocalD;//tinyInfoLocalD not used for the moment
540 vector<mcIdType> tinyInfoLocal;
541 vector<string> tinyInfoLocalS;
542 mesh->getTinySerializationInformation(tinyInfoLocalD,tinyInfoLocal,tinyInfoLocalS);
543 const MPI_Comm *comm=getCommunicator();
546 lgth[0]=ToIdType(tinyInfoLocal.size());
547 lgth[1]=idsToSend->getNbOfElems();
548 comInterface.send(&lgth,2,MPI_ID_TYPE,procId,START_TAG_MESH_XCH,*_comm);
549 comInterface.send(&tinyInfoLocal[0],(int)tinyInfoLocal.size(),MPI_ID_TYPE,procId,START_TAG_MESH_XCH+1,*comm);
551 DataArrayIdType *v1Local=0;
552 DataArrayDouble *v2Local=0;
553 mesh->serialize(v1Local,v2Local);
554 comInterface.send(v1Local->getPointer(),(int)v1Local->getNbOfElems(),MPI_ID_TYPE,procId,START_TAG_MESH_XCH+2,*comm);
555 comInterface.send(v2Local->getPointer(),(int)v2Local->getNbOfElems(),MPI_DOUBLE,procId,START_TAG_MESH_XCH+3,*comm);
556 //finished for mesh, ids now
557 comInterface.send(const_cast<mcIdType *>(idsToSend->getConstPointer()),(int)lgth[1],MPI_ID_TYPE,procId,START_TAG_MESH_XCH+4,*comm);
563 void OverlapElementLocator::receiveMesh(int procId, MEDCouplingPointSet* &mesh, DataArrayIdType *&ids) const
567 const MPI_Comm *comm=getCommunicator();
568 CommInterface comInterface=_group.getCommInterface();
569 comInterface.recv(lgth,2,MPI_ID_TYPE,procId,START_TAG_MESH_XCH,*_comm,&status);
570 std::vector<mcIdType> tinyInfoDistant(lgth[0]);
571 ids=DataArrayIdType::New();
572 ids->alloc(lgth[1],1);
573 comInterface.recv(&tinyInfoDistant[0],(int)lgth[0],MPI_ID_TYPE,procId,START_TAG_MESH_XCH+1,*comm,&status);
574 mesh=MEDCouplingPointSet::BuildInstanceFromMeshType((MEDCouplingMeshType)tinyInfoDistant[0]);
575 std::vector<std::string> unusedTinyDistantSts;
576 vector<double> tinyInfoDistantD(1);//tinyInfoDistantD not used for the moment
577 DataArrayIdType *v1Distant=DataArrayIdType::New();
578 DataArrayDouble *v2Distant=DataArrayDouble::New();
579 mesh->resizeForUnserialization(tinyInfoDistant,v1Distant,v2Distant,unusedTinyDistantSts);
580 comInterface.recv(v1Distant->getPointer(),(int)v1Distant->getNbOfElems(),MPI_ID_TYPE,procId,START_TAG_MESH_XCH+2,*comm,&status);
581 comInterface.recv(v2Distant->getPointer(),(int)v2Distant->getNbOfElems(),MPI_DOUBLE,procId,START_TAG_MESH_XCH+3,*comm,&status);
582 mesh->unserialization(tinyInfoDistantD,tinyInfoDistant,v1Distant,v2Distant,unusedTinyDistantSts);
583 //finished for mesh, ids now
584 comInterface.recv(ids->getPointer(),(int)lgth[1],MPI_ID_TYPE,procId,1144,*comm,&status);
586 v1Distant->decrRef();
587 v2Distant->decrRef();