src/ParaMEDMEM/ElementLocator.cxx

   1 //  Copyright (C) 2007-2008  CEA/DEN, EDF R&D
   2 //
   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.
   7 //
   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.
  12 //
  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
  16 //
  17 //  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 #include <mpi.h>
  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"
  27
  28 #include <map>
  29 #include <set>
  30 #include <limits>
  31
  32 using namespace std;
  33
  34 namespace ParaMEDMEM
  35 {
  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)
  43   {
  44     _union_group = _local_group.fuse(distant_group);
  45     _computeBoundingBoxes();
  46   }
  47
  48   ElementLocator::~ElementLocator()
  49   {
  50     delete _union_group;
  51     delete [] _domain_bounding_boxes;
  52   }
  53
  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
  58   //  the distant mesh
  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,
  64                                     int*& distant_ids)
  65   {
  66     int dim  = _local_cell_mesh->getSpaceDimension();
  67     int rank = _union_group->translateRank(&_distant_group,idistantrank);
  68
  69     if (find(_distant_proc_ids.begin(), _distant_proc_ids.end(),rank)==_distant_proc_ids.end())
  70       {
  71         return;
  72       }
  73
  74     set <int> elems;
  75     double* distant_bb =  _domain_bounding_boxes+rank*2*dim;
  76     double* elem_bb=new double[2*dim];
  77
  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() ;
  82
  83     for ( int ielem=0; ielem<_local_cell_mesh->getNumberOfCells() ; ielem++)
  84       {
  85         for (int i=0; i<dim; i++)
  86           {
  87             elem_bb[i*2]=std::numeric_limits<double>::max();
  88             elem_bb[i*2+1]=-std::numeric_limits<double>::max();
  89           }
  90
  91         for (int inode=conn_index[ielem]+1; inode<conn_index[ielem+1]; inode++)//+1 due to offset of cell type.
  92           {
  93             int node= conn[inode];
  94
  95             for (int idim=0; idim<dim; idim++)
  96               {
  97                 if ( coords[node*dim+idim] < elem_bb[idim*2] )
  98                   {
  99                     elem_bb[idim*2] = coords[node*dim+idim] ;
 100                   }
 101                 if ( coords[node*dim+idim] > elem_bb[idim*2+1] )
 102                   {
 103                     elem_bb[idim*2+1] = coords[node*dim+idim] ;
 104                   }
 105               }
 106           }
 107         if (_intersectsBoundingBox(elem_bb, distant_bb, dim))
 108           {
 109             elems.insert(ielem);
 110           }
 111       }
 112     //send_mesh contains null pointer if elems is empty
 113     MEDCouplingUMesh* send_mesh = _meshFromElems(elems);
 114
 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
 119
 120     int* distant_ids_send=0;
 121     if (elems.size()>0)
 122       {
 123         distant_ids_send = new int[elems.size()];
 124         int index=0;
 125         for (std::set<int>::const_iterator iter = elems.begin(); iter!= elems.end(); iter++)
 126           {
 127             distant_ids_send[index]=*iter;
 128             index++;
 129           }
 130       }
 131     _exchangeMesh(send_mesh, distant_mesh, idistantrank, distant_ids_send, distant_ids);
 132     delete[] distant_ids_send;
 133     delete[] elem_bb;
 134     send_mesh->decrRef();
 135   }
 136
 137   void ElementLocator::exchangeMethod(const std::string& sourceMeth, int idistantrank, std::string& targetMeth)
 138   {
 139     CommInterface comm_interface=_union_group->getCommInterface();
 140     MPIProcessorGroup* group=static_cast<MPIProcessorGroup*> (_union_group);
 141     const MPI_Comm* comm=(group->getComm());
 142     MPI_Status status;
 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,
 150                             *comm, &status);
 151     targetMeth=recv_buffer;
 152     delete [] recv_buffer;
 153   }
 154
 155
 156   // ======================
 157   // Compute bounding boxes
 158   // ======================
 159
 160   void ElementLocator::_computeBoundingBoxes()
 161   {
 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() ;
 166
 167     int nbnodes =  _local_cell_mesh->getNumberOfNodes();
 168     double * minmax=new double [2*dim];
 169     for (int idim=0; idim<dim; idim++)
 170       {
 171         minmax[idim*2]=std::numeric_limits<double>::max();
 172         minmax[idim*2+1]=-std::numeric_limits<double>::max();
 173       }
 174
 175     for (int i=0; i<nbnodes; i++)
 176       {
 177         for (int idim=0; idim<dim;idim++)
 178           {
 179             if ( minmax[idim*2] > coords[i*dim+idim] )
 180               {
 181                 minmax[idim*2] = coords[i*dim+idim] ;
 182               }
 183             if ( minmax[idim*2+1] < coords[i*dim+idim] )
 184               {
 185                 minmax[idim*2+1] = coords[i*dim+idim] ;
 186               }
 187           }
 188       }
 189
 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,
 194                              *comm);
 195
 196     for (int i=0; i< _distant_group.size(); i++)
 197       {
 198         int rank= _union_group->translateRank(&_distant_group,i);
 199
 200         if (_intersectsBoundingBox(rank))
 201           {
 202             _distant_proc_ids.push_back(rank);
 203           }
 204       }
 205     delete[] minmax;
 206   }
 207
 208
 209   // =============================================
 210   // Intersect Bounding Box (with a given "irank")
 211   // =============================================
 212   bool ElementLocator::_intersectsBoundingBox(int irank)
 213   {
 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;
 217
 218     for (int idim=0; idim < _local_cell_mesh->getSpaceDimension(); idim++)
 219       {
 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;
 224       }
 225     return true;
 226   }
 227
 228   // =============================================
 229   // Intersect Bounding Box given 2 Bounding Boxes
 230   // =============================================
 231   bool ElementLocator::_intersectsBoundingBox(double* bb1, double* bb2, int dim)
 232   {
 233     double bbtemp[2*dim];
 234     double deltamax=0.0;
 235     double adjustment_eps=getBoundingBoxAdjustment();
 236
 237     for (int i=0; i< dim; i++)
 238       {
 239         double delta = bb1[2*i+1]-bb1[2*i];
 240         if ( delta > deltamax )
 241           {
 242             deltamax = delta ;
 243           }
 244         //    deltamax = (delta>deltamax)?delta:deltamax;
 245       }
 246     for (int i=0; i<dim; i++)
 247       {
 248         bbtemp[i*2]=bb1[i*2]-deltamax*adjustment_eps;
 249         bbtemp[i*2+1]=bb1[i*2+1]+deltamax*adjustment_eps;
 250       }
 251
 252     for (int idim=0; idim < dim; idim++)
 253       {
 254         bool intersects = (bbtemp[idim*2]<bb2[idim*2+1])
 255           && (bb2[idim*2]<bbtemp[idim*2+1]) ;
 256         if (!intersects) return false;
 257       }
 258     return true;
 259   }
 260
 261
 262   // ======================
 263   // Exchanging meshes data
 264   // ======================
 265   void ElementLocator::_exchangeMesh( MEDCouplingUMesh* local_mesh,
 266                                       MEDCouplingUMesh*& distant_mesh,
 267                                       int iproc_distant,
 268                                       const int* distant_ids_send,
 269                                       int*& distant_ids_recv)
 270   {
 271     CommInterface comm_interface=_union_group->getCommInterface();
 272
 273     // First stage : exchanging sizes
 274     // ------------------------------
 275
 276     int* send_buffer = new int[5];
 277     int* recv_buffer = new int[5];
 278
 279     //treatment for non-empty mesh
 280     int nbconn=0;
 281     int nbelems=0;
 282
 283     if (local_mesh !=0)
 284       {
 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;
 292       }
 293     else
 294       {
 295         for (int i=0; i<5; i++)
 296           {
 297             send_buffer[i]=0;
 298           }
 299       }
 300
 301     MPIProcessorGroup* group=static_cast<MPIProcessorGroup*> (_union_group);
 302     const MPI_Comm* comm=(group->getComm());
 303     MPI_Status status;
 304
 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,
 308                                                      iproc_distant);
 309
 310     comm_interface.sendRecv(send_buffer, 5, MPI_INT, iprocdistant_in_union, 1112,
 311                             recv_buffer, 5, MPI_INT,iprocdistant_in_union,1112,
 312                             *comm, &status);
 313
 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];
 319
 320     delete[] send_buffer;
 321     delete[] recv_buffer;
 322
 323     // Second stage : exchanging connectivity buffers
 324     // ----------------------------------------------
 325
 326     int nb_integers = nbconn + 2*nbelems + 1;
 327     send_buffer     = new int[nb_integers];
 328     const int* conn = 0;
 329     const int* global_numbering=0;
 330     int* ptr_buffer = send_buffer;
 331
 332     if (local_mesh != 0)
 333       {
 334         conn = local_mesh->getNodalConnectivity()->getPointer();
 335
 336         global_numbering = local_mesh->getNodalConnectivityIndex()->getPointer() ;
 337
 338         //copying the data in the integer buffer
 339
 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));
 345       }
 346
 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];
 350
 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,
 356                             *comm, &status);
 357
 358     if ( nb_integers>0 )
 359       {
 360         delete[] send_buffer;
 361       }
 362
 363     // Third stage : exchanging coordinates
 364     // ------------------------------------
 365
 366     int nb_recv_floats = distant_space_dim*distant_nnodes;
 367     int nb_send_floats = 0;
 368     double* coords=0;
 369
 370     if ( local_mesh!=0 )
 371       {
 372         nb_send_floats = local_mesh->getSpaceDimension()
 373           * local_mesh->getNumberOfNodes();
 374         coords = local_mesh->getCoords()->getPointer();
 375       }
 376
 377     DataArrayDouble* myCoords=DataArrayDouble::New();
 378     myCoords->alloc(distant_nnodes,distant_space_dim);
 379
 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);
 385
 386
 387     // Reconstructing an image of the distant mesh locally
 388
 389     if ( nb_recv_integers>0 && distant_space_dim !=0 )
 390       {
 391         MEDCouplingUMesh* meshing = MEDCouplingUMesh::New() ;
 392
 393         // Coordinates
 394         meshing->setCoords(myCoords) ;
 395         myCoords->decrRef();
 396         // Connectivity
 397
 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 ;
 403
 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) ;
 409         myConnec->decrRef();
 410         myConnecIndex->decrRef();
 411
 412         // correspondence between the distant ids and the ids of
 413         // the local reconstruction
 414
 415         distant_ids_recv=new int [distant_nb_elems];
 416         for (int i=0; i<distant_nb_elems; i++)
 417           {
 418             distant_ids_recv[i]=*work++;
 419           }
 420
 421         // Mesh dimension
 422         meshing->setMeshDimension(distant_mesh_dim);
 423
 424         distant_mesh=meshing;
 425         delete[] recv_buffer;
 426       }
 427
 428   }
 429
 430
 431   // ==============
 432   // _meshFromElems
 433   // ==============
 434
 435   MEDCouplingUMesh* ElementLocator::_meshFromElems(set<int>& elems)
 436   {
 437     //returns null pointer if there are no elems in the mesh
 438     if ( elems.size()==0 ) return 0;
 439
 440     // Defining pointers
 441     const int* conn_mesh =
 442       const_cast<int*> (_local_cell_mesh->getNodalConnectivity()->getPointer());
 443
 444     const int* conn_index =
 445       const_cast<int*> (_local_cell_mesh->getNodalConnectivityIndex()->getPointer());
 446
 447     const double* coords =
 448       const_cast<double*> ( _local_cell_mesh->getCoords()->getPointer());
 449
 450     set<int> nodes;
 451     int nbconn=0;
 452     for (set<int>::const_iterator iter=elems.begin(); iter!=elems.end(); iter++)
 453       {
 454         // Conn_index : C-like Addresses
 455         for (int inode=conn_index[*iter]+1; inode<conn_index[*iter+1]; inode++)
 456           {
 457             nodes.insert(conn_mesh[inode]);
 458             nbconn++ ;
 459           }
 460       }
 461
 462     map<int,int> big2small;
 463     int i=0;
 464     for (set<int>::const_iterator iter=nodes.begin(); iter!=nodes.end(); iter++)
 465       {
 466         big2small[*iter]=i;
 467         i++;
 468       }
 469
 470     // Memory allocate
 471     DataArrayInt *conn = DataArrayInt::New() ;
 472     conn->alloc(nbconn+elems.size(),1) ;
 473     int *connPtr=conn->getPointer();
 474
 475     DataArrayInt * connIndex = DataArrayInt::New() ;
 476     connIndex->alloc(elems.size()+1,1) ;
 477     int* connIndexPtr=connIndex->getPointer();
 478
 479     DataArrayDouble *new_coords = DataArrayDouble::New() ;
 480     new_coords->alloc(nodes.size(), _local_cell_mesh->getSpaceDimension()) ;
 481     double *new_coords_ptr = new_coords->getPointer();
 482
 483     // New connectivity table
 484     int index=0;
 485     int mainIndex=0;
 486     for (set<int>::const_iterator iter=elems.begin(); iter!=elems.end(); iter++,mainIndex++)
 487       {
 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++)
 491           {
 492             connPtr[index]=big2small[conn_mesh[inode]] ; // C-like number
 493             index++;
 494           }
 495       }
 496     connIndexPtr[mainIndex]=index;
 497     // Coordinates
 498     index=0;
 499     for (set<int>::const_iterator iter=nodes.begin(); iter!=nodes.end(); iter++)
 500       {
 501         int dim = _local_cell_mesh->getSpaceDimension();
 502         for (int i=0; i<dim;i++)
 503           {
 504             new_coords_ptr[index]=coords[(*iter)*dim+i];
 505             index++;
 506           }
 507       }
 508
 509     // Initialize
 510     MEDCouplingUMesh* meshing = MEDCouplingUMesh::New() ;
 511     meshing->setCoords(new_coords) ;
 512     new_coords->decrRef();
 513     meshing->setConnectivity(conn, connIndex) ;
 514     conn->decrRef();
 515     connIndex->decrRef();
 516     meshing->setMeshDimension(_local_cell_mesh->getMeshDimension());
 517
 518     return meshing;
 519   }
 520 }