src/MEDSPLITTER/MEDSPLITTER_MESHCollection.cxx

   1 // Copyright (C) 2007-2013  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
  20 #include "MEDMEM_ConnectZone.hxx"
  21 #include "MEDMEM_DriversDef.hxx"
  22 #include "MEDMEM_Mesh.hxx"
  23 #include "MEDMEM_Meshing.hxx"
  24 #include "MEDMEM_GaussLocalization.hxx"
  25 #include "MEDMEM_Field.hxx"
  26 #include "MEDMEM_CellModel.hxx"
  27 #include "MEDMEM_Group.hxx"
  28 #include "MEDMEM_MeshFuse.hxx"
  29
  30 #include "MEDMEM_Exception.hxx"
  31
  32 #include "MEDSPLITTER_utils.hxx"
  33
  34 #include "MEDSPLITTER_Graph.hxx"
  35
  36 #include "MEDSPLITTER_Topology.hxx"
  37 #include "MEDSPLITTER_ParallelTopology.hxx"
  38 #include "MEDSPLITTER_SequentialTopology.hxx"
  39 #include "MEDSPLITTER_ParaDomainSelector.hxx"
  40 #include "MEDSPLITTER_MeshSendReceive.hxx"
  41 #include "MEDSPLITTER_JointExchangeData.hxx"
  42
  43 #include "MEDSPLITTER_MESHCollection.hxx"
  44 #include "MEDSPLITTER_MESHCollectionDriver.hxx"
  45 #include "MEDSPLITTER_MESHCollectionMedXMLDriver.hxx"
  46 #include "MEDSPLITTER_MESHCollectionMedAsciiDriver.hxx"
  47
  48 #include "MEDSPLITTER_UserGraph.hxx"
  49
  50 #if defined(MED_ENABLE_METIS) || defined(MED_ENABLE_PARMETIS)
  51 #include "MEDSPLITTER_METISGraph.hxx"
  52 #endif
  53 #ifdef MED_ENABLE_SCOTCH
  54 #include "MEDSPLITTER_SCOTCHGraph.hxx"
  55 #endif
  56
  57 #include "InterpKernelHashMap.hxx"
  58
  59 #include <vector>
  60 #include <string>
  61 #include <set>
  62
  63 #include <iostream>
  64 #include <fstream>
  65
  66 using namespace MEDSPLITTER;
  67
  68 using namespace INTERP_KERNEL;
  69
  70 //template inclusion
  71 #include "MEDSPLITTER_MESHCollection.H"
  72
  73 MESHCollection::MESHCollection()
  74   : _topology(0),
  75     _owns_topology(false),
  76     _driver(0),
  77     _domain_selector( 0 ),
  78     _i_non_empty_mesh(-1),
  79     _driver_type(MEDSPLITTER::MedXML),
  80     _subdomain_boundary_creates(false),
  81     _family_splitting(false),
  82     _create_empty_groups(false)
  83 {
  84 }
  85
  86 namespace
  87 {
  88   //================================================================================
  89   /*!
  90    * \brief Creates a new domain mesh
  91    */
  92   //================================================================================
  93
  94   MEDMEM::MESH* newMesh(const std::string& name, int dim, int space, MEDMEM::MESH* meshToDelete=0)
  95   {
  96     delete meshToDelete;
  97     MEDMEM::MESHING* mesh = new MEDMEM::MeshFuse;
  98     mesh->setName( name );
  99     //mesh->setMeshDimension ( dim );
 100     //mesh->setSpaceDimension( space );
 101     return mesh;
 102   }
 103 }
 104
 105 /*!constructor creating a new mesh collection (mesh series + topology)
 106  *from an old collection and a new topology
 107  *
 108  * On output, the constructor has built the meshes corresponding to the new mesh collection.
 109  * The new topology has been updated so that face and node mappings are included.
 110  * The families have been cast to their projections in the new topology.
 111  *
 112  * \param initial_collection collection from which the data (coordinates, connectivity) are taken
 113  * \param topology topology containing the cell mappings
 114  */
 115
 116 MESHCollection::MESHCollection(const MESHCollection& initial_collection, Topology* topology, bool family_splitting, bool create_empty_groups)
 117   : _topology(topology),
 118     _owns_topology(false),
 119     _cell_graph(topology->getGraph()),
 120     _driver(0),
 121     _domain_selector( initial_collection._domain_selector ),
 122     _i_non_empty_mesh(-1),
 123     _name(initial_collection._name),
 124     _driver_type(MEDSPLITTER::MedXML),
 125     _subdomain_boundary_creates(false),
 126     _family_splitting(family_splitting),
 127     _create_empty_groups(create_empty_groups)
 128 {
 129   string mesh_name = initial_collection.getName();
 130   _mesh.resize(_topology->nbDomain());
 131
 132   int space_dim = initial_collection.getSpaceDimension();
 133   int mesh_dim  = initial_collection.getMeshDimension();
 134   if ( mesh_dim < 1 )
 135     space_dim = mesh_dim = initial_collection._driver->readMeshDimension();
 136
 137   for (int idomain=0; idomain < _topology->nbDomain(); idomain++)
 138   {
 139     //creating the new mesh
 140     _mesh[idomain]= newMesh( MEDMEM::STRING(mesh_name)<<"_"<<idomain+1, mesh_dim, space_dim );
 141
 142     createNodalConnectivity(initial_collection,idomain, MED_EN::MED_CELL);
 143
 144     if ( _mesh[idomain]->getNumberOfNodes() > 0 )
 145       _i_non_empty_mesh = idomain;
 146   }
 147
 148   _topology->createFaceMapping(initial_collection, *this );
 149   for (int idomain=0; idomain<_topology->nbDomain(); idomain++)
 150   {
 151     switch (mesh_dim)
 152     {
 153     case 3:
 154       createNodalConnectivity(initial_collection,idomain, MED_EN::MED_FACE);
 155       break;
 156     case 2:
 157       createNodalConnectivity(initial_collection,idomain, MED_EN::MED_EDGE);
 158       break;
 159     default :
 160       if ( !isParallelMode() || _domain_selector->isMyDomain( idomain ))
 161         cerr<<"MEDSPLITTER : Mesh dimension must be 2 or 3"<<endl;
 162     }
 163   }
 164
 165   castFamilies(initial_collection);
 166
 167   // Exchange domain parts
 168
 169   if ( isParallelMode() )
 170   {
 171     _domain_selector->setNbDomains( _topology->nbDomain() );
 172
 173     vector< MeshSendReceive > mesh_sender( _topology->nbDomain() );
 174     list<int> domainsToClear; // sent domains
 175     bool isSent;
 176     // first, send domains
 177     for (int idomain=0; idomain<_topology->nbDomain(); idomain++)
 178     {
 179       // get node numbers global over all procs
 180       vector<int> glob_nodes_all_proc( _topology->getNodeNumber( idomain )); // to fill in
 181       vector<int> glob_cells_all_proc( _topology->getCellNumber( idomain ));
 182       vector<int> glob_faces_all_proc( _topology->getFaceNumber( idomain ));
 183       if ( !glob_cells_all_proc.empty() )
 184       {
 185         // get ids global on this proc
 186         _topology->getNodeList( idomain, & glob_nodes_all_proc[0] );
 187         _topology->getCellList( idomain, & glob_cells_all_proc[0] );
 188         _topology->getFaceList( idomain, & glob_faces_all_proc[0] );
 189         // convert cell ids to ids global over all procs
 190         int cell_shift = _domain_selector->getProcShift();
 191         for ( int i = 0; i < glob_cells_all_proc.size(); ++i )
 192           glob_cells_all_proc[i] += cell_shift;
 193       }
 194       if ( _domain_selector->isMyDomain( idomain ))
 195       {
 196         // prepare to receiving other parts of the domain
 197         ((MEDMEM::MeshFuse*) _mesh[idomain])->setNodeNumbers( glob_nodes_all_proc );
 198         _topology->getFusedCellNumbers( idomain ) = glob_cells_all_proc;
 199         _topology->getFusedFaceNumbers( idomain ) = glob_faces_all_proc;
 200       }
 201       else
 202       {
 203         // sending
 204         int target_proc = _domain_selector->getProccessorID( idomain );
 205         mesh_sender[ idomain ].send( target_proc, idomain, _mesh[idomain],
 206                                      glob_cells_all_proc,
 207                                      glob_faces_all_proc,
 208                                      glob_nodes_all_proc );
 209         if ( !glob_nodes_all_proc.empty() )
 210           domainsToClear.push_back( idomain );
 211       }
 212       // clear just sent domain meshes
 213       for ( list<int>::iterator dom = domainsToClear.begin(); dom != domainsToClear.end(); )
 214       {
 215         if (( isSent = mesh_sender[ *dom ].isSent() ))
 216           _mesh[*dom] = newMesh( _mesh[*dom]->getName(), mesh_dim, space_dim, _mesh[*dom]);
 217         dom = isSent ? domainsToClear.erase( dom ) : ++dom;
 218       }
 219     }
 220
 221     // then, receive domains
 222     MeshSendReceive mesh_receiver;
 223     int this_proc = _domain_selector->rank();
 224     for (int idomain=0; idomain<_topology->nbDomain(); idomain++)
 225     {
 226       if ( _domain_selector->isMyDomain( idomain ))
 227       {
 228         for (int iproc = 0; iproc < _domain_selector->nbProcs(); ++iproc)
 229         {
 230           if ( iproc == this_proc ) continue;
 231           vector<int>  nodes_other_proc, cells_other_proc, faces_other_proc;
 232           mesh_receiver.recv( iproc, idomain, cells_other_proc, faces_other_proc,nodes_other_proc);
 233           if ( MEDMEM::MESH* received_mesh = mesh_receiver.getMesh() )
 234           {
 235             // unite meshes and global node numbers stored in MeshFuse
 236             MEDMEM::MeshFuse* fuse = (MEDMEM::MeshFuse*) _mesh[idomain];
 237             fuse->concatenate( received_mesh, nodes_other_proc );
 238             delete received_mesh;
 239
 240             // unite global element numbers
 241             fuse->append( MED_EN::MED_CELL,
 242                           _topology->getFusedCellNumbers( idomain ), cells_other_proc );
 243
 244             fuse->append( mesh_dim==3 ? MED_EN::MED_FACE : MED_EN::MED_EDGE,
 245                           _topology->getFusedFaceNumbers( idomain ), faces_other_proc );
 246
 247             if ( _mesh[idomain]->getNumberOfNodes() > 0 )
 248               _i_non_empty_mesh = idomain;
 249           }
 250         }
 251       }
 252       // clear just sent domain meshes
 253       for ( list<int>::iterator dom = domainsToClear.begin(); dom != domainsToClear.end(); )
 254       {
 255         if (( isSent = mesh_sender[ *dom ].isSent() ))
 256           _mesh[*dom] = newMesh( _mesh[*dom]->getName(), mesh_dim, space_dim,_mesh[*dom]);
 257         dom = isSent ? domainsToClear.erase( dom ) : ++dom;
 258       }
 259     }
 260     // clear sent domain meshes
 261     mesh_sender.clear();
 262     for ( list<int>::iterator dom = domainsToClear.begin(); dom != domainsToClear.end(); ++dom)
 263       _mesh[*dom] = newMesh( _mesh[*dom]->getName(), mesh_dim, space_dim,_mesh[*dom]);
 264
 265     _topology->recreateMappingAfterFusion( getMesh() );
 266   }
 267   if ( _i_non_empty_mesh < 0 ) // non of domains resides on this proc,
 268     _i_non_empty_mesh = 0; // in this case we need only dimension that is set to all meshes
 269
 270 }
 271
 272 /*! constructing the MESH collection from a distributed file
 273  *
 274  * \param filename name of the master file containing the list of all the MED files
 275  */
 276 MESHCollection::MESHCollection(const string& filename)
 277   : _topology(0),
 278     _owns_topology(true),
 279     _driver(0),
 280     _domain_selector( 0 ),
 281     _i_non_empty_mesh(-1),
 282     _driver_type(MEDSPLITTER::Undefined),
 283     _subdomain_boundary_creates(false),
 284     _family_splitting(false),
 285     _create_empty_groups(false)
 286 {
 287   char filenamechar[256];
 288   strcpy(filenamechar,filename.c_str());
 289   try
 290   {
 291     _driver=new MESHCollectionMedXMLDriver(this);
 292     _driver->read (filenamechar);
 293     _driver_type = MedXML;
 294
 295   }
 296   catch(MEDEXCEPTION&){
 297     delete _driver;
 298     try
 299     {
 300       _driver=new MESHCollectionMedAsciiDriver(this);
 301       _driver->read (filenamechar);
 302       _driver_type=MedAscii;
 303     }
 304     catch(MEDEXCEPTION&)
 305     {
 306       delete _driver;
 307       throw MEDEXCEPTION("file does not comply with any recognized format");
 308     }
 309   }
 310   for ( int idomain = 0; idomain < _mesh.size(); ++idomain )
 311     if ( _mesh[idomain] && _mesh[idomain]->getNumberOfNodes() > 0 )
 312       _i_non_empty_mesh = idomain;
 313 }
 314
 315 /*! Constructing the MESH collection from selected domains of a distributed file
 316  *
 317  * \param filename  - name of the master file containing the list of all the MED files
 318  * \param domainSelector - selector of domains to load
 319  */
 320 MESHCollection::MESHCollection(const string& filename, ParaDomainSelector& domainSelector)
 321   : _topology(0),
 322     _owns_topology(true),
 323     _driver(0),
 324     _domain_selector( domainSelector.nbProcs() > 1 ? & domainSelector : 0 ),
 325     _i_non_empty_mesh(-1),
 326     _driver_type(MEDSPLITTER::Undefined),
 327     _subdomain_boundary_creates(false),
 328     _family_splitting(false),
 329     _create_empty_groups(false)
 330 {
 331   try
 332   {
 333     _driver=new MESHCollectionMedXMLDriver(this);
 334     _driver->read ( (char*)filename.c_str(), _domain_selector );
 335     _driver_type = MedXML;
 336   }
 337   catch(MEDEXCEPTION&)
 338   {
 339     delete _driver;
 340     try
 341     {
 342       _driver=new MESHCollectionMedAsciiDriver(this);
 343       _driver->read ( (char*)filename.c_str(), _domain_selector );
 344       _driver_type=MedAscii;
 345     }
 346     catch(MEDEXCEPTION&)
 347     {
 348       delete _driver;
 349       throw MEDEXCEPTION("file does not comply with any recognized format");
 350     }
 351   }
 352   if ( isParallelMode() )
 353     // to know nb of cells on each proc to compute global cell ids from locally global
 354     _domain_selector->gatherNbOf( MED_EN::MED_CELL, getMesh() );
 355
 356   // find non-empty domain mesh
 357   for ( int idomain = 0; idomain < _mesh.size(); ++idomain )
 358     if ( _mesh[idomain] && _mesh[idomain]->getNumberOfNodes() > 0 )
 359       _i_non_empty_mesh = idomain;
 360 }
 361
 362 /*! constructing the MESH collection from a sequential MED-file
 363  *
 364  * \param filename MED file
 365  * \param meshname name of the mesh that is to be read
 366  */
 367 MESHCollection::MESHCollection(const string& filename, const string& meshname)
 368   : _topology(0),
 369     _owns_topology(true),
 370     _driver(0),
 371     _domain_selector( 0 ),
 372     _i_non_empty_mesh(-1),
 373     _name(meshname),
 374     _driver_type(MEDSPLITTER::MedXML),
 375     _subdomain_boundary_creates(false),
 376     _family_splitting(false),
 377     _create_empty_groups(false)
 378 {
 379   char filenamechar[256];
 380   char meshnamechar[256];
 381   strcpy(filenamechar,filename.c_str());
 382   strcpy(meshnamechar,meshname.c_str());
 383   try // avoid memory leak in case of inexistent filename
 384   {
 385     retrieveDriver()->readSeq (filenamechar,meshnamechar);
 386   }
 387   catch ( MED_EXCEPTION& e )
 388   {
 389     if ( _driver ) delete _driver; _driver=0;
 390     throw e;
 391   }
 392   if ( _mesh[0] && _mesh[0]->getNumberOfNodes() > 0 )
 393     _i_non_empty_mesh = 0;
 394 }
 395
 396 MESHCollection::~MESHCollection()
 397 {
 398   for (unsigned i=0; i<_mesh.size();i++)
 399     if (_mesh[i]!=0) {/*delete*/ _mesh[i]->removeReference(); }
 400   for (unsigned i=0; i<_connect_zones.size();i++)
 401     if (_connect_zones[i]!=0) {delete _connect_zones[i];}
 402   if (_driver !=0) {delete _driver; _driver=0;}
 403   if (_topology!=0 && _owns_topology) {delete _topology; _topology=0;}
 404 }
 405
 406 /*!gets the connectivity for a certain type
 407  *
 408  * The output array type_connectivity should have been allocated
 409  * at dimension nbnode_per_type* nb_cells before the call
 410  *
 411  * \param cell_list list of elements (global cell numbers) for which the connectivity is required
 412  * \param nb_cells number of elements
 413  * \param entity type of entity for which the nodal connectivity is required
 414  * \param type type of the elements for which the connectivity is required
 415  * \param type_connectivity on output contains the connectivity of all the elements of the list
 416  * */
 417
 418 void MESHCollection::getNodeConnectivity(const int* cell_list,int nb_cells,MED_EN::medEntityMesh entity,
 419                                          MED_EN::medGeometryElement type, int* type_connectivity) const
 420 {
 421   int *local=new int[nb_cells];
 422   int *ip=new int[nb_cells];
 423   switch (entity)
 424   {
 425   case MED_EN::MED_CELL:
 426     _topology->convertGlobalCellList(cell_list,nb_cells,local,ip);
 427     break;
 428   case MED_EN::MED_FACE:
 429   case MED_EN::MED_EDGE:
 430     _topology->convertGlobalFaceList(cell_list,nb_cells,local,ip);
 431     break;
 432   }
 433
 434
 435   //  int nbnode_per_type=(int)type%100;
 436   //  vector<int> number_of_types_array(_topology->nbDomain(),0);
 437   //  for (int i=0; i<_topology->nbDomain(); i++)
 438   //    number_of_types_array[i]=_mesh[i]->getNumberOfTypes(entity);
 439
 440   //defining a connectivity table for different domains
 441   vector  <const int*> conn_ip(_topology->nbDomain());
 442   vector  <const int*> conn_index_ip(_topology->nbDomain());
 443
 444
 445   vector< map <MED_EN::medGeometryElement, int> > offset;
 446   //  offset.resize(_topology->nbDomain());
 447
 448   for (int i=0; i<_topology->nbDomain();i++)
 449   {
 450     if ( !_mesh[i] ) continue;
 451     int nb_elem = _mesh[i]->getNumberOfElements(entity,type);
 452     if (nb_elem>0)
 453     {
 454       conn_ip[i]=_mesh[i]->getConnectivity(MED_EN::MED_NODAL,entity,MED_EN::MED_ALL_ELEMENTS);
 455       conn_index_ip[i] = _mesh[i]->getConnectivityIndex(MED_EN::MED_NODAL,entity);
 456       //       global_index= _mesh[i]->getGlobalNumberingIndex(entity);
 457     }
 458     else
 459     {
 460       conn_ip[i]=0;
 461       conn_index_ip[i]=0;
 462     }
 463     //      int number_of_types = number_of_types_array[i];
 464     //      const MEDMEM::CELLMODEL* types =  _mesh[ip[icell]]->getCellsTypes(entity);
 465     //      for (int itype=0; itype<number_of_types; itype++)
 466     //        offset[i][types[itype].getType()]=global_index[itype]-1;
 467   }
 468
 469   int* type_connectivity_ptr=type_connectivity;
 470   for (int icell=0; icell<nb_cells; icell++)
 471   {
 472     //      int type_offset = offset[ip[icell]][type];
 473     const int* conn=conn_ip[ip[icell]];
 474     const int* conn_index=conn_index_ip[ip[icell]];
 475     for (int inode=conn_index[local[icell]-1]; inode<conn_index[local[icell]]; inode++)
 476     {
 477       *type_connectivity_ptr=
 478         _topology->convertNodeToGlobal(ip[icell],conn[inode-1]);
 479       type_connectivity_ptr++;
 480     }
 481   }
 482
 483   delete[]local;
 484   delete[]ip;
 485 }
 486
 487 /*!gets the connectivity for MED_POLYGON type
 488  *
 489  * \param cell_list list of elements (global cell numbers) for which the connectivity is required
 490  * \param nb_cells number of elements
 491  * \param entity type of entity for which the nodal connectivity is required
 492  * \param type_connectivity on output contains the connectivity of all the elements of the list
 493  * \param connectivity_index on output contains the connectivity index for all polygons
 494  * */
 495
 496 void MESHCollection::getPolygonNodeConnectivity(const int* cell_list,int nb_cells,MED_EN::medEntityMesh entity,
 497                                                 vector<int>& type_connectivity, vector<int>& connectivity_index) const
 498 {
 499
 500   int *local=new int[nb_cells];
 501   int *ip=new int[nb_cells];
 502   switch (entity)
 503   {
 504   case MED_EN::MED_CELL:
 505     _topology->convertGlobalCellList(cell_list,nb_cells,local,ip);
 506     break;
 507   case MED_EN::MED_FACE:
 508   case MED_EN::MED_EDGE:
 509     _topology->convertGlobalFaceList(cell_list,nb_cells,local,ip);
 510     break;
 511   }
 512
 513
 514   //defining a connectivity table for different domains
 515   vector  <const int*> conn_ip(_topology->nbDomain());
 516   vector  <const int*> conn_index_ip(_topology->nbDomain());
 517   vector <const int* > conn_face_index(_topology->nbDomain());
 518   vector<int> nb_plain_elems(_topology->nbDomain());
 519
 520   vector< map <MED_EN::medGeometryElement, int> > offset;
 521
 522   for (int i=0; i<_topology->nbDomain();i++)
 523   {
 524     int nb_elem = _mesh[i]->getNumberOfElements(entity,MED_EN::MED_POLYGON);
 525     if (nb_elem>0)
 526     {
 527       conn_ip[i]=_mesh[i]->getConnectivity(MED_EN::MED_NODAL,entity,MED_EN::MED_ALL_ELEMENTS);
 528       conn_index_ip[i] = _mesh[i]->getConnectivityIndex(MED_EN::MED_NODAL,entity);
 529     }
 530     else
 531     {
 532       conn_ip[i]=0;
 533       conn_index_ip[i]=0;
 534     }
 535   }
 536
 537   connectivity_index.resize(nb_cells+1);
 538   connectivity_index[0]=1;
 539   for (int icell=0; icell<nb_cells; icell++)
 540   {
 541     //int nb_plain= nb_plain_elems[ip[icell]];
 542     const int* conn=conn_ip[ip[icell]];
 543     const int* conn_index=conn_index_ip[ip[icell]];
 544     for (int inode=conn_index[local[icell]-1/*-nb_plain*/]; inode<conn_index[local[icell]/*-nb_plain*/]; inode++)
 545     {
 546       type_connectivity.push_back(
 547                                   _topology->convertNodeToGlobal(ip[icell],conn[inode-1]));
 548     }
 549     connectivity_index[icell+1]=connectivity_index[icell]
 550       -conn_index[local[icell]-1/*-nb_plain*/]+conn_index[local[icell]/*-nb_plain*/];
 551   }
 552
 553   delete[]local;
 554   delete[]ip;
 555 }
 556
 557
 558 /*!gets the connectivity for MED_POLYHEDRA type
 559  *
 560  * \param cell_list list of elements (global cell numbers) for which the connectivity is required
 561  * \param nb_cells number of elements
 562  * \param entity type of entity for which the nodal connectivity is required
 563  * \param type_connectivity on output contains the connectivity of all the elements of the list
 564  * \param connectivity_index on output contains the connectivity index for all polygons
 565  * */
 566
 567 void MESHCollection::getPolyhedraNodeConnectivity(const int* cell_list,int nb_cells,MED_EN::medEntityMesh entity,
 568                                                   vector<int>& type_connectivity, vector<int>& connectivity_index/*, vector<int>& face_connectivity_index*/) const
 569 {
 570
 571   int *local=new int[nb_cells];
 572   int *ip=new int[nb_cells];
 573   switch (entity)
 574   {
 575   case MED_EN::MED_CELL:
 576     _topology->convertGlobalCellList(cell_list,nb_cells,local,ip);
 577     break;
 578   case MED_EN::MED_FACE:
 579   case MED_EN::MED_EDGE:
 580     _topology->convertGlobalFaceList(cell_list,nb_cells,local,ip);
 581     break;
 582   }
 583
 584
 585   //defining a connectivity table for different domains
 586   vector  <const int*> conn_ip(_topology->nbDomain());
 587   vector  <const int*> conn_index_ip(_topology->nbDomain());
 588   vector<int> nb_plain_elems(_topology->nbDomain());
 589
 590   vector< map <MED_EN::medGeometryElement, int> > offset;
 591
 592   for (int i=0; i<_topology->nbDomain();i++)
 593   {
 594     nb_plain_elems[i] = _mesh[i]->getNumberOfElements(entity, MED_EN::MED_ALL_ELEMENTS);
 595     int nb_elem = _mesh[i]->getNumberOfElements(entity,MED_EN::MED_POLYHEDRA);
 596     if (nb_elem>0)
 597     {
 598       conn_ip[i]=_mesh[i]->getConnectivity(MED_EN::MED_NODAL,MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
 599       conn_index_ip[i] = _mesh[i]->getConnectivityIndex(MED_EN::MED_NODAL,MED_EN::MED_CELL);
 600     }
 601     else
 602     {
 603       conn_ip[i]=0;
 604       conn_index_ip[i]=0;
 605     }
 606   }
 607
 608   connectivity_index.resize(nb_cells+1);
 609   connectivity_index[0]=1;
 610   for (int icell=0; icell<nb_cells; icell++)
 611   {
 612     const int* conn=conn_ip[ip[icell]];
 613     const int* conn_index=conn_index_ip[ip[icell]];
 614     connectivity_index[icell+1]=connectivity_index[icell]+
 615       conn_index[local[icell]]-conn_index[local[icell]-1];
 616
 617     for (int inode=conn_index[local[icell]-1]; inode<conn_index[local[icell]]; inode++)
 618     {
 619       if ( conn[inode-1] == -1 )
 620         type_connectivity.push_back( -1 );
 621       else
 622         type_connectivity.push_back(_topology->convertNodeToGlobal(ip[icell],conn[inode-1]));
 623     }
 624
 625   }
 626
 627   delete[]local;
 628   delete[]ip;
 629 }
 630
 631 /*! constructing the MESH collection from a file
 632  *
 633  * The method creates as many MED-files as there are domains in the
 634  * collection. It also creates a master file that lists all the MED files.
 635  * The MED files created in ths manner contain joints that describe the
 636  * connectivity between subdomains.
 637  *
 638  * \param filename name of the master file that will contain the list of the MED files
 639  *
 640  */
 641 void MESHCollection::write(const string& filename)
 642 {
 643   //building the connect zones necessary for writing joints
 644   cout<<"Building Connect Zones"<<endl;
 645   if (_topology->nbDomain()>1)
 646     buildConnectZones();
 647   cout <<"End of connect zones building"<<endl;
 648   //suppresses link with driver so that it can be changed for writing
 649   if (_driver!=0)delete _driver;
 650   _driver=0;
 651
 652   char filenamechar[256];
 653   strcpy(filenamechar,filename.c_str());
 654   retrieveDriver()->write (filenamechar, _domain_selector);
 655 }
 656
 657 /*! creates or gets the link to the collection driver
 658  */
 659 MESHCollectionDriver* MESHCollection::retrieveDriver()
 660 {
 661   if (_driver==0)
 662   {
 663     switch(_driver_type)
 664     {
 665     case MedXML:
 666       _driver=new MESHCollectionMedXMLDriver(this);
 667       break;
 668     case MedAscii:
 669       _driver=new MESHCollectionMedAsciiDriver(this);
 670       break;
 671     default:
 672       throw MEDEXCEPTION("Unrecognized driver");
 673     }
 674   }
 675
 676   return _driver;
 677 }
 678
 679
 680 /*! gets an existing driver
 681  *
 682  */
 683 MESHCollectionDriver* MESHCollection::getDriver() const
 684 {
 685   return _driver;
 686 }
 687
 688
 689 /*! gets the list of types for global numbers cell_list
 690  *
 691  * \param cell_list list of global numbers
 692  * \param entity entity type
 693  * \param type_list on output, list of types for the cells given in cell_list
 694  */
 695 void MESHCollection::getTypeList(int* cell_list,int nb_cells,
 696                                  MED_EN::medEntityMesh entity,
 697                                  MED_EN::medGeometryElement* type_list) const
 698 {
 699   MESSAGE_MED (" Beginning of getTypeList with entity "<<entity);
 700   int *local=new int[nb_cells];
 701   int *ip=new int[nb_cells];
 702   switch (entity)
 703   {
 704   case MED_EN::MED_CELL:
 705     _topology->convertGlobalCellList(cell_list,nb_cells,local,ip);
 706     break;
 707   case MED_EN::MED_FACE:
 708   case MED_EN::MED_EDGE:
 709     _topology->convertGlobalFaceList(cell_list,nb_cells,local,ip);
 710     break;
 711   }
 712
 713   for (int icell=0; icell<nb_cells; icell++)
 714   {
 715     type_list[icell]=_mesh[ip[icell]]->getElementType(entity,local[icell]);
 716   }
 717   delete[]local;
 718   delete[]ip;
 719   MESSAGE_MED("end of getTypeList");
 720 }
 721
 722
 723
 724 /*!gets the descending connectivity for a certain type
 725  *
 726  * The output array type_connectivity should have been allocated
 727  * at dimension nbnode_per_type* nb_cells before the call
 728  *
 729  * \param cell_list list of elements (global cell numbers) for which the connectivity is required
 730  * \param nb_cells number of elements
 731  * \param entity type of entity for which the nodal connectivity is required
 732  * \param type type of the elements for which the connectivity is required
 733  * \param type_connectivity on output contains the connectivity of all the elements of the list
 734  * */
 735
 736 void MESHCollection::getFaceConnectivity(const int* cell_list,int nb_cells,MED_EN::medEntityMesh entity,
 737                                          MED_EN::medGeometryElement type, int* type_connectivity) const
 738 {
 739   int *local=new int[nb_cells];
 740   int *ip=new int[nb_cells];
 741   switch (entity)
 742   {
 743   case MED_EN::MED_CELL:
 744     _topology->convertGlobalCellList(cell_list,nb_cells,local,ip);
 745     break;
 746   case MED_EN::MED_FACE:
 747   case MED_EN::MED_EDGE:
 748     _topology->convertGlobalFaceList(cell_list,nb_cells,local,ip);
 749     break;
 750   }
 751
 752
 753   int nbface_per_type;
 754   switch (type){
 755   case 308:
 756     nbface_per_type=6;
 757     break;
 758   case 304:
 759     nbface_per_type=4;
 760     break;
 761   case 306:
 762     nbface_per_type=5;
 763     break;
 764   }
 765
 766   vector<int> number_of_types_array(_topology->nbDomain(),0);
 767   for (int i=0; i<_topology->nbDomain(); i++)
 768     number_of_types_array[i]=_mesh[i]->getNumberOfTypes(entity);
 769
 770   //defining a connectivity table for different domains
 771   vector  <const int*> conn_ip(_topology->nbDomain());
 772   for (int i=0; i<_topology->nbDomain();i++)
 773   {
 774     int nb_elem = _mesh[i]->getNumberOfElements(entity,type);
 775     if (nb_elem>0)
 776       conn_ip[i]=_mesh[i]->getConnectivity(MED_EN::MED_DESCENDING,entity,type);
 777     else
 778       conn_ip[i]=0;
 779   }
 780
 781   for (int icell=0; icell<nb_cells; icell++)
 782   {
 783     int number_of_types = number_of_types_array[ip[icell]];
 784     const MEDMEM::CELLMODEL* types =  _mesh[ip[icell]]->getCellsTypes(entity);
 785     int type_offset=0;
 786     for (int itype=0; itype< number_of_types; itype++)
 787     {
 788       if (types[itype].getType() < type)
 789         type_offset += _mesh[ip[icell]]->getNumberOfElements(entity,types[itype].getType());
 790     }
 791     const int* conn=conn_ip[ip[icell]];
 792     for (int iface=0; iface<nbface_per_type; iface++)
 793     {
 794       type_connectivity[icell*nbface_per_type+iface] = _topology->convertFaceToGlobal
 795         (ip[icell], abs(conn[(local[icell] - type_offset - 1) * nbface_per_type + iface]));
 796     }
 797   }
 798
 799   delete[]local;
 800   delete[]ip;
 801 }
 802
 803 /*! gets the list of coordinates for a given list of global node numbers
 804  *
 805  * The vector containing the coordinates on output should
 806  * have been allocated at a dimension _space_dimension * nb_nodes
 807  * before the call
 808  *
 809  * \param node_list list of global node numbers
 810  * \param nb_nodes number of nodes in the list
 811  * \param coordinates on output, contains the coordinates
 812  */
 813
 814 void MESHCollection::getCoordinates(int* node_list,int nb_nodes, double* coordinates) const
 815 {
 816   int* local=new int[nb_nodes];
 817   int* ip=new int[nb_nodes];
 818   int space_dimension= getSpaceDimension();
 819   _topology->convertGlobalNodeList(node_list,nb_nodes,local,ip);
 820   for (int i=0; i< nb_nodes; i++)
 821   {
 822     const double* coord=_mesh[ip[i]]->getCoordinates(MED_EN::MED_FULL_INTERLACE);
 823     for (int icoord=0; icoord<space_dimension; icoord++)
 824       coordinates[i*space_dimension+icoord]=coord[(local[i]-1)*space_dimension+icoord];
 825   }
 826   delete[]local;
 827   delete[] ip;
 828 }
 829 /*! returns constituent entity*/
 830 MED_EN::medEntityMesh MESHCollection::getSubEntity() const
 831 {
 832   return getMeshDimension()==3 ? MED_EN::MED_FACE : MED_EN::MED_EDGE;
 833 }
 834
 835 /*! retrieves the space dimension*/
 836 int MESHCollection::getSpaceDimension() const
 837 {
 838   return _i_non_empty_mesh < 0 ? -1 : _mesh[_i_non_empty_mesh]->getSpaceDimension();
 839 }
 840 /*! retrieves the mesh dimension*/
 841 int MESHCollection::getMeshDimension() const
 842 {
 843   return _i_non_empty_mesh < 0 ? -1 : _mesh[_i_non_empty_mesh]->getMeshDimension();
 844 }
 845
 846 /*! retrieves the type of coordinates system*/
 847 string MESHCollection::getSystem() const
 848 {
 849   return _i_non_empty_mesh < 0 ? "" : _mesh[_i_non_empty_mesh]->getCoordinatesSystem();
 850 }
 851
 852 /*!retrieves the name of the mesh*/
 853 string MESHCollection::getMeshName() const
 854 {
 855   return _i_non_empty_mesh < 0 ? (_mesh[0] ? _mesh[0]->getName() : "") : _mesh[_i_non_empty_mesh]->getName();
 856 }
 857
 858 vector<MEDMEM::MESH*>& MESHCollection::getMesh()
 859 {
 860   return _mesh;
 861 }
 862
 863 MEDMEM::MESH* MESHCollection::getMesh(int idomain) const
 864 {
 865   return _mesh[idomain];
 866 }
 867
 868 vector<MEDMEM::CONNECTZONE*>& MESHCollection::getCZ()
 869 {
 870   return _connect_zones;
 871 }
 872
 873 Topology* MESHCollection::getTopology() const
 874 {
 875   return _topology;
 876 }
 877
 878 void MESHCollection::setTopology(Topology* topo)
 879 {
 880   if (_topology!=0)
 881   {
 882     throw MED_EXCEPTION(STRING("Erreur : topology is already set"));
 883   }
 884   else
 885     _topology = topo;
 886 }
 887
 888 void MESHCollection::setIndivisibleGroup(const string& name)
 889 {
 890   _indivisible_regions.push_back(name);
 891
 892 }
 893
 894 /*! Browses the domains and the regions that have
 895  * been marked as indivisible in order to create a vector
 896  * the dimlension of which is the total number of cells, and
 897  * that contains 0 if the cell belongs to no indivisible group
 898  * and that contains an integer corresponding to the group otherwise.
 899  *
 900  * \param   indivisible_tag on input is an int* allocated as int[nbcells]
 901  *        on output contains the tags
 902  */
 903
 904
 905 void MESHCollection::treatIndivisibleRegions(int* indivisible_tag)
 906 {
 907   //tag 0 is positioned on all the cells that are not affected by these tags
 908   for (int i=0; i<_topology->nbCells(); i++)
 909     indivisible_tag[i]=0;
 910
 911   //treating cell groups
 912   for (int idomain=0; idomain<_topology->nbDomain();idomain++)
 913     for (int igroup=0; igroup<_mesh[idomain]->getNumberOfGroups(MED_EN::MED_CELL); igroup++)
 914       for (unsigned i=0; i<_indivisible_regions.size(); i++)
 915       {
 916         const MEDMEM::GROUP* group = _mesh[idomain]->getGroup(MED_EN::MED_CELL,igroup+1);
 917         string groupname = group->getName();
 918         if (trim(groupname)==trim(_indivisible_regions[i]))
 919         {
 920           int nbcells=group->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
 921           const int* numbers=group->getNumber(MED_EN::MED_ALL_ELEMENTS);
 922           int* global=new int[nbcells];
 923           _topology->convertCellToGlobal(idomain,numbers,nbcells,global);
 924           for (int icell=0; icell<nbcells; icell++)
 925             indivisible_tag[global[icell]-1]=i+1;
 926           delete[] global;
 927         }
 928       }
 929 }
 930
 931 //================================================================================
 932 /*!
 933  * \brief Create cell->node and node->cell connectivities for domains
 934  */
 935 //================================================================================
 936
 937 template<class TID2CONN>
 938 void MESHCollection::fillGlobalConnectivity(TID2CONN & node2cell, TID2CONN & cell2node )
 939 {
 940   for (int idomain=0; idomain<_topology->nbDomain(); idomain++)
 941   {
 942     if ( !_mesh[idomain] ) continue;
 943     //  MED_EN::medGeometryElement* type_array;
 944     int nb_cells= _topology->nbCells(idomain);
 945     int* cell_list = new int[nb_cells];
 946
 947     //retrieving global id list
 948     _topology->getCellList(idomain, cell_list);
 949
 950     int nb_plain_cells = _mesh[idomain]->getNumberOfElements(MED_EN::MED_CELL,
 951                                                              MED_EN::MED_ALL_ELEMENTS);
 952     if (nb_plain_cells >0)
 953     {
 954       const int* conn_index = _mesh[idomain]->getConnectivityIndex(MED_EN::MED_NODAL,
 955                                                                    MED_EN::MED_CELL);
 956
 957       const int* conn =   _mesh[idomain]->getConnectivity(MED_EN::MED_NODAL,
 958                                                           MED_EN::MED_CELL,
 959                                                           MED_EN::MED_ALL_ELEMENTS);
 960       int nbnodes = conn_index[nb_plain_cells]-1;
 961       int* global_nodes =new int [nbnodes];
 962       _topology->convertNodeToGlobal(idomain, conn, nbnodes, global_nodes);
 963       for (int icell=0; icell< nb_plain_cells; icell++)
 964       {
 965         for (int inode=conn_index[icell]; inode < conn_index[icell+1]; inode++)
 966         {
 967           int node_global_id = global_nodes[inode-1];
 968           if ( node_global_id > 0 )
 969           {
 970             int cell_global_id = cell_list[icell];
 971             cell2node [cell_global_id].push_back(node_global_id);
 972             node2cell [node_global_id].push_back(cell_global_id);
 973           }
 974         }
 975       }
 976       delete[] global_nodes;
 977     }
 978
 979     delete[] cell_list;
 980   }
 981 }
 982
 983 /*! Method creating the cell graph
 984  *
 985  * \param array returns the pointer to the structure that contains the graph
 986  * \param edgeweight returns the pointer to the table that contains the edgeweights
 987  *        (only used if indivisible regions are required)
 988  */
 989
 990 void MESHCollection::buildCellGraph(MEDMEM::MEDSKYLINEARRAY* & array,int *& edgeweights )
 991 {
 992
 993   int cell_number=1;
 994   int node_number=1;
 995   for (int i=0; i<_topology->nbDomain(); i++)
 996   {
 997     cell_number+=_topology->getCellNumber(i);
 998     node_number+=_topology->getNodeNumber(i);
 999   }
1000   //list of cells for a given node
1001   //vector< vector<int> > node2cell(node_number);
1002   map< int, vector<int> > node2cell;
1003
1004   //list of nodes for a given cell
1005   //vector< vector <int> > cell2node(cell_number);
1006   map< int, vector <int> > cell2node;
1007
1008   //  map<MED_EN::medGeometryElement,int*> type_cell_list;
1009
1010   //tagging for the indivisible regions
1011   int* indivisible_tag=0;
1012   bool has_indivisible_regions=false;
1013   if (!_indivisible_regions.empty())
1014   {
1015     has_indivisible_regions=true;
1016     indivisible_tag=new int[_topology->nbCells()];
1017     treatIndivisibleRegions(indivisible_tag);
1018   }
1019
1020   fillGlobalConnectivity(node2cell, cell2node );
1021
1022   cout << "beginning of skyline creation"<<endl;
1023   //creating the MEDMEMSKYLINEARRAY containing the graph
1024
1025   int* size = new int[_topology->nbCells()];
1026   int** temp=new int*[_topology->nbCells()];
1027   int** temp_edgeweight=0;
1028   if (has_indivisible_regions)
1029     temp_edgeweight=new int*[_topology->nbCells()];
1030
1031   int cell_glob_shift = 0;
1032
1033   // Get connection to cells on other procs
1034   multimap< int, int > loc2dist; // global cell ids on this proc -> other proc cells
1035   if ( isParallelMode() )
1036   {
1037     cell_glob_shift = _domain_selector->getProcShift();
1038
1039     set<int> loc_domains; // domains on this proc
1040     for ( int idom = 0; idom < _mesh.size(); ++idom )
1041       if ( _mesh[ idom ] )
1042         loc_domains.insert( idom );
1043
1044     for ( int idom = 0; idom < _mesh.size(); ++idom )
1045     {
1046       if ( !_mesh[idom] ) continue;
1047       vector<int> loc2glob_corr; // pairs of corresponding cells (loc_loc & glob_dist)
1048       retrieveDriver()->readLoc2GlobCellConnect(idom, loc_domains, _domain_selector, loc2glob_corr);
1049       //MEDMEM::STRING out;
1050       for ( int i = 0; i < loc2glob_corr.size(); i += 2 )
1051       {
1052         int glob_here  = _topology->convertCellToGlobal(idom,loc2glob_corr[i]);
1053         int glob_there = loc2glob_corr[i+1];
1054         loc2dist.insert ( make_pair( glob_here, glob_there));
1055         //out << glob_here << "-" << glob_there << " ";
1056       }
1057       //cout << "\nRank " << _domain_selector->rank() << ": BndCZ: " << out << endl;
1058     }
1059   }
1060
1061   //going across all cells
1062
1063   map<int,int> cells_neighbours;
1064   for (int i=0; i< _topology->nbCells(); i++)
1065   {
1066
1067
1068     vector<int> cells(50);
1069
1070     for (unsigned inode=0; inode< cell2node[i+1].size(); inode++)
1071     {
1072       int nodeid=cell2node[i+1][inode];
1073       for (unsigned icell=0; icell<node2cell[nodeid].size();icell++)
1074         cells_neighbours[node2cell[nodeid][icell]]++;
1075     }
1076     size[i]=0;
1077     int dimension = getMeshDimension();
1078     cells.clear();
1079
1080     for (map<int,int>::const_iterator iter=cells_neighbours.begin(); iter != cells_neighbours.end(); iter++)
1081     {
1082       if (iter->second >= dimension && iter->first != i+1)
1083       {
1084         cells.push_back(iter->first + cell_glob_shift);
1085         //       cells[isize++]=iter->first;
1086       }
1087     }
1088     // add neighbour cells from distant domains
1089     multimap< int, int >::iterator loc_dist = loc2dist.find( i+1 );
1090     for (; loc_dist!=loc2dist.end() && loc_dist->first==( i+1 ); ++loc_dist )
1091       cells.push_back( loc_dist->second );
1092
1093     size[i]=cells.size();
1094
1095     temp[i]=new int[size[i]];
1096     if (has_indivisible_regions)
1097       temp_edgeweight[i]=new int[size[i]];
1098     //
1099     int itemp=0;
1100
1101     for (vector<int>::const_iterator iter=cells.begin(); iter!=cells.end();iter++)
1102     {
1103       temp[i][itemp]=*iter;
1104       if (has_indivisible_regions)
1105       {
1106         int tag1 = indivisible_tag[(i+1)-1];
1107         int tag2 = indivisible_tag[*iter-1];
1108         if (tag1==tag2 && tag1!=0)
1109           temp_edgeweight[i][itemp]=_topology->nbCells()*100000;
1110         else
1111           temp_edgeweight[i][itemp]=1;
1112       }
1113       itemp++;
1114     }
1115     cells_neighbours.clear();
1116   }
1117   cout <<"end of graph definition"<<endl;
1118   int* index=new int[_topology->nbCells()+1];
1119   index[0]=1;
1120   for (int i=0; i<_topology->nbCells(); i++)
1121     index[i+1]=index[i]+size[i];
1122
1123   node2cell.clear();
1124   cell2node.clear();
1125   if (indivisible_tag!=0) delete [] indivisible_tag;
1126
1127   //SKYLINEARRAY structure holding the cell graph
1128   array= new MEDMEM::MEDSKYLINEARRAY(_topology->nbCells(),index[_topology->nbCells()]-index[0]);
1129   array->setIndex(index);
1130
1131   for (int i=0; i<_topology->nbCells(); i++)
1132   {
1133     array->setI(i+1,temp[i]);
1134     delete[]temp[i];
1135   }
1136
1137   if (has_indivisible_regions)
1138   {
1139     edgeweights=new int[array->getLength()];
1140     for (int i=0; i<_topology->nbCells(); i++)
1141     {
1142       for (int j=index[i]; j<index[i+1];j++)
1143         edgeweights[j-1]=temp_edgeweight[i][j-index[i]];
1144       delete[] temp_edgeweight[i];
1145     }
1146     delete[]temp_edgeweight;
1147   }
1148   delete[] index;
1149   delete[] temp;
1150   delete[] size;
1151
1152   cout<< "end of graph creation"<<endl;
1153 }
1154
1155 /*! Creates the partition corresponding to the cell graph and the partition number
1156  *
1157  * \param nbdomain number of subdomains for the newly created graph
1158  *
1159  * returns a topology based on the new graph
1160  */
1161 Topology* MESHCollection::createPartition(int nbdomain,
1162                                           Graph::splitter_type split,
1163                                           const string& options_string,
1164                                           int* user_edge_weights,
1165                                           int* user_vertices_weights)
1166 {
1167   if (nbdomain <1) throw MEDEXCEPTION("Number of subdomains must be >0");
1168   MEDMEM::MEDSKYLINEARRAY* array=0;
1169   int* edgeweights=0;
1170
1171   MESSAGE_MED("Building cell graph");
1172   buildCellGraph(array,edgeweights);
1173
1174   switch (split)
1175   {
1176   case Graph::METIS:
1177 #if defined(MED_ENABLE_METIS) || defined(MED_ENABLE_PARMETIS)
1178     _cell_graph=boost::shared_ptr<Graph>(new METISGraph(array,edgeweights));
1179 #else
1180     throw MEDEXCEPTION("METIS Graph is not available. Check your products, please.");
1181 #endif
1182     break;
1183   case Graph::SCOTCH:
1184 #ifdef MED_ENABLE_SCOTCH
1185     _cell_graph=boost::shared_ptr<Graph>(new SCOTCHGraph(array,edgeweights));
1186 #else
1187     throw MEDEXCEPTION("SCOTCH Graph is not available. Check your products, please.");
1188 #endif
1189     break;
1190   }
1191
1192   //!user-defined weights
1193   if (user_edge_weights!=0)
1194     _cell_graph->setEdgesWeights(user_edge_weights);
1195   if (user_vertices_weights!=0)
1196     _cell_graph->setVerticesWeights(user_vertices_weights);
1197
1198   MESSAGE_MED("Partitioning graph");
1199   _cell_graph->partGraph(nbdomain,options_string,_domain_selector);
1200
1201   // DEBUG
1202 //   MEDMEM::STRING out("RESULT GRAPH #");
1203 //   out << (_domain_selector?_domain_selector->rank():0) << ": ";
1204 //   const int* part = _cell_graph->getPart();
1205 //   int n = _cell_graph->nbVertices();
1206 //   for ( int e=0; e < n; ++e )
1207 //     out << part[e] <<" ";
1208 //   cout << out << endl;
1209
1210
1211   MESSAGE_MED("Building new topology");
1212   //_cell_graph is a shared pointer
1213   Topology* topology = new ParallelTopology (_cell_graph, nbdomain, getMeshDimension());
1214
1215   //cleaning
1216   if (edgeweights!=0) delete[] edgeweights;
1217   //if (array!=0) delete array;
1218   MESSAGE_MED("End of partition creation");
1219   return topology;
1220 }
1221
1222 /*! Creates a topology for a partition specified by the user
1223  *
1224  * \param table user-specified partition (for each cell contains the domain number from 0 to n-1)
1225  *
1226  * returns a topology based on the new partition
1227  */
1228 Topology* MESHCollection::createPartition(const int* partition)
1229 {
1230   MEDMEM::MEDSKYLINEARRAY* array=0;
1231   int* edgeweights=0;
1232
1233   buildCellGraph(array,edgeweights);
1234
1235   set<int> domains;
1236   for (int i=0; i<_topology->nbCells(); i++)
1237   {
1238     domains.insert(partition[i]);
1239   }
1240   int nbdomain=domains.size();
1241
1242   _cell_graph=boost::shared_ptr<Graph>(new UserGraph(array, partition, _topology->nbCells()));
1243
1244   //_cell_graph is a shared pointer
1245   Topology* topology = new ParallelTopology (_cell_graph, nbdomain, getMeshDimension());
1246
1247   //if (array!=0) delete array;
1248   return topology;
1249 }
1250
1251
1252 /*! building Connect Zones for storing the informations
1253  * of the connectivity
1254  *
1255  * The connect zones are created for every domain that has common nodes with
1256  * domain \a idomain
1257  *
1258  * \param idomain domain number for which the connect zones are created
1259  * */
1260
1261 // void MESHCollection::buildConnectZones(int idomain)
1262 // {
1263 //   // constructing node/node correspondencies
1264 //   vector<MEDMEM::MEDSKYLINEARRAY*> node_node_correspondency;
1265 //   node_node_correspondency.resize(_topology->nbDomain());
1266
1267 //   cout << "Computing node/node corresp"<<endl;
1268
1269 //   _topology->computeNodeNodeCorrespondencies(idomain, node_node_correspondency );
1270
1271 //   for (int idistant=0; idistant< _topology->nbDomain(); idistant++)
1272 //   {
1273 //     // on regarde si une correspondance noeud/noeud a été trouvée
1274 //     // entre idomain et idistant
1275 //     // si oui, on crée une connectzone
1276 //     if (node_node_correspondency[idistant]!=0)
1277 //     {
1278 //       MEDMEM::CONNECTZONE* cz= new MEDMEM::CONNECTZONE();
1279 //       cz->setLocalMesh(_mesh[idomain]);
1280 //       cz->setDistantMesh(_mesh[idistant]);
1281 //       cz->setLocalDomainNumber(idomain);
1282 //       cz->setDistantDomainNumber(idistant);
1283 //       cz-> setName ("Connect zone defined by SPLITTER");
1284 //       cz->setNodeCorresp(node_node_correspondency[idistant]);
1285 //       _connect_zones.push_back(cz);
1286 //     }
1287 //   }
1288 //   cout << "Computing node/node corresp"<<endl;
1289
1290 //   vector<MEDMEM::MEDSKYLINEARRAY*> cell_cell_correspondency;
1291 //   cell_cell_correspondency.resize(_topology->nbDomain());
1292 //   _topology->computeCellCellCorrespondencies(idomain, cell_cell_correspondency, _cell_graph.get());
1293
1294 //   for (int idistant=0; idistant< _topology->nbDomain(); idistant++)
1295 //   {
1296 //     //the connect zone has been created by the node/node computation
1297 //     if (cell_cell_correspondency[idistant]!=0)
1298 //     {
1299 //       MEDMEM::CONNECTZONE* cz=0;
1300 //       for (int icz=0; icz<_connect_zones.size();icz++)
1301 //         if (_connect_zones[icz]->getLocalDomainNumber()==idomain &&
1302 //             _connect_zones[icz]->getDistantDomainNumber()==idistant)
1303 //           cz = _connect_zones[icz];
1304 //       if (cz!=0)
1305 //         cz->setEntityCorresp(MED_EN::MED_CELL,MED_EN::MED_CELL, cell_cell_correspondency[idistant]);
1306 //       else
1307 //         throw MEDEXCEPTION("MESHCollection::buildConnectZones() -A connect zone should exist");
1308 //       //delete cell_cell_correspondency[idistant];
1309 //     }
1310
1311 //   }
1312 // }
1313
1314 //================================================================================
1315 /*!
1316  * \brief Adds a group of joint faces
1317  *  \param loc_face_ids - local numbers of faces
1318  *  \param idomian - domain index where faces are local
1319  *  \param idistant - the other domain index
1320  */
1321 //================================================================================
1322
1323 void MESHCollection::addJointGroup(const std::vector<int>& loc_face_ids, int idomain, int idistant)
1324 {
1325   MEDMEM::MESHING* meshing = dynamic_cast<MEDMEM::MESHING*> (_mesh[idomain]);
1326   MED_EN::medEntityMesh constituent_entity = getSubEntity();
1327
1328   MEDMEM::STRING jointname("joint_");
1329   jointname<<idistant+1;
1330
1331   MEDMEM::GROUP * tmp_grp = new GROUP, * joint_group = tmp_grp;
1332   // try to find already present group with such a name
1333   //  vector<MEDMEM::GROUP*> groups = meshing->getGroups( constituent_entity );
1334   //  for ( int g = 0; g < groups.size(); ++g )
1335   //    if ( groups[g]->getName() == jointname.str() )
1336   //    {
1337   //      joint_group = groups[g];
1338   //      break;
1339   //    }
1340   // assure uniqueness of group name
1341   bool unique = false;
1342   vector<MEDMEM::GROUP*> groups = meshing->getGroups( constituent_entity );
1343   do
1344   {
1345     unique = true;
1346     for ( int g = 0; unique && g < groups.size(); ++g )
1347       unique = ( groups[g]->getName() != jointname );
1348     if ( !unique )
1349       jointname << "_" << idomain+1;
1350   }
1351   while ( !unique );
1352   joint_group->setName(jointname);
1353   joint_group->setMesh(meshing);
1354   joint_group->setEntity(constituent_entity);
1355   map<MED_EN::medGeometryElement, vector<int> > joint_types;
1356
1357   int nbfaces = loc_face_ids.size();
1358   for (int i=0; i<nbfaces; i++)
1359   {
1360     MED_EN::medGeometryElement type = meshing->getElementType(constituent_entity,loc_face_ids[i]);
1361     joint_types[type].push_back(loc_face_ids[i]);
1362   }
1363   joint_group->setNumberOfGeometricType(joint_types.size());
1364   MED_EN::medGeometryElement* types=new MED_EN::medGeometryElement[joint_types.size()];
1365   int* nb_in_types=new int[joint_types.size()];
1366   int* group_index=new int[joint_types.size()+1];
1367
1368   group_index[0]=1;
1369   int itype=0;
1370   int iface =0;
1371   int* group_value=new int[nbfaces];
1372   for (map<MED_EN::medGeometryElement, vector<int> >::const_iterator iterj=joint_types.begin();
1373        iterj != joint_types.end();
1374        iterj++)
1375   {
1376     nb_in_types[itype]=(iterj->second).size();
1377     types[itype]=iterj->first;
1378     itype++;
1379     group_index[itype]=group_index[itype-1]+(iterj->second).size();
1380     for (int i=0; i<  (iterj->second).size(); i++)
1381       group_value[iface++]=(iterj->second)[i];
1382   }
1383   joint_group->setGeometricType(types);
1384   joint_group->setNumberOfElements(nb_in_types);
1385   joint_group->setNumber(group_index, group_value, /*shallowCopy=*/true);
1386   delete[] types;
1387   delete[] nb_in_types;
1388
1389   if ( joint_group == tmp_grp )
1390     meshing->addGroup(*tmp_grp);
1391   tmp_grp->removeReference();
1392 }
1393
1394 /*! building Connect Zones for storing the informations
1395  * of the connectivity
1396  * */
1397
1398 void MESHCollection::buildConnectZones()
1399 {
1400   vector <map <MED_EN::medGeometryElement, vector<MEDSPLITTER_FaceModel*> > > face_map(_topology->nbDomain());
1401   map< pair<int,int>, MEDMEM::MEDSKYLINEARRAY*> cell_corresp_here;
1402
1403   MED_EN::medEntityMesh constituent_entity = getSubEntity();
1404
1405   if ( isParallelMode() )
1406   {
1407     buildConnectZonesBetweenProcs(face_map, cell_corresp_here);
1408   }
1409
1410   cout << "Computing node/node corresp"<<endl;
1411
1412   //Creating nodes
1413   for (int idomain=0; idomain<_topology->nbDomain(); idomain++)
1414   {
1415
1416     // constructing node/node correspondencies
1417     vector<MEDMEM::MEDSKYLINEARRAY*> node_node_correspondency(_topology->nbDomain());
1418     _topology->computeNodeNodeCorrespondencies(idomain, node_node_correspondency );
1419
1420     for (int idistant=0; idistant< _topology->nbDomain(); idistant++)
1421     {
1422       // on regarde si une correspondance noeud/noeud a été trouvée
1423       // entre idomain et idistant
1424       // si oui, on crée une connectzone
1425       if (node_node_correspondency[idistant]!=0)
1426       {
1427         MEDMEM::CONNECTZONE* cz= new MEDMEM::CONNECTZONE();
1428         cz->setLocalMesh(_mesh[idomain]);
1429         cz->setDistantMesh(_mesh[idistant]);
1430         cz->setLocalDomainNumber(idomain);
1431         cz->setDistantDomainNumber(idistant);
1432         cz-> setName ("Connect zone defined by SPLITTER");
1433         cz->setNodeCorresp(node_node_correspondency[idistant]);
1434         _connect_zones.push_back(cz);
1435       }
1436     }
1437   }
1438   cout << "Computing face corresp"<<endl;
1439
1440   //creating faces if required
1441   if (_subdomain_boundary_creates)
1442   {
1443     int global_face_id = _topology->getFaceNumber()+1;
1444     //int global_face_id = _topology->getMaxGlobalFace()+1;
1445
1446     map <pair<int,int>, vector<int> > faces_in_joint;
1447
1448     if ( !isParallelMode() )
1449       // taking faces that are already present in the mesh into account
1450       for (int idomain=0; idomain<_topology->nbDomain(); idomain++)
1451       {
1452         getFaces(idomain,face_map[idomain]);
1453       }
1454
1455     // creating faces that are located at the interface between
1456     // subdomains
1457
1458     vector <int> nb_added_groups( _topology->nbDomain(), 0 );
1459
1460     for (int idomain=0; idomain<_topology->nbDomain(); idomain++)
1461     {
1462       vector<MEDMEM::MEDSKYLINEARRAY*> cell_cell_correspondency( _topology->nbDomain() );
1463       if ( !isParallelMode() )
1464         _topology->computeCellCellCorrespondencies(idomain, cell_cell_correspondency, _cell_graph.get());
1465
1466       for (int idistant=0; idistant< _topology->nbDomain(); idistant++)
1467       {
1468         if (idistant <= idomain) continue;
1469
1470         MEDMEM::MEDSKYLINEARRAY* cell_correspondency = 0;
1471         if ( isParallelMode() )
1472           cell_correspondency = cell_corresp_here[ make_pair (idomain,idistant)];
1473         else
1474           cell_correspondency = cell_cell_correspondency[idistant];
1475
1476         //the connect zone has been created by the node/node computation
1477
1478         if ( cell_correspondency )
1479         {
1480           int nbcells      = cell_correspondency->getNumberOf();
1481           const int* index = cell_correspondency->getIndex();
1482           const int* value = cell_correspondency->getValue();
1483           if ( isParallelMode() )
1484             global_face_id = _domain_selector->getFisrtGlobalIdOfSubentity( idomain, idistant );
1485
1486           for (int ilocal=0; ilocal<nbcells; ilocal++)
1487           {
1488             for (int icelldistant = index[ilocal]; icelldistant < index[ilocal+1]; icelldistant++)
1489             {
1490               int distant_id = value[icelldistant-1];
1491               MEDSPLITTER_FaceModel* face = getCommonFace(idomain,ilocal+1,idistant,distant_id,global_face_id);
1492               face_map[idomain][face->getType()].push_back(face);
1493               MEDSPLITTER_FaceModel* face2 = getCommonFace(idistant,distant_id,idomain, ilocal+1,global_face_id);
1494               face_map[idistant][face->getType()].push_back(face2);
1495               faces_in_joint[make_pair(idomain,idistant)].push_back(global_face_id);
1496               global_face_id++;
1497             }
1498           }
1499         }
1500
1501       }
1502       //cleaning up
1503       for (int idistant=0; idistant< _topology->nbDomain(); idistant++)
1504         delete cell_cell_correspondency[idistant];
1505     }
1506
1507
1508     _topology->recreateFaceMapping(face_map);
1509
1510     //transforming the face_map into a constituent entity connectivity
1511     for (int idomain=0; idomain< _topology->nbDomain();idomain++)
1512     {
1513       int nbtypes = face_map[idomain].size();
1514       vector<medGeometryElement> types;
1515       vector <int> nb_elems;
1516       vector <int*> conn;
1517
1518       MEDMEM::MESHING* meshing = dynamic_cast<MEDMEM::MESHING*> (_mesh[idomain]);
1519       if ( !meshing->getConnectivityptr() )
1520         continue; // no cells in idomain
1521
1522       for (map <medGeometryElement, vector<MEDSPLITTER_FaceModel*> >::const_iterator iter= face_map[idomain].begin();
1523            iter != face_map[idomain].end(); iter ++)
1524       {
1525         types.push_back(iter->first);
1526         int nb_elem_in_type = (iter->second).size();
1527         nb_elems.push_back(nb_elem_in_type);
1528         int nb_node_per_type=(iter->first)%100;
1529         int* connectivity= new int [nb_node_per_type*nb_elem_in_type];
1530         for (int ielem=0; ielem<nb_elem_in_type; ielem++)
1531         {
1532           for (int inode=0;  inode<nb_node_per_type; inode++)
1533             connectivity[ielem*nb_node_per_type+inode]=(*(iter->second)[ielem])[inode];
1534         }
1535         conn.push_back(connectivity);
1536
1537       }
1538       //setting the faces in the mesh
1539       meshing->setNumberOfTypes(nbtypes,constituent_entity);
1540       meshing->setTypes(&types[0],constituent_entity);
1541       meshing->setNumberOfElements(&nb_elems[0],constituent_entity);
1542
1543       for (int itype=0; itype<nbtypes; itype++)
1544       {
1545         meshing->setConnectivity( constituent_entity, types[itype], conn[itype] );
1546         delete[]conn[itype];
1547       }
1548       for (int idistant =0; idistant<_topology->nbDomain(); idistant++)
1549       {
1550         map <pair<int,int>, vector<int> >::iterator iter;
1551         iter = faces_in_joint.find(make_pair(idomain,idistant));
1552         if (iter == faces_in_joint.end())
1553         {
1554           iter = faces_in_joint.find (make_pair(idistant,idomain));
1555           if (iter == faces_in_joint.end())
1556             continue;
1557         }
1558
1559         int nbfaces = (iter->second).size();
1560         vector<int> face_joint(nbfaces*2);
1561         MEDMEM::CONNECTZONE* cz=0;
1562         for (unsigned icz=0; icz<_connect_zones.size();icz++)
1563           if (_connect_zones[icz]->getLocalDomainNumber()==idomain &&
1564               _connect_zones[icz]->getDistantDomainNumber()==idistant)
1565             cz = _connect_zones[icz];
1566
1567         int nbtotalfaces= _topology->getFaceNumber(idomain);
1568
1569         //creating arrays for the MEDSKYLINEARRAY structure containing the joint
1570         int* index =new int[nbtotalfaces+1];
1571         for (int i=0; i<nbtotalfaces+1;i++)
1572           index[i]=0;
1573         int*value=new int[nbfaces];
1574
1575         map<int,int> faces;
1576         vector<int> local_faces( nbfaces );
1577         for (int iface=0; iface<nbfaces; iface++)
1578         {
1579           int iglobal = (iter->second)[iface];
1580           int localid=_topology->convertGlobalFace(iglobal,idomain);
1581           int distantid=_topology->convertGlobalFace(iglobal,idistant);
1582           faces.insert(make_pair(localid,distantid));
1583           local_faces[iface]=localid;
1584         }
1585
1586         int iloc=0;
1587         index[0]=1;
1588         for (map<int,int>::const_iterator iter=faces.begin();
1589              iter != faces.end();
1590              iter++)
1591         {
1592           index[iter->first]=1;
1593           value[iloc++]=iter->second;
1594         }
1595
1596         for (int i=0; i<nbtotalfaces;i++)
1597           index[i+1]+=index[i];
1598         bool shallowcopy=true;
1599         MEDMEM::MEDSKYLINEARRAY* skarray=new MEDMEM::MEDSKYLINEARRAY(nbtotalfaces,nbfaces,index,value,shallowcopy);
1600
1601         if (cz!=0)
1602           cz->setEntityCorresp(constituent_entity,constituent_entity,skarray);
1603         else
1604           throw MEDEXCEPTION("MESHCollection::buildConnectZones() -A connect zone should exist");
1605         // Creating a group of the faces constituting the joint
1606         addJointGroup( local_faces, idomain, idistant );
1607         nb_added_groups[ idomain ]++;
1608       }
1609     }
1610
1611     if ( isParallelMode() )
1612     {
1613       // Now all faces have got local ids and we can receive local ids from other procs.
1614       // Set face/face data to zones with other procs and create a group
1615       for (int icz=0; icz<_connect_zones.size();icz++)
1616       {
1617         MEDMEM::CONNECTZONE* cz=_connect_zones[icz];
1618         if ( _domain_selector->isMyDomain( cz->getDistantDomainNumber()) ) continue;
1619
1620         int glob_id = _domain_selector->getFisrtGlobalIdOfSubentity( cz->getLocalDomainNumber(),
1621                                                                      cz->getDistantDomainNumber());
1622         int nb_cz_faces = _domain_selector->getNbCellPairs( cz->getDistantDomainNumber(),
1623                                                             cz->getLocalDomainNumber());
1624         vector< int > loc_ids_here( nb_cz_faces );
1625         for ( int i = 0; i < nb_cz_faces; ++i )
1626           loc_ids_here[i] = _topology->convertGlobalFace(glob_id++,cz->getLocalDomainNumber());
1627
1628         int* loc_ids_dist = _domain_selector->exchangeSubentityIds( cz->getLocalDomainNumber(),
1629                                                                     cz->getDistantDomainNumber(),
1630                                                                     loc_ids_here );
1631         int nb_faces_here= _topology->getFaceNumber(cz->getLocalDomainNumber());
1632         int* face_index = new int[ nb_faces_here+1 ];
1633         face_index[0]=1;
1634         for ( int loc_id = 0, i = 0; loc_id < nb_faces_here; ++loc_id)
1635         {
1636           face_index[ loc_id+1 ] = face_index[ loc_id ];
1637           if ( i < loc_ids_here.size() && loc_ids_here[i] == loc_id+1 )
1638           {
1639             face_index[ loc_id+1 ]++;
1640             i++;
1641           }
1642         }
1643         MEDMEM::MEDSKYLINEARRAY* skarray=
1644           new MEDMEM::MEDSKYLINEARRAY(nb_faces_here, nb_cz_faces, face_index, loc_ids_dist, true);
1645         cz->setEntityCorresp(constituent_entity,constituent_entity,skarray);
1646
1647         addJointGroup( loc_ids_here, cz->getLocalDomainNumber(), cz->getDistantDomainNumber());
1648         nb_added_groups[ cz->getLocalDomainNumber() ]++;
1649       }
1650     }
1651
1652     for (int idomain=0; idomain<_topology->nbDomain(); idomain++)
1653     {
1654       // delete face_map
1655       for (map <medGeometryElement, vector<MEDSPLITTER_FaceModel*> >::const_iterator iter= face_map[idomain].begin();
1656            iter != face_map[idomain].end(); iter ++)
1657         for (unsigned i=0; i<(iter->second).size();i++)
1658           delete (iter->second)[i];
1659
1660       if ( nb_added_groups[ idomain ] > 0 &&
1661            _mesh[idomain]->getNumberOfFamilies( constituent_entity ) > 0 )
1662         // needed because if there were face families before, driver won't
1663         // create families from just added groups (see MEDMEM_MedMeshDriver.cxx:3330),
1664         // actually it is a bug of driver - it must check presence of groups in families
1665         _mesh[idomain]->createFamilies();
1666     }
1667   }
1668
1669   if ( isParallelMode() )
1670     // Excange info on types of constituent_entity needed while writing joints
1671     // to get ids local in geom type for distant procs
1672     _domain_selector->gatherEntityTypesInfo( _mesh, constituent_entity );
1673
1674   cout << "Computing cell/cell corresp"<<endl;
1675
1676   //Creating cell/cell correspondencies
1677   for (int idomain=0;idomain<_topology->nbDomain();idomain++)
1678   {
1679     vector<MEDMEM::MEDSKYLINEARRAY*> cell_cell_correspondency( _topology->nbDomain() );
1680     if ( !isParallelMode() )
1681       _topology->computeCellCellCorrespondencies(idomain,cell_cell_correspondency,_cell_graph.get());
1682
1683     for (int idistant=0; idistant< _topology->nbDomain(); idistant++)
1684     {
1685       MEDMEM::MEDSKYLINEARRAY* cell_correspondency = 0;
1686       if ( isParallelMode() )
1687         cell_correspondency = cell_corresp_here[ make_pair (idomain,idistant)];
1688       else
1689         cell_correspondency = cell_cell_correspondency[idistant];
1690
1691       //the connect zone has been created by the node/node computation
1692       if ( cell_correspondency )
1693       {
1694         MEDMEM::CONNECTZONE* cz=0;
1695         for (unsigned icz=0; icz<_connect_zones.size();icz++)
1696           if (_connect_zones[icz]->getLocalDomainNumber()==idomain &&
1697               _connect_zones[icz]->getDistantDomainNumber()==idistant)
1698             cz = _connect_zones[icz];
1699         if (cz!=0)
1700           cz->setEntityCorresp(MED_EN::MED_CELL,MED_EN::MED_CELL, cell_correspondency);
1701         else
1702           throw MEDEXCEPTION("MESHCollection::buildConnectZones() -A connect zone should exist");
1703       }
1704     }
1705   }
1706 }
1707
1708 /*! building Connect Zones for storing the informations
1709  * of the connectivity in the parallel mode
1710  * */
1711
1712 void MESHCollection::buildConnectZonesBetweenProcs(TGeom2FacesByDomian & face_map,
1713                                                    map< pair<int,int>, MEDMEM::MEDSKYLINEARRAY*> & cell_cell_correspondency_here)
1714 {
1715   using namespace MED_EN;
1716
1717   // graph over all procs
1718   auto_ptr<Graph> global_graph( _domain_selector->gatherGraph( _cell_graph.get() ));
1719
1720   vector< vector< JointExchangeData > > joints_of_domain( _topology->nbDomain() );
1721
1722   for (int idomain=0; idomain<_topology->nbDomain(); idomain++)
1723   {
1724     if ( !_domain_selector->isMyDomain( idomain )) continue;
1725
1726     vector< JointExchangeData > & joints = joints_of_domain[ idomain ];
1727     joints.resize( _topology->nbDomain() );
1728
1729     // Find corresponding cells on other procs
1730
1731     const int* gra_index = global_graph->getGraph()->getIndex();
1732     const int* gra_value = global_graph->getGraph()->getValue();
1733     const int* partition = global_graph->getPart();
1734     const int dj = gra_index[0];
1735
1736     vector< int > glob_cells_here( _topology->getCellNumber( idomain ));
1737     _topology->getCellList( idomain, & glob_cells_here[0]);
1738     for ( int loc_here = 0; loc_here < glob_cells_here.size(); ++loc_here )
1739     {
1740       int glob_here = glob_cells_here[ loc_here ];
1741       for ( int j = gra_index[ glob_here-1 ]; j < gra_index[ glob_here ]; ++j )
1742       {
1743         int glob_neighbor = gra_value[ j-dj ];
1744         int neighbor_dom = partition[ glob_neighbor-1 ];
1745         if ( neighbor_dom == idomain ) continue;
1746
1747         if ( _domain_selector->isMyDomain( neighbor_dom ))
1748         {
1749           joints[ neighbor_dom ].addCellCorrespondence
1750             (_mesh[idomain], neighbor_dom, idomain, glob_neighbor, glob_here, loc_here + 1,
1751              _topology->convertGlobalCell(glob_neighbor).second );
1752         }
1753         else
1754         {
1755           joints[ neighbor_dom ].addCellCorrespondence
1756             (_mesh[idomain], neighbor_dom, idomain, glob_neighbor, glob_here, loc_here + 1 );
1757         }
1758       }
1759     }
1760   }
1761   global_graph.reset(); // free memory
1762
1763   // set joints in a queue to exchange
1764   typedef map< int, JointExchangeData* > TOrderedJoints;
1765   TOrderedJoints queue;
1766   for (int idomain=0; idomain<_topology->nbDomain(); idomain++)
1767   {
1768     if ( !_domain_selector->isMyDomain( idomain )) continue;
1769
1770     vector< JointExchangeData > & joints = joints_of_domain[ idomain ];
1771     for (int idist=0; idist<_topology->nbDomain(); ++idist )
1772     {
1773       JointExchangeData& joint = joints[idist];
1774
1775       int nb_cell_pairs = joint.nbCellPairs();
1776       if ( nb_cell_pairs == 0 )
1777         continue;
1778       else
1779         _domain_selector->setNbCellPairs( nb_cell_pairs, idist, idomain );
1780
1781       joint.setMeshes( idist, _mesh[idist], idomain, _mesh[idomain] );
1782
1783       if ( _domain_selector->isMyDomain( idist ))
1784       {
1785         // a joint on this proc
1786         cell_cell_correspondency_here[ make_pair( idomain, idist )] = joint.makeCellCorrespArray();
1787       }
1788       else
1789       {
1790         // a joint with distant proc
1791         joint.setConnectivity( & ((MEDMEM::MeshFuse*)_mesh[idomain])->getNodeNumbers()[0] );
1792         int order = _domain_selector->jointId( idomain, idist );
1793         queue[ order ] = & joint;
1794       }
1795     }
1796   }
1797   // gather info on cell geom types needed to exchange joints
1798   _domain_selector->gatherEntityTypesInfo( _mesh, MED_EN::MED_CELL );
1799
1800   // gather info on nb of sub-entities to compute their global numbers for joints
1801   _domain_selector->gatherNbOf( getSubEntity(), _mesh );
1802   _domain_selector->gatherNbCellPairs();
1803   if ( _subdomain_boundary_creates )
1804   {
1805     // taking faces that are already present in the mesh into account
1806     for (int idomain=0; idomain<_topology->nbDomain(); idomain++)
1807       if ( _domain_selector->isMyDomain( idomain ))
1808         getFaces(idomain,face_map[idomain]);
1809   }
1810   else
1811   {
1812     face_map.clear(); // mark for the joint not to create faces
1813   }
1814
1815   // exchange joint data with other procs and make CONNECTZONEs
1816   TOrderedJoints::iterator ord_joint = queue.begin();
1817   for ( ; ord_joint != queue.end(); ++ord_joint )
1818   {
1819     JointExchangeData* joint = ord_joint->second;
1820
1821     _domain_selector->exchangeJoint( joint );
1822     if ( _subdomain_boundary_creates )
1823     {
1824       int first_sub_id = _domain_selector->getFisrtGlobalIdOfSubentity( joint->localDomain(),
1825                                                                         joint->distantDomain() );
1826       joint->setFisrtGlobalIdOfSubentity( first_sub_id );
1827     }
1828     _connect_zones.push_back ( joint->makeConnectZone( face_map ));
1829   }
1830 }
1831
1832 /*! projects old collection families on new collection families
1833  */
1834 void MESHCollection::castFamilies(const MESHCollection& old_collection)
1835 {
1836   vector <list<int> > element_array  (_topology->nbDomain());
1837
1838   //loop on old domains to create groups out of the existing families
1839   if (_family_splitting)
1840     for (int idomain=0; idomain < old_collection._topology->nbDomain(); idomain++)
1841       old_collection.getMesh(idomain)->createGroups();
1842
1843   //definition of the entities array which
1844   //defines the entities over which the information is cast
1845   MED_EN::medEntityMesh entities[3];
1846   entities[0]=MED_EN::MED_NODE;
1847   entities[1]=getSubEntity();
1848   entities[2]=MED_EN::MED_CELL;
1849
1850   for (int ientity=0; ientity<=2;ientity++)
1851   {
1852
1853     //int nbgroups = old_collection.getMesh(0)->getNumberOfGroups(entities[ientity]);
1854
1855     map <string, set<int> > group_map;
1856     for (int idomain=0; idomain < old_collection._topology->nbDomain(); idomain++)
1857     {
1858       if ( !old_collection.getMesh(idomain) ) continue;
1859       for (int igroup=0; igroup<old_collection.getMesh(idomain)->getNumberOfGroups(entities[ientity]); igroup++)
1860       {
1861         //retrieves a group
1862         MEDMEM::GROUP* group = (old_collection.getMesh(idomain)->getGroups(entities[ientity]))[igroup];
1863         //increments the number of groups if it is a new group
1864         //if (group_map.find(group->getName())==group_map.end())
1865
1866         group_map[group->getName()].insert(idomain);
1867         //   group_map.insert(make_pair(group->getName(), idomain);
1868
1869       }
1870     }
1871     int nbgroups=group_map.size();
1872     vector <int> igroupold(old_collection._topology->nbDomain(),0);
1873     map<string,set<int> >::const_iterator iter=group_map.begin();
1874
1875     for (int igroup=0; igroup<nbgroups; igroup++)
1876     {
1877       vector <const MEDMEM::SUPPORT*> old_supports(old_collection._topology->nbDomain());
1878       string group_name = iter->first;
1879       iter++;
1880
1881       //parameters stored for passing group description
1882       // from the old meshes to the new ones
1883
1884       for (int idomain=0; idomain < old_collection._topology->nbDomain(); idomain++)
1885       {
1886         //                for (set<int>::iterator iter=group_map[group_name].begin(); iter!=group_map[group_name].end(); iter++)
1887         //                cout << *iter<<" ";
1888         //                cout <<endl;
1889         if (group_map[group_name].find(idomain)==group_map[group_name].end()) continue;
1890
1891         //retrieves the group igroup on domain idomain
1892         MEDMEM::GROUP* group = (old_collection.getMesh(idomain)->getGroups(entities[ientity]))[igroupold[idomain]];
1893         old_supports[idomain] = static_cast<const MEDMEM::SUPPORT*> (group);
1894         igroupold[idomain]++;
1895       }
1896
1897       vector <MEDMEM::GROUP*>new_groups(_topology->nbDomain());
1898       vector <MEDMEM::SUPPORT*> new_supports(_topology->nbDomain());
1899       for (int i=0; i<_topology->nbDomain(); i++)
1900       {
1901         new_groups[i]=new MEDMEM::GROUP();
1902         new_supports[i]=static_cast<MEDMEM::SUPPORT*>(new_groups[i]);
1903       }
1904       castSupport(old_collection,old_supports,new_supports);
1905
1906       //creating new groups from the previous list of elements
1907       for (int idomain=0; idomain <_topology->nbDomain(); idomain++)
1908       {
1909         MEDMEM::MESHING* mesh_builder=static_cast<MEDMEM::MESHING*> (_mesh[idomain]);
1910         if ( new_supports[idomain] )
1911           mesh_builder->addGroup(*new_groups[idomain]);
1912       }
1913       //groups are copied by the addGroup method,
1914       //so they can be safely deleted here
1915       for (int i=0; i<_topology->nbDomain(); i++)
1916       {
1917         if ( new_supports[i] ) new_groups[i]->removeReference();
1918       }
1919
1920     }// on groups
1921   }//on entities
1922 }
1923
1924
1925 void MESHCollection::castSupport(const MESHCollection& old_collection, vector<const MEDMEM::SUPPORT*>& old_support, vector<MEDMEM::SUPPORT*>& new_support)
1926 {
1927
1928   if (old_collection._topology->nbDomain() != (int)old_support.size())
1929   {
1930     throw MED_EXCEPTION(STRING("Error : wrong call to MESHCollection::castSupport"));
1931   }
1932   vector <list<int> > element_array  (_topology->nbDomain());
1933
1934   //parameters stored for passing description
1935   // from the old meshes to the new ones
1936   string name;
1937   string description;
1938   MED_EN::medEntityMesh entity;
1939   vector <string> support_name(1);
1940   support_name[0]="support";
1941   for (int inew=0; inew< _topology->nbDomain(); inew++)
1942     element_array[inew].clear();
1943
1944   for (int idomain=0; idomain < old_collection._topology->nbDomain(); idomain++)
1945   {
1946     //retrieves the group igroup on domain idomain
1947     const MEDMEM::SUPPORT* support = old_support[idomain];
1948     if (old_support[idomain]==0) continue;
1949     name = support->getName();
1950     description=support->getDescription();
1951     int nbelem = support->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
1952     if (nbelem==0 && !_create_empty_groups) continue;
1953
1954     int* list_of_elems;
1955     if (support->isOnAllElements())
1956     {
1957       list_of_elems = new int[nbelem];
1958       for (int i=0; i<nbelem;i++)
1959         list_of_elems[i]=i+1;
1960     }
1961     else
1962       list_of_elems = const_cast<int*> (support->getNumber(MED_EN::MED_ALL_ELEMENTS));
1963
1964     int* array=new int[nbelem];
1965     int* ip=0;
1966     int* local=0;
1967     int* full_array=0;
1968     entity = support->getEntity();
1969     int size;
1970
1971     switch (entity)
1972     {
1973     case MED_EN::MED_CELL :
1974       ip=new int[nbelem];
1975       local= new int[nbelem];
1976       size=nbelem;
1977       old_collection.getTopology()->convertCellToGlobal(idomain,list_of_elems,nbelem,array);
1978       _topology->convertGlobalCellList(array,nbelem,local,ip);
1979       for (int i=0; i<nbelem; i++)
1980         //              cell_arrays[ip[i]][local[i]]=id;
1981       {
1982         //          cout <<"(glob,ip,iloc)/nbelem"<<array[i]<<" "<<ip[i]<<" "<<local[i]<<"/"<<nbelem<<endl;
1983         element_array[ip[i]].push_back(local[i]);
1984       }
1985       break;
1986     case MED_EN::MED_FACE :
1987     case MED_EN::MED_EDGE :
1988       old_collection.getTopology()->convertFaceToGlobal(idomain,list_of_elems,nbelem,array);
1989       _topology->convertGlobalFaceListWithTwins(array,nbelem,local,ip,full_array,size);
1990       for (int i=0; i<size; i++)
1991         element_array[ip[i]].push_back(local[i]);
1992       delete[] full_array;
1993       break;
1994     case MED_EN::MED_NODE :
1995       old_collection.getTopology()->convertNodeToGlobal(idomain,list_of_elems,nbelem,array);
1996       _topology->convertGlobalNodeListWithTwins(array,nbelem,local,ip,full_array,size);
1997       for (int i=0; i<size; i++)
1998         element_array[ip[i]].push_back(local[i]);
1999       delete[] full_array;
2000       break;
2001
2002     }
2003     delete[] ip;
2004     delete[] local;
2005     delete[] array;
2006
2007     if (support->isOnAllElements()) delete[] list_of_elems;
2008   }
2009
2010   //creating new groups from the previous list of elements
2011   for (int idomain=0; idomain <_topology->nbDomain(); idomain++)
2012   {
2013     if ( _mesh[idomain]->getNumberOfNodes() < 1 ||
2014          (element_array[idomain].empty() && !_create_empty_groups))
2015     {
2016       new_support[idomain]->removeReference();
2017       new_support[idomain]=0;
2018       continue;
2019     }
2020     MEDMEM::SUPPORT* support= new_support[idomain];
2021     support->setName(name);
2022     support->setMesh(_mesh[idomain]);
2023     support->setDescription(description);
2024     support->setEntity(entity);
2025
2026     if ( !element_array[idomain].empty() ) /* if() was added for issue 0021576
2027                                               to prevent creation of faces */
2028       {
2029         element_array[idomain].sort();
2030         element_array[idomain].unique();
2031
2032         if (entity != MED_EN::MED_NODE)
2033           support->fillFromElementList(element_array[idomain]);
2034         else
2035           support->fillFromNodeList(element_array[idomain]);
2036       }
2037   }
2038 }
2039
2040 void MESHCollection::castField(const MESHCollection& old_collection, const string& fieldname, int itnumber, int ordernumber)
2041 {
2042   int type=old_collection.getDriver()->getFieldType(fieldname);
2043   char field_char[80];
2044   strcpy(field_char,fieldname.c_str());
2045
2046   if (type ==0)
2047     castFields<int>(old_collection, field_char, itnumber, ordernumber);
2048   else
2049     castFields<double>(old_collection, field_char, itnumber, ordernumber);
2050 }
2051
2052 void MESHCollection::castAllFields(const MESHCollection& initial_collection)
2053 {
2054   vector <string> field_names;
2055   vector <int> iternumber;
2056   vector <int> ordernumber;
2057   vector <int> types;
2058   initial_collection.getDriver()->readFileStruct(field_names,iternumber,ordernumber,types);
2059
2060   for (unsigned i=0; i<field_names.size(); i++)
2061   {
2062     char field_char[80];
2063     strcpy(field_char,field_names[i].c_str());
2064
2065     // choosing whether the field is of int or double type
2066     if (types[i] ==0)
2067       castFields<int>(initial_collection, field_char, iternumber[i], ordernumber[i]);
2068     else
2069       castFields<double>(initial_collection, field_char, iternumber[i], ordernumber[i]);
2070   }
2071 }
2072
2073 void MESHCollection::createNodalConnectivity(const MESHCollection& initial_collection,int idomain, MED_EN::medEntityMesh entity)
2074 {
2075   MESSAGE_MED ("beginning of createNodalConnectivity for entity "<<entity);
2076   int dimension=0;
2077   int nb_elems=0;
2078   MEDMEM::MESHING* mesh_builder = static_cast<MEDMEM::MESHING*>(_mesh[idomain]);
2079
2080
2081   //number of elements per type
2082   std::map<MED_EN::medGeometryElement,int> type_numbers;
2083
2084   //creating arrays for storing global numbers and cell types
2085   switch (entity)
2086   {
2087   case MED_EN::MED_CELL:
2088     dimension=initial_collection.getMeshDimension();
2089     nb_elems=_topology->getCellNumber(idomain);
2090     break;
2091   case MED_EN::MED_EDGE:
2092   case MED_EN::MED_FACE:
2093     dimension=initial_collection.getMeshDimension()-1;
2094     nb_elems=_topology->getFaceNumber(idomain);
2095     break;
2096   default:
2097     nb_elems=0;
2098     break;
2099   }
2100
2101   if (nb_elems == 0) return;
2102   SCRUTE_MED(nb_elems);
2103
2104
2105   int *list= new int[nb_elems];
2106   MED_EN::medGeometryElement *cell_type_list= new MED_EN::medGeometryElement[nb_elems];
2107
2108
2109   //  cout << "Beginning of retrieval "<<endl;
2110   //retrieving global id list
2111   switch (entity)
2112   {
2113   case MED_EN::MED_CELL:
2114     _topology->getCellList(idomain,list);
2115     break;
2116   case MED_EN::MED_EDGE:
2117   case MED_EN::MED_FACE:
2118     _topology->getFaceList(idomain,list);
2119     break;
2120   default:
2121
2122     break;
2123   }
2124
2125   //retrieving cell_types
2126   initial_collection.getTypeList(list,nb_elems,entity,cell_type_list);
2127   //  cout <<"end of type retrieval"<<endl;
2128   //vector containing the number of cells per type
2129   type_numbers.clear();
2130   for (int icell=0; icell<nb_elems; icell++)
2131   {
2132     map<MED_EN::medGeometryElement,int>::iterator iter= type_numbers.find(cell_type_list[icell]);
2133     if (iter!=type_numbers.end())
2134       (iter->second)++;
2135     else
2136       type_numbers[cell_type_list[icell]]=1;
2137
2138   }
2139   //cout << "Nombre de tetras"<<type_numbers[304]<<endl;
2140   int nb_present_types=type_numbers.size();
2141
2142   //setting the list of cells for each type
2143   map<MED_EN::medGeometryElement,int> index;
2144
2145   map<MED_EN::medGeometryElement,int*> type_cell_list;
2146
2147   MED_EN::MESH_ENTITIES::const_iterator currentEntity;
2148   std::map<MED_EN::medGeometryElement,int>::const_iterator iter;
2149   //currentEntity  = MED_EN::meshEntities.find(entity);
2150   for (iter = type_numbers.begin();iter != type_numbers.end(); iter++)
2151   {
2152     MED_EN::medGeometryElement type = iter->first;
2153     if (!isDimensionOK(type,dimension)) continue;
2154     //if (iter->second==0) continue;
2155     index[type]=0;
2156     type_cell_list[type]=new int[type_numbers[type]];
2157     // cout << "type :"<<type<<" nb:"<<type_numbers[type]<<endl;
2158   }
2159
2160   for (int icell=0; icell<nb_elems; icell++)
2161   {
2162     type_cell_list[cell_type_list[icell]][index[cell_type_list[icell]]++]=list[icell];
2163   }
2164
2165   delete[]list;
2166   delete[]cell_type_list;
2167
2168   //setting the list of present types
2169   int* present_type_numbers=new int[nb_present_types];
2170   MED_EN::medGeometryElement* type_array = new MED_EN::medGeometryElement[nb_present_types];
2171   MESSAGE_MED("Nb de types presents "<<nb_present_types);
2172   int itype=0;
2173   for (iter = type_numbers.begin();iter != type_numbers.end(); iter++)
2174   {
2175     MED_EN::medGeometryElement type = iter->first;
2176     if (!isDimensionOK(type,dimension)) continue;
2177
2178     type_array[itype]=type;
2179
2180     present_type_numbers[itype]=type_numbers[type];
2181
2182     MESSAGE_MED("Nombre d'elements de type "<<type<<" : "<<type_numbers[type]);
2183     itype++;
2184   }
2185
2186   //retrieving connectivity in global numbering for each type
2187   map<MED_EN::medGeometryElement,int*> type_connectivity;
2188   vector<int> polygon_conn;
2189   vector<int> polygon_conn_index;
2190   vector<int> polyhedron_conn;
2191   vector<int> polyhedron_conn_index;
2192   vector<int> polyhedron_face_index;
2193
2194   //Treating nodes
2195
2196
2197   for (iter = type_numbers.begin();iter != type_numbers.end(); iter++)
2198   {
2199     MED_EN::medGeometryElement type = iter->first;
2200
2201
2202     if (!isDimensionOK(type,dimension)) continue;
2203     if (type != MED_EN::MED_POLYGON && type != MED_EN::MED_POLYHEDRA)
2204     {
2205       int nbnode_per_type = (int)type%100;
2206       type_connectivity[type]=new int[type_numbers[type]*nbnode_per_type];
2207       initial_collection.getNodeConnectivity(type_cell_list[type],type_numbers[type],entity,type,type_connectivity[type]);
2208     }
2209     else if (type == MED_EN::MED_POLYGON && dimension==2)
2210     {
2211       initial_collection.getPolygonNodeConnectivity(type_cell_list[type],type_numbers[type],entity,polygon_conn,polygon_conn_index);
2212     }
2213     else if (type == MED_EN::MED_POLYHEDRA && dimension==3)
2214     {
2215       initial_collection.getPolyhedraNodeConnectivity(type_cell_list[type],type_numbers[type],entity,polyhedron_conn,polyhedron_conn_index);
2216     }
2217     delete[] type_cell_list[type];
2218   }
2219
2220   //creating node mapping
2221   //!TODO : compute the total number of nodes
2222   if (entity==MED_EN::MED_CELL)
2223   {
2224     _topology->createNodeMapping(type_connectivity,type_numbers,polygon_conn,polygon_conn_index,
2225                                  polyhedron_conn,polyhedron_conn_index,polyhedron_face_index,idomain);
2226   }
2227
2228   //converting node global numberings to local numberings
2229   //for (iter = (*currentEntity).second.begin();iter != (*currentEntity).second.end(); iter++)
2230   for (iter = type_numbers.begin();iter != type_numbers.end(); iter++)
2231   {
2232     MED_EN::medGeometryElement type = iter->first;
2233
2234     if (!isDimensionOK(type, dimension)) continue;
2235     if (type_numbers[type]==0) continue;
2236     if (type != MED_EN::MED_POLYGON && type != MED_EN::MED_POLYHEDRA)
2237     {
2238       int nbnode_per_type = (int)type%100;
2239       _topology->convertToLocal2ndVersion(type_connectivity[type],type_numbers[type]*nbnode_per_type,idomain);
2240     }
2241     else if (type == MED_EN::MED_POLYGON && dimension==2)
2242     {
2243       int nbpoly = type_numbers[type];
2244       _topology->convertToLocal2ndVersion(&polygon_conn[0], polygon_conn_index[nbpoly]-1, idomain);
2245     }
2246     else if (type == MED_EN::MED_POLYHEDRA && dimension==3)
2247     {
2248       int nbpoly = type_numbers[type];
2249       _topology->convertToLocal2ndVersion(&polyhedron_conn[0], polyhedron_face_index[polyhedron_conn_index[nbpoly]-1]-1, idomain);
2250     }
2251
2252   }
2253
2254
2255   //writing coordinates
2256   if (entity==MED_EN::MED_CELL)
2257   {
2258     //setting coordinates from initial_collection coordinates
2259     int nbnode=_topology->getNodeNumber(idomain);
2260     MESSAGE_MED("Number of nodes on domain "<< idomain <<" : "<<nbnode);
2261
2262     double* coordinates=new double[initial_collection.getSpaceDimension()*nbnode];
2263     int* node_list=new int[nbnode];
2264     _topology->getNodeList(idomain,node_list);
2265     initial_collection.getCoordinates(node_list,nbnode,coordinates);
2266     delete[] node_list;
2267
2268     // redundant specification of number of nodes is required!! MED imperfection, sorry...
2269
2270     //TODO : change MEDMEM so that it accepts a direct setting of coordinates
2271     // (in the present version, it is deep-copied)
2272     mesh_builder->setCoordinates(initial_collection.getSpaceDimension(),
2273                                  nbnode, coordinates, initial_collection.getSystem(),
2274                                  MED_EN::MED_FULL_INTERLACE);
2275     delete [] coordinates;
2276   }
2277
2278   int nb_plain_types=0;
2279   for (iter = type_numbers.begin();iter != type_numbers.end(); iter++)
2280   {
2281     MED_EN::medGeometryElement type = iter->first;
2282
2283     if (!isDimensionOK(type, dimension)) continue;
2284     if (type_numbers[type]==0) continue;
2285     nb_plain_types++;
2286   }
2287   mesh_builder->setNumberOfTypes(nb_plain_types,entity);
2288   mesh_builder->setTypes(type_array,entity);
2289   mesh_builder->setNumberOfElements(present_type_numbers,entity);
2290
2291   delete[]present_type_numbers;
2292   delete[]type_array;
2293   //setting node connectivities
2294   for (iter = type_numbers.begin();iter != type_numbers.end(); iter++)
2295   {
2296     MED_EN::medGeometryElement type = iter->first;
2297
2298     if (!isDimensionOK(type,dimension)) continue;
2299     if (type_numbers[type]==0) continue;
2300
2301     if (type != MED_EN::MED_POLYHEDRA && type != MED_EN::MED_POLYGON)
2302     {
2303       mesh_builder->setConnectivity(entity,type,type_connectivity[type]);
2304       delete[] type_connectivity[type];
2305     }
2306     else if (type == MED_EN::MED_POLYGON && dimension ==2)
2307     {
2308       mesh_builder->setConnectivity(entity,type,&polygon_conn[0],&polygon_conn_index[0]);
2309     }
2310     else if (type == MED_EN::MED_POLYHEDRA && dimension ==3)
2311     {
2312       mesh_builder->setConnectivity(entity,type,&polyhedron_conn[0],&polyhedron_conn_index[0]);
2313     }
2314   }
2315   MESSAGE_MED("end of createNodalConnectivity");
2316 }
2317
2318
2319 /*! retrieves the faces that are present in a mesh and stores them in a
2320  * dynamic structure made of a map of MEDSPLITTER_FaceModel
2321  *
2322  * \param idomain domain id on which the faces are collected
2323  * \param face_map container storing the faces
2324  */
2325 void MESHCollection::getFaces(int idomain,
2326                               map<MED_EN::medGeometryElement, vector<MEDSPLITTER_FaceModel*> >& face_map)
2327 {
2328   MED_EN::medEntityMesh constituent_entity = getSubEntity();
2329   const medGeometryElement* types;
2330   try
2331   {
2332     types = _mesh[idomain]->getTypes(constituent_entity);
2333     if ( !types ) return;
2334   }
2335   catch(MEDEXCEPTION&){ return;}
2336
2337   int nbtypes  = _mesh[idomain]->getNumberOfTypes(constituent_entity);
2338   const int* global_numbering= _mesh[idomain]->getGlobalNumberingIndex(constituent_entity);
2339   int* conn = const_cast<int*> (_mesh[idomain]->getConnectivity(MED_EN::MED_NODAL,constituent_entity, MED_EN::MED_ALL_ELEMENTS));
2340   for (int itype=0; itype<nbtypes; itype++)
2341   {
2342     for (int iface=global_numbering[itype]; iface<global_numbering[itype+1]; iface++)
2343     {
2344       MEDSPLITTER_FaceModel* face_model = new MEDSPLITTER_FaceModel();
2345       MED_EN::medGeometryElement type =  types[itype];
2346       face_model->setType(type);
2347       int nbnodes = type%100;
2348       face_model->setNbNodes(nbnodes);
2349       face_model->setGlobal(_topology->convertFaceToGlobal(idomain,iface));
2350       for (int i=0; i<nbnodes; i++)
2351       {
2352         (*face_model)[i]=*conn++;
2353       }
2354       face_map[type].push_back(face_model);
2355     }
2356   }
2357 }
2358
2359 /*! retrieves the face that is common to two cells located on two different processors
2360  *
2361  * \param ip1 domain id for cell 1
2362  * \param ilocal1 cell id for cell 1
2363  * \param ip2 domain id for cell 2
2364  * \param ilocal2 cell id for cell 2
2365  * \param face_index global index for the newly created face
2366  */
2367 MEDSPLITTER_FaceModel* MESHCollection::getCommonFace(int ip1,int ilocal1,int ip2,int ilocal2,int face_index)
2368 {
2369   MED_EN::medGeometryElement type1 = _mesh[ip1]->getElementType(MED_EN::MED_CELL,ilocal1);
2370   MEDMEM::CELLMODEL celltype1 (type1);
2371
2372   const int* conn_index1 =  _mesh[ip1]->getConnectivityIndex(MED_EN::MED_NODAL,MED_EN::MED_CELL);
2373   const int* conn1 = _mesh[ip1]->getConnectivity(MED_EN::MED_NODAL,MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
2374
2375   // MED_EN::medGeometryElement type2 = _mesh[ip2]->getElementType(MED_EN::MED_CELL,ilocal2);
2376   //MEDMEM::CELLTYPE celltype2 (type2);
2377   const int* conn_index2 =  _mesh[ip2]->getConnectivityIndex(MED_EN::MED_NODAL,MED_EN::MED_CELL);
2378   const int* conn2 = _mesh[ip2]->getConnectivity(MED_EN::MED_NODAL,MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
2379
2380   vector<int> nodes1, nodes1_local;
2381   vector<int> nodes2;
2382   for (int i=  conn_index1[ilocal1-1]; i<conn_index1[ilocal1]; i++)
2383   {
2384     nodes1.push_back(_topology->convertNodeToGlobal(ip1,*(conn1+i-1)));
2385     nodes1_local.push_back( conn1[i-1] );
2386   }
2387   for (int i=  conn_index2[ilocal2-1]; i<conn_index2[ilocal2]; i++)
2388     nodes2.push_back(_topology->convertNodeToGlobal(ip2,*(conn2+i-1)));
2389
2390   return MEDSPLITTER_FaceModel::getCommonFace( &nodes1[0], &nodes1_local[0], celltype1,
2391                                                &nodes2[0], nodes2.size(),  face_index);
2392 }
2393
2394 //================================================================================
2395 /*!
2396  * \brief Makes a face common for two given cells
2397  *  \param nodes1 - globl nodes of the first cell
2398  *  \param nodes1_local - local nodes of the first cell
2399  *  \param celltype1 - cell model of the first cell
2400  *  \param nodes2 - globl nodes of the second cell
2401  *  \param nb_nodes2 - nb of nodes of the second cell
2402  *  \param global_id - id of the new common face
2403  */
2404 //================================================================================
2405
2406 MEDSPLITTER_FaceModel*
2407 MEDSPLITTER_FaceModel::getCommonFace(const int*               nodes1,
2408                                      const int*               nodes1_local,
2409                                      const MEDMEM::CELLMODEL& celltype1,
2410                                      const int*               nodes2,
2411                                      int                      nb_nodes2,
2412                                      int                      global_id)
2413 {
2414   int nbfaces= celltype1.getNumberOfConstituents(1);
2415   int ** faces = celltype1.getConstituents(1);
2416   MED_EN::medGeometryElement* types = celltype1.getConstituentsType(1);
2417   int iface=0;
2418   int dimension=celltype1.getDimension();
2419
2420   while (iface<nbfaces)
2421   {
2422     //SCRUTE_MED (iface);
2423     int nbnodes= types[iface]%100;
2424     const int* nodes = celltype1.getNodesConstituent(1,iface+1);
2425     int common_nodes=0;
2426     for (int i=0; i<nbnodes;i++)
2427     {
2428       for (int i2=0; i2<nb_nodes2; i2++)
2429       {
2430         if (nodes1[nodes[i]-1]==nodes2[i2]) common_nodes++;
2431       }
2432     }
2433     if (common_nodes>=dimension) break;
2434     iface++;
2435   }
2436
2437   if (iface==nbfaces)
2438     throw MEDEXCEPTION("MEDSPLITTER::getCommonFace - No common face found !");
2439
2440   MEDSPLITTER_FaceModel* face_model = new MEDSPLITTER_FaceModel;
2441   face_model->setType(types[iface]);
2442   int nbnodes = types[iface]%100;
2443   face_model->setNbNodes(nbnodes);
2444   face_model->setGlobal(global_id);
2445   for (int i=0; i<nbnodes; i++)
2446     (*face_model)[i]=nodes1_local[faces[iface][i]-1];
2447
2448   return face_model;
2449 }