1 // Copyright (C) 2007-2008 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
20 #include "CommInterface.hxx"
21 #include "ElementLocator.hxx"
22 #include "Topology.hxx"
23 #include "BlockTopology.hxx"
24 #include "ParaMESH.hxx"
25 #include "ProcessorGroup.hxx"
26 #include "MPIProcessorGroup.hxx"
36 ElementLocator::ElementLocator(const ParaMESH& sourceMesh,
37 const ProcessorGroup& distant_group)
38 : _local_para_mesh(sourceMesh),
39 _local_cell_mesh(sourceMesh.getCellMesh()),
40 _local_face_mesh(sourceMesh.getFaceMesh()),
41 _local_group(*sourceMesh.getBlockTopology()->getProcGroup()),
42 _distant_group(distant_group)
44 _union_group = _local_group.fuse(distant_group);
45 _computeBoundingBoxes();
48 ElementLocator::~ElementLocator()
51 delete [] _domain_bounding_boxes;
54 // ==========================================================================
55 // Procedure for exchanging mesh between a distant proc and a local processor
56 // param idistantrank proc id on distant group
57 // param distant_mesh on return , points to a local reconstruction of
59 // param distant_ids on return, contains a vector defining a correspondence
60 // between the distant ids and the ids of the local reconstruction
61 // ==========================================================================
62 void ElementLocator::exchangeMesh(int idistantrank,
63 MEDCouplingUMesh*& distant_mesh,
66 int dim = _local_cell_mesh->getSpaceDimension();
67 int rank = _union_group->translateRank(&_distant_group,idistantrank);
69 if (find(_distant_proc_ids.begin(), _distant_proc_ids.end(),rank)==_distant_proc_ids.end())
75 double* distant_bb = _domain_bounding_boxes+rank*2*dim;
76 double* elem_bb=new double[2*dim];
78 //defining pointers to med
79 const int* conn = _local_cell_mesh->getNodalConnectivity()->getPointer() ;
80 const int* conn_index= _local_cell_mesh->getNodalConnectivityIndex()->getPointer();
81 const double* coords = _local_cell_mesh->getCoords()->getPointer() ;
83 for ( int ielem=0; ielem<_local_cell_mesh->getNumberOfCells() ; ielem++)
85 for (int i=0; i<dim; i++)
87 elem_bb[i*2]=std::numeric_limits<double>::max();
88 elem_bb[i*2+1]=-std::numeric_limits<double>::max();
91 for (int inode=conn_index[ielem]+1; inode<conn_index[ielem+1]; inode++)//+1 due to offset of cell type.
93 int node= conn[inode];
95 for (int idim=0; idim<dim; idim++)
97 if ( coords[node*dim+idim] < elem_bb[idim*2] )
99 elem_bb[idim*2] = coords[node*dim+idim] ;
101 if ( coords[node*dim+idim] > elem_bb[idim*2+1] )
103 elem_bb[idim*2+1] = coords[node*dim+idim] ;
107 if (_intersectsBoundingBox(elem_bb, distant_bb, dim))
112 //send_mesh contains null pointer if elems is empty
113 MEDCouplingUMesh* send_mesh = _meshFromElems(elems);
115 // Constituting an array containing the ids of the elements that are
116 // going to be sent to the distant subdomain.
117 // This array enables the correct redistribution of the data when the
118 // interpolated field is transmitted to the target array
120 int* distant_ids_send=0;
123 distant_ids_send = new int[elems.size()];
125 for (std::set<int>::const_iterator iter = elems.begin(); iter!= elems.end(); iter++)
127 distant_ids_send[index]=*iter;
131 _exchangeMesh(send_mesh, distant_mesh, idistantrank, distant_ids_send, distant_ids);
132 delete[] distant_ids_send;
134 send_mesh->decrRef();
137 void ElementLocator::exchangeMethod(const std::string& sourceMeth, int idistantrank, std::string& targetMeth)
139 CommInterface comm_interface=_union_group->getCommInterface();
140 MPIProcessorGroup* group=static_cast<MPIProcessorGroup*> (_union_group);
141 const MPI_Comm* comm=(group->getComm());
143 // it must be converted to union numbering before communication
144 int idistRankInUnion = group->translateRank(&_distant_group,idistantrank);
145 char *recv_buffer=new char[4];
146 std::vector<char> send_buffer(4);
147 std::copy(sourceMeth.begin(),sourceMeth.end(),send_buffer.begin());
148 comm_interface.sendRecv(&send_buffer[0], 4, MPI_CHAR,idistRankInUnion, 1112,
149 recv_buffer, 4, MPI_CHAR,idistRankInUnion, 1112,
151 targetMeth=recv_buffer;
152 delete [] recv_buffer;
156 // ======================
157 // Compute bounding boxes
158 // ======================
160 void ElementLocator::_computeBoundingBoxes()
162 CommInterface comm_interface =_union_group->getCommInterface();
163 int dim = _local_cell_mesh->getSpaceDimension();
164 _domain_bounding_boxes = new double[2*dim*_union_group->size()];
165 const double* coords = _local_cell_mesh->getCoords()->getPointer() ;
167 int nbnodes = _local_cell_mesh->getNumberOfNodes();
168 double * minmax=new double [2*dim];
169 for (int idim=0; idim<dim; idim++)
171 minmax[idim*2]=std::numeric_limits<double>::max();
172 minmax[idim*2+1]=-std::numeric_limits<double>::max();
175 for (int i=0; i<nbnodes; i++)
177 for (int idim=0; idim<dim;idim++)
179 if ( minmax[idim*2] > coords[i*dim+idim] )
181 minmax[idim*2] = coords[i*dim+idim] ;
183 if ( minmax[idim*2+1] < coords[i*dim+idim] )
185 minmax[idim*2+1] = coords[i*dim+idim] ;
190 MPIProcessorGroup* group=static_cast<MPIProcessorGroup*> (_union_group);
191 const MPI_Comm* comm = group->getComm();
192 comm_interface.allGather(minmax, 2*dim, MPI_DOUBLE,
193 _domain_bounding_boxes,2*dim, MPI_DOUBLE,
196 for (int i=0; i< _distant_group.size(); i++)
198 int rank= _union_group->translateRank(&_distant_group,i);
200 if (_intersectsBoundingBox(rank))
202 _distant_proc_ids.push_back(rank);
209 // =============================================
210 // Intersect Bounding Box (with a given "irank")
211 // =============================================
212 bool ElementLocator::_intersectsBoundingBox(int irank)
214 int dim=_local_cell_mesh->getSpaceDimension();
215 double* local_bb = _domain_bounding_boxes+_union_group->myRank()*2*dim;
216 double* distant_bb = _domain_bounding_boxes+irank*2*dim;
218 for (int idim=0; idim < _local_cell_mesh->getSpaceDimension(); idim++)
220 const double eps = 1e-12;
221 bool intersects = (distant_bb[idim*2]<local_bb[idim*2+1]+eps)
222 && (local_bb[idim*2]<distant_bb[idim*2+1]+eps);
223 if (!intersects) return false;
228 // =============================================
229 // Intersect Bounding Box given 2 Bounding Boxes
230 // =============================================
231 bool ElementLocator::_intersectsBoundingBox(double* bb1, double* bb2, int dim)
233 double bbtemp[2*dim];
235 double adjustment_eps=getBoundingBoxAdjustment();
237 for (int i=0; i< dim; i++)
239 double delta = bb1[2*i+1]-bb1[2*i];
240 if ( delta > deltamax )
244 // deltamax = (delta>deltamax)?delta:deltamax;
246 for (int i=0; i<dim; i++)
248 bbtemp[i*2]=bb1[i*2]-deltamax*adjustment_eps;
249 bbtemp[i*2+1]=bb1[i*2+1]+deltamax*adjustment_eps;
252 for (int idim=0; idim < dim; idim++)
254 bool intersects = (bbtemp[idim*2]<bb2[idim*2+1])
255 && (bb2[idim*2]<bbtemp[idim*2+1]) ;
256 if (!intersects) return false;
262 // ======================
263 // Exchanging meshes data
264 // ======================
265 void ElementLocator::_exchangeMesh( MEDCouplingUMesh* local_mesh,
266 MEDCouplingUMesh*& distant_mesh,
268 const int* distant_ids_send,
269 int*& distant_ids_recv)
271 CommInterface comm_interface=_union_group->getCommInterface();
273 // First stage : exchanging sizes
274 // ------------------------------
276 int* send_buffer = new int[5];
277 int* recv_buffer = new int[5];
279 //treatment for non-empty mesh
285 nbelems = local_mesh->getNumberOfCells();
286 nbconn = local_mesh->getMeshLength();
287 send_buffer[0] = local_mesh->getSpaceDimension();
288 send_buffer[1] = local_mesh->getMeshDimension();
289 send_buffer[2] = local_mesh->getNumberOfNodes();
290 send_buffer[3] = nbelems;
291 send_buffer[4] = nbconn;
295 for (int i=0; i<5; i++)
301 MPIProcessorGroup* group=static_cast<MPIProcessorGroup*> (_union_group);
302 const MPI_Comm* comm=(group->getComm());
305 // iproc_distant is the number of proc in distant group
306 // it must be converted to union numbering before communication
307 int iprocdistant_in_union = group->translateRank(&_distant_group,
310 comm_interface.sendRecv(send_buffer, 5, MPI_INT, iprocdistant_in_union, 1112,
311 recv_buffer, 5, MPI_INT,iprocdistant_in_union,1112,
314 int distant_space_dim = recv_buffer[0];
315 int distant_mesh_dim = recv_buffer[1];
316 int distant_nnodes = recv_buffer[2];
317 int distant_nb_elems = recv_buffer[3];
318 int distant_nb_conn = recv_buffer[4];
320 delete[] send_buffer;
321 delete[] recv_buffer;
323 // Second stage : exchanging connectivity buffers
324 // ----------------------------------------------
326 int nb_integers = nbconn + 2*nbelems + 1;
327 send_buffer = new int[nb_integers];
329 const int* global_numbering=0;
330 int* ptr_buffer = send_buffer;
334 conn = local_mesh->getNodalConnectivity()->getPointer();
336 global_numbering = local_mesh->getNodalConnectivityIndex()->getPointer() ;
338 //copying the data in the integer buffer
340 memcpy(ptr_buffer, global_numbering, (nbelems+1)*sizeof(int));
341 ptr_buffer += nbelems+1;
342 memcpy(ptr_buffer,conn, nbconn*sizeof(int));
343 ptr_buffer += nbconn;
344 memcpy(ptr_buffer, distant_ids_send, nbelems*sizeof(int));
347 // Preparing the recv buffers
348 int nb_recv_integers = distant_nb_conn + 2*distant_nb_elems + 1 ;
349 recv_buffer=new int[nb_recv_integers];
351 // Exchanging integer buffer
352 comm_interface.sendRecv(send_buffer, nb_integers, MPI_INT,
353 iprocdistant_in_union, 1111,
354 recv_buffer, nb_recv_integers, MPI_INT,
355 iprocdistant_in_union,1111,
360 delete[] send_buffer;
363 // Third stage : exchanging coordinates
364 // ------------------------------------
366 int nb_recv_floats = distant_space_dim*distant_nnodes;
367 int nb_send_floats = 0;
372 nb_send_floats = local_mesh->getSpaceDimension()
373 * local_mesh->getNumberOfNodes();
374 coords = local_mesh->getCoords()->getPointer();
377 DataArrayDouble* myCoords=DataArrayDouble::New();
378 myCoords->alloc(distant_nnodes,distant_space_dim);
380 comm_interface.sendRecv(coords, nb_send_floats, MPI_DOUBLE,
381 iprocdistant_in_union, 1112,
382 myCoords->getPointer(), nb_recv_floats, MPI_DOUBLE,
383 iprocdistant_in_union, 1112,
384 *group->getComm(), &status);
387 // Reconstructing an image of the distant mesh locally
389 if ( nb_recv_integers>0 && distant_space_dim !=0 )
391 MEDCouplingUMesh* meshing = MEDCouplingUMesh::New() ;
394 meshing->setCoords(myCoords) ;
398 int *work=recv_buffer;
399 DataArrayInt* myConnecIndex=DataArrayInt::New();
400 myConnecIndex->alloc(distant_nb_elems+1,1);
401 memcpy(myConnecIndex->getPointer(), work, (distant_nb_elems+1)*sizeof(int));
402 work += distant_nb_elems + 1 ;
404 DataArrayInt* myConnec=DataArrayInt::New();
405 myConnec->alloc(distant_nb_conn,1);
406 memcpy(myConnec->getPointer(), work, (distant_nb_conn)*sizeof(int));
407 work+=distant_nb_conn;
408 meshing->setConnectivity(myConnec, myConnecIndex) ;
410 myConnecIndex->decrRef();
412 // correspondence between the distant ids and the ids of
413 // the local reconstruction
415 distant_ids_recv=new int [distant_nb_elems];
416 for (int i=0; i<distant_nb_elems; i++)
418 distant_ids_recv[i]=*work++;
422 meshing->setMeshDimension(distant_mesh_dim);
424 distant_mesh=meshing;
425 delete[] recv_buffer;
435 MEDCouplingUMesh* ElementLocator::_meshFromElems(set<int>& elems)
437 //returns null pointer if there are no elems in the mesh
438 if ( elems.size()==0 ) return 0;
441 const int* conn_mesh =
442 const_cast<int*> (_local_cell_mesh->getNodalConnectivity()->getPointer());
444 const int* conn_index =
445 const_cast<int*> (_local_cell_mesh->getNodalConnectivityIndex()->getPointer());
447 const double* coords =
448 const_cast<double*> ( _local_cell_mesh->getCoords()->getPointer());
452 for (set<int>::const_iterator iter=elems.begin(); iter!=elems.end(); iter++)
454 // Conn_index : C-like Addresses
455 for (int inode=conn_index[*iter]+1; inode<conn_index[*iter+1]; inode++)
457 nodes.insert(conn_mesh[inode]);
462 map<int,int> big2small;
464 for (set<int>::const_iterator iter=nodes.begin(); iter!=nodes.end(); iter++)
471 DataArrayInt *conn = DataArrayInt::New() ;
472 conn->alloc(nbconn+elems.size(),1) ;
473 int *connPtr=conn->getPointer();
475 DataArrayInt * connIndex = DataArrayInt::New() ;
476 connIndex->alloc(elems.size()+1,1) ;
477 int* connIndexPtr=connIndex->getPointer();
479 DataArrayDouble *new_coords = DataArrayDouble::New() ;
480 new_coords->alloc(nodes.size(), _local_cell_mesh->getSpaceDimension()) ;
481 double *new_coords_ptr = new_coords->getPointer();
483 // New connectivity table
486 for (set<int>::const_iterator iter=elems.begin(); iter!=elems.end(); iter++,mainIndex++)
488 connIndexPtr[mainIndex]=index;
489 connPtr[index++]=conn_mesh[conn_index[*iter]];
490 for (int inode = conn_index[*iter]+1; inode < conn_index[*iter+1]; inode++)
492 connPtr[index]=big2small[conn_mesh[inode]] ; // C-like number
496 connIndexPtr[mainIndex]=index;
499 for (set<int>::const_iterator iter=nodes.begin(); iter!=nodes.end(); iter++)
501 int dim = _local_cell_mesh->getSpaceDimension();
502 for (int i=0; i<dim;i++)
504 new_coords_ptr[index]=coords[(*iter)*dim+i];
510 MEDCouplingUMesh* meshing = MEDCouplingUMesh::New() ;
511 meshing->setCoords(new_coords) ;
512 new_coords->decrRef();
513 meshing->setConnectivity(conn, connIndex) ;
515 connIndex->decrRef();
516 meshing->setMeshDimension(_local_cell_mesh->getMeshDimension());