src/MEDSPLITTER/MEDSPLITTER_ParallelTopology.cxx

   1 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 #include <set>
  20 #include <map>
  21 #include <vector>
  22
  23 #include "InterpKernelHashMap.hxx"
  24
  25 #include "MEDMEM_CellModel.hxx"
  26 #include "MEDMEM_ConnectZone.hxx"
  27 #include "MEDMEM_DriversDef.hxx"
  28 #include "MEDMEM_Exception.hxx"
  29 #include "MEDMEM_Mesh.hxx"
  30 #include "MEDMEM_MeshFuse.hxx"
  31 #include "MEDMEM_SkyLineArray.hxx"
  32 #include "MEDMEM_Utilities.hxx"
  33
  34 #include "MEDSPLITTER_MESHCollection.hxx"
  35 #include "MEDSPLITTER_Topology.hxx"
  36 #include "MEDSPLITTER_Graph.hxx"
  37 #include "MEDSPLITTER_ParallelTopology.hxx"
  38
  39 using namespace INTERP_KERNEL;
  40
  41 using namespace MEDSPLITTER;
  42
  43 //empty constructor
  44 ParallelTopology::ParallelTopology():m_nb_domain(0),m_mesh_dimension(0)
  45 {}
  46
  47 //!constructing topology according to mesh collection
  48 ParallelTopology::ParallelTopology(const vector<MEDMEM::MESH*>& meshes,
  49                                    const vector<MEDMEM::CONNECTZONE*>& cz,
  50                                    vector<int*>& cellglobal,
  51                                    vector<int*>& nodeglobal,
  52                                    vector<int*>& faceglobal):m_nb_domain(meshes.size())/*,m_mesh_dimension(meshes[0]->getMeshDimension())*/
  53 {
  54
  55   int index_global=0;
  56   int index_node_global=0;
  57   int index_face_global=0;
  58
  59   m_nb_cells.resize(m_nb_domain);
  60   m_nb_nodes.resize(m_nb_domain);
  61   m_nb_faces.resize(m_nb_domain);
  62
  63   m_loc_to_glob.resize(m_nb_domain);
  64   m_node_loc_to_glob.resize(m_nb_domain);
  65   m_face_loc_to_glob.resize(m_nb_domain);
  66
  67   MED_EN::medEntityMesh constituent_entity;
  68
  69   bool parallel_mode = false;
  70   for (int idomain=0; !parallel_mode && idomain<m_nb_domain; idomain++)
  71     parallel_mode = (!meshes[idomain]);
  72
  73   for (int idomain=0; idomain<m_nb_domain; idomain++)
  74   {
  75     if ( !meshes[idomain] )
  76       continue;
  77     m_mesh_dimension = meshes[idomain]->getMeshDimension();
  78     constituent_entity = (m_mesh_dimension == 3 ? MED_EN::MED_FACE : MED_EN::MED_EDGE );
  79
  80     //creating cell maps
  81     m_nb_cells[idomain]=meshes[idomain]->getNumberOfElements(MED_EN::MED_CELL, MED_EN::MED_ALL_ELEMENTS);
  82     //              cout << "Nb cells (domain "<<idomain<<") = "<<m_nb_cells[idomain];
  83     m_loc_to_glob[idomain].resize(m_nb_cells[idomain]);
  84
  85     if (cellglobal[idomain]==0 || parallel_mode)
  86     {
  87       MESSAGE_MED("Creating global numbering");
  88       //creating global numbering from scratch
  89       for (int i=0; i<m_nb_cells[idomain]; i++)
  90       {
  91         index_global++;
  92         m_glob_to_loc[index_global]=make_pair(idomain,i+1);
  93         //m_loc_to_glob[make_pair(idomain,i+1)]=index_global;
  94         m_loc_to_glob[idomain][i]=index_global;
  95         //        cout<<"glob:"<<index_global<<" --> ("<<idomain<<","<<i+1<<")"<<endl;
  96       }
  97     }
  98     //using global numbering coming from a previous numbering
  99     else
 100     {
 101       MESSAGE_MED("Using former global numbering");
 102       for (int i=0; i<m_nb_cells[idomain]; i++)
 103       {
 104         int global=cellglobal[idomain][i];
 105         m_glob_to_loc[global]=make_pair(idomain,i+1);
 106         //m_loc_to_glob[make_pair(idomain,i+1)]=global;
 107         m_loc_to_glob[idomain][i]=global;
 108         index_global++;
 109         //        cout<<"glob:"<<global<<" --> ("<<idomain<<","<<i+1<<")"<<endl;
 110       }
 111     }
 112
 113     //cas sequentiel
 114     if (m_nb_domain==1)
 115     {
 116       m_nb_total_cells=index_global;
 117       m_nb_cells[0]=index_global;
 118       m_node_loc_to_glob[idomain].resize(meshes[idomain]->getNumberOfNodes());
 119       for (int i=0; i<meshes[idomain]->getNumberOfNodes(); i++)
 120       {
 121         m_node_glob_to_loc.insert(make_pair(i+1,make_pair(0,i+1)));
 122         //m_node_loc_to_glob.insert(make_pair(make_pair(0,i+1), i+1));
 123         m_node_loc_to_glob[0][i]=i+1;
 124       }
 125       m_nb_total_nodes=meshes[idomain]->getNumberOfNodes();
 126       m_nb_nodes[0]=m_nb_total_nodes;
 127
 128       //                      meshes[idomain]->getConnectivity( MED_EN::MED_DESCENDING, MED_EN::MED_CELL, MED_EN::MED_ALL_ELEMENTS);
 129       int nbfaces=meshes[idomain]->getNumberOfElements(constituent_entity,MED_EN::MED_ALL_ELEMENTS);
 130       m_face_loc_to_glob[idomain].resize(nbfaces);
 131       for (int i=0; i<nbfaces; i++)
 132       {
 133         m_face_glob_to_loc.insert(make_pair(i+1,make_pair(0,i+1)));
 134         //m_face_loc_to_glob.insert(make_pair(make_pair(0,i+1), i+1));
 135         m_face_loc_to_glob[0][i]=i+1;
 136       }
 137       m_nb_total_faces=nbfaces;
 138       m_nb_faces[0]=nbfaces;
 139       MESSAGE_MED ("nb total cells "<< m_nb_total_cells);
 140       MESSAGE_MED("nb total nodes "<< m_nb_total_nodes);
 141       MESSAGE_MED("nb total faces "<< m_nb_total_faces);
 142       return;
 143     }
 144
 145     //creating node maps
 146     m_nb_nodes[idomain]=meshes[idomain]->getNumberOfNodes();
 147     INTERP_KERNEL::HashMap <int,pair<int,int> > local2distant;
 148     m_node_loc_to_glob[idomain].resize(m_nb_nodes[idomain]);
 149     for (unsigned icz=0; icz<cz.size(); icz++)
 150     {
 151       if (cz[icz]->getLocalDomainNumber() == idomain &&
 152           cz[icz]->getLocalDomainNumber()>cz[icz]->getDistantDomainNumber())
 153       {
 154         int nb_node= cz[icz]->getNodeNumber();
 155         const int* node_corresp=cz[icz]->getNodeCorrespValue();
 156         int distant_ip = cz[icz]->getDistantDomainNumber();
 157         for (int i=0; i< nb_node; i++)
 158         {
 159           int local= node_corresp[i*2];
 160           int distant = node_corresp[i*2+1];
 161           local2distant.insert(make_pair(local, make_pair(distant_ip,distant)));
 162         }
 163       }
 164     }
 165     // setting mappings for all nodes
 166     if (nodeglobal[idomain]==0)
 167     {
 168       for (int inode=0; inode<m_nb_nodes[idomain]; inode++)
 169       {
 170         if (local2distant.find(inode+1)==local2distant.end())
 171         {
 172           index_node_global++;
 173           m_node_glob_to_loc.insert(make_pair(index_node_global,make_pair(idomain,inode+1)));
 174           //m_node_loc_to_glob[make_pair(idomain,inode+1)]=index_node_global;
 175           m_node_loc_to_glob[idomain][inode]=index_node_global;
 176         }
 177         else
 178         {
 179           int ip = (local2distant.find(inode+1)->second).first;
 180           int distant = (local2distant.find(inode+1)->second).second;
 181           //int global_number=m_loc_to_glob[make_pair(ip,distant)];
 182           int global_number=m_loc_to_glob[ip][distant-1];
 183           m_node_glob_to_loc.insert(make_pair(global_number,make_pair(idomain,inode+1)));
 184           //m_node_loc_to_glob[make_pair(idomain,inode+1)]=global_number;
 185           m_node_loc_to_glob[idomain][inode]=global_number;
 186         }
 187       }
 188     }
 189     //using former node numbering
 190     else
 191     {//       cout << "("<<idomain<<","<<i+1<<")->"<<i+1<<endl;
 192       for (int inode=0; inode<m_nb_nodes[idomain]; inode++)
 193       {
 194         int global_number=nodeglobal[idomain][inode];
 195         //        cout << "global_number "<<global_number<<endl;
 196         m_node_glob_to_loc.insert(make_pair(global_number,make_pair(idomain,inode+1)));
 197         //m_node_loc_to_glob[make_pair(idomain,inode+1)]=global_number;
 198         m_node_loc_to_glob[idomain][inode]=global_number;
 199       }
 200     }
 201
 202
 203     //creating  dimension d-1 component mappings
 204
 205     //              meshes[idomain]->getConnectivity( MED_EN::MED_DESCENDING, MED_EN::MED_CELL, MED_EN::MED_ALL_ELEMENTS);
 206     m_nb_faces[idomain]=meshes[idomain]->getNumberOfElements(constituent_entity,MED_EN::MED_ALL_ELEMENTS);
 207     MESSAGE_MED ("Nb faces domain " << idomain<<m_nb_faces[idomain]);
 208     m_face_loc_to_glob[idomain].resize(m_nb_faces[idomain]);
 209     local2distant.clear();
 210     for (unsigned icz=0; icz<cz.size(); icz++)
 211     {
 212       if (cz[icz]->getLocalDomainNumber() == idomain &&
 213           cz[icz]->getLocalDomainNumber()>cz[icz]->getDistantDomainNumber())
 214       {
 215         int nb_face= cz[icz]->getFaceNumber();
 216         const int* face_corresp=cz[icz]->getFaceCorrespValue();
 217         int distant_ip = cz[icz]->getDistantDomainNumber();
 218         for (int i=0; i< nb_face; i++)
 219         {
 220           int local= face_corresp[i*2];
 221           int distant = face_corresp[i*2+1];
 222           local2distant.insert(make_pair(local, make_pair(distant_ip,distant)));
 223         }
 224       }
 225     }
 226     // setting mappings for all faces
 227     if (faceglobal[idomain]==0)
 228     {
 229       for (int iface=0; iface<m_nb_faces[idomain]; iface++)
 230       {
 231         if (local2distant.find(iface+1)==local2distant.end())
 232         {
 233           index_face_global++;
 234           m_face_glob_to_loc.insert(make_pair(index_face_global,make_pair(idomain,iface+1)));
 235           //m_face_loc_to_glob[make_pair(idomain,iface+1)]=index_face_global;
 236           m_face_loc_to_glob[idomain][iface]=index_face_global;
 237         }
 238         else
 239         {
 240           int ip = (local2distant.find(iface+1)->second).first;
 241           int distant = (local2distant.find(iface+1)->second).second;
 242           //int global_number=m_loc_to_glob[make_pair(ip,distant)];
 243           int global_number=m_loc_to_glob[ip][distant-1];
 244           m_face_glob_to_loc.insert(make_pair(global_number,make_pair(idomain,iface+1)));
 245           //m_face_loc_to_glob[make_pair(idomain,iface+1)]=global_number;
 246           m_face_loc_to_glob[idomain][iface]=global_number;
 247         }
 248       }
 249     }
 250     //using former face numbering
 251     else
 252     {
 253       for (int iface=0; iface<m_nb_faces[idomain]; iface++)
 254       {
 255         int global_number=faceglobal[idomain][iface];
 256         //SCRUTE_MED(global_number);
 257         m_face_glob_to_loc.insert(make_pair(global_number,make_pair(idomain,iface+1)));
 258         //m_face_loc_to_glob[make_pair(idomain,iface+1)]=global_number;
 259         m_face_loc_to_glob[idomain][iface]=global_number;
 260       }
 261     }
 262   }
 263
 264   m_nb_total_cells=index_global;
 265   m_nb_total_nodes=index_node_global;
 266   m_nb_total_faces=index_face_global;
 267   SCRUTE_MED(m_nb_total_cells);
 268   SCRUTE_MED(m_nb_total_faces);
 269   SCRUTE_MED(m_nb_total_nodes);
 270
 271 }
 272
 273
 274 //!constructing ParallelTopology from an old topology and a graph
 275 ParallelTopology::ParallelTopology(boost::shared_ptr<Graph> graph, int nb_domain, int mesh_dimension):
 276   m_nb_domain(nb_domain),
 277   m_mesh_dimension(mesh_dimension),
 278   m_nb_cells(graph->nbVertices()),
 279   m_graph(graph)
 280 {
 281   m_nb_cells.resize(m_nb_domain);
 282   m_nb_nodes.resize(m_nb_domain);
 283   m_nb_faces.resize(m_nb_domain);
 284
 285   m_loc_to_glob.resize(m_nb_domain);
 286   m_node_loc_to_glob.resize(m_nb_domain);
 287   m_face_loc_to_glob.resize(m_nb_domain);
 288
 289   // used in parallel mode only
 290   m_cell_loc_to_glob_fuse.resize(m_nb_domain);
 291   m_face_loc_to_glob_fuse.resize(m_nb_domain);
 292
 293   for (int i=0; i<m_nb_domain; i++)
 294     m_nb_cells[i]=0;
 295
 296   const int* part = graph-> getPart();
 297   m_nb_total_cells= graph->nbVertices();
 298
 299   for (int icell=0; icell<m_nb_total_cells; icell++)
 300   {
 301     int idomain = part[icell];
 302     m_nb_cells[idomain]++;
 303     //m_loc_to_glob[make_pair(idomain,m_nb_cells[idomain])]=icell+1;
 304     m_loc_to_glob[idomain].push_back(icell+1);
 305     m_glob_to_loc[icell+1]=make_pair(idomain,m_nb_cells[idomain]);
 306
 307   }
 308   for (int idomain=0; idomain<m_nb_domain; idomain++)
 309     MESSAGE_MED("Nombre de cellules dans le domaine "<< idomain <<" : "<<m_nb_cells[idomain]);
 310
 311   SCRUTE_MED(m_nb_total_cells);
 312
 313 }
 314
 315 ParallelTopology::~ParallelTopology()
 316 {
 317
 318 }
 319
 320 /*!Converts a list of global node numbers
 321  * to a distributed array with local cell numbers.
 322  *
 323  * If a node in the list is represented on several domains,
 324  * only the first value is returned
 325  * */
 326 void ParallelTopology::convertGlobalNodeList(const int* node_list, int nbnode, int* local, int* ip)
 327 {
 328   if (m_node_glob_to_loc.empty())
 329     throw MEDMEM::MEDEXCEPTION("convertGlobalNodeList - Node mapping has not yet been built");
 330   for (int i=0; i< nbnode; i++)
 331   {
 332     pair<int,int> local_node = m_node_glob_to_loc.find(node_list[i])->second;
 333     ip[i]=local_node.first;
 334     local[i]=local_node.second;
 335   }
 336 }
 337
 338 /*!Converts a list of global node numbers on domain ip
 339  * to a distributed array with local cell numbers.
 340  *
 341  * If a node in the list is represented on several domains,
 342  * only the value with domain ip is returned
 343  *
 344  * */
 345 void ParallelTopology::convertGlobalNodeList(const int* node_list, int nbnode, int* local, int ip)
 346 {
 347   if (m_node_glob_to_loc.empty())
 348     throw MEDMEM::MEDEXCEPTION("convertGlobalNodeList - Node mapping has not yet been built");
 349
 350   for (int i=0; i< nbnode; i++)
 351   {
 352     typedef INTERP_KERNEL::HashMultiMap<int,pair<int,int> >::iterator mmiter;
 353     pair<mmiter,mmiter> range=m_node_glob_to_loc.equal_range(node_list[i]);
 354     for (mmiter it=range.first; it !=range.second; it++)
 355     {
 356       int ipfound=(it->second).first;
 357       if (ipfound==ip)
 358         local[i]=(it->second).second;
 359     }
 360   }
 361 }
 362
 363 /*!Converts a list of global node numbers
 364  * to a distributed array with local cell numbers.
 365  *
 366  * If a node in the list is represented on several domains,
 367  * all the values are put in the array
 368  * */
 369 void ParallelTopology::convertGlobalNodeListWithTwins(const int* node_list, int nbnode, int*& local, int*& ip,int*& full_array, int& size)
 370 {
 371   if (m_node_glob_to_loc.empty())
 372     throw MEDMEM::MEDEXCEPTION("convertGlobalNodeList - Node mapping has not yet been built");
 373
 374   size=0;
 375   for (int i=0; i< nbnode; i++)
 376   {
 377     int count= m_node_glob_to_loc.count(node_list[i]);
 378     //      if (count > 1)
 379     //        cout << "noeud " << node_list[i]<< " doublon d'ordre " << count<<endl;
 380     size+=count;
 381   }
 382   int index=0;
 383   ip=new int[size];
 384   local=new int[size];
 385   full_array=new int[size];
 386   for (int i=0; i< nbnode; i++)
 387   {
 388     typedef INTERP_KERNEL::HashMultiMap<int,pair<int,int> >::iterator mmiter;
 389     pair<mmiter,mmiter> range=m_node_glob_to_loc.equal_range(node_list[i]);
 390     for (mmiter it=range.first; it !=range.second; it++)
 391     {
 392       ip[index]=(it->second).first;
 393       local[index]=(it->second).second;
 394       full_array [index]=node_list[i];
 395       index++;
 396     }
 397
 398   }
 399 }
 400
 401 /*!Converts a list of global face numbers
 402  * to a distributed array with local face numbers.
 403  *
 404  * If a face in the list is represented on several domains,
 405  * all the values are put in the array
 406  * */
 407 void ParallelTopology::convertGlobalFaceListWithTwins(const int* face_list, int nbface, int*& local, int*& ip, int*& full_array,int& size)
 408 {
 409   size=0;
 410   for (int i=0; i< nbface; i++)
 411   {
 412     //int count = m_face_glob_to_loc.count(face_list[i]);
 413     //if (count >1) MESSAGE_MED("face en doublon "<<face_list[i]);
 414     size+= m_face_glob_to_loc.count(face_list[i]);
 415   }
 416   int index=0;
 417   ip=new int[size];
 418   local=new int[size];
 419   full_array=new int[size];
 420   for (int i=0; i< nbface; i++)
 421   {
 422     typedef INTERP_KERNEL::HashMultiMap<int,pair<int,int> >::iterator mmiter;
 423     pair<mmiter,mmiter> range=m_face_glob_to_loc.equal_range(face_list[i]);
 424     for (mmiter it=range.first; it !=range.second; it++)
 425     {
 426       ip[index]=(it->second).first;
 427       local[index]=(it->second).second;
 428       full_array[index]=face_list[i];
 429       index++;
 430     }
 431
 432   }
 433 }
 434
 435 //!converts a list of global cell numbers
 436 //!to a distributed array with local cell numbers
 437 void ParallelTopology::convertGlobalCellList(const int* cell_list, int nbcell, int* local, int* ip)
 438 {
 439   for (int i=0; i< nbcell; i++)
 440   {
 441     INTERP_KERNEL::HashMap<int, pair<int,int> >::const_iterator iter = m_glob_to_loc.find(cell_list[i]);
 442     ip[i]=(iter->second).first;
 443     local[i]=(iter->second).second;
 444   }
 445 }
 446
 447 /*!Converts a list of global face numbers
 448  * to a distributed array with local face numbers
 449  */
 450 void ParallelTopology::convertGlobalFaceList(const int* face_list, int nbface, int* local, int* ip)
 451 {
 452   for (int i=0; i< nbface; i++)
 453   {
 454     INTERP_KERNEL::HashMap<int, pair<int,int> >::const_iterator iter = m_face_glob_to_loc.find(face_list[i]);
 455     if (iter == m_face_glob_to_loc.end())
 456     {
 457       throw MED_EXCEPTION("convertGlobalFaceList - Face  not found");
 458     }
 459     ip[i]=(iter->second).first;
 460     local[i]=(iter->second).second;
 461     //    cout << " in convertGlobalFAceList face global "<<face_list[i]<<" -> ("<<ip[i]<<","<<local[i]<<")"<<endl;
 462   }
 463 }
 464
 465 /*!Converts a list of global node numbers on domain ip
 466  * to a distributed array with local cell numbers.
 467  *
 468  * If a node in the list is represented on several domains,
 469  * only the value with domain ip is returned
 470  *
 471  * */
 472 void ParallelTopology::convertGlobalFaceList(const int* face_list, int nbface, int* local, int ip)
 473 {
 474   for (int i=0; i< nbface; i++)
 475   {
 476     typedef INTERP_KERNEL::HashMultiMap<int,pair<int,int> >::iterator mmiter;
 477     pair<mmiter,mmiter> range=m_face_glob_to_loc.equal_range(face_list[i]);
 478     for (mmiter it=range.first; it !=range.second; it++)
 479     {
 480       int ipfound=(it->second).first;
 481       if (ipfound==ip)
 482         local[i]=(it->second).second;
 483
 484     }
 485   }
 486 }
 487
 488
 489 ////creating node mapping
 490 void ParallelTopology::createNodeMapping(map<MED_EN::medGeometryElement,int*>& type_connectivity,
 491                                          map<MED_EN::medGeometryElement,int>& present_type_numbers,
 492                                          vector<int>& polygon_conn,
 493                                          vector<int>& polygon_conn_index,
 494                                          vector<int>& polyhedron_conn,
 495                                          vector<int>& polyhedron_conn_index,
 496                                          vector<int>& polyhedron_face_index,
 497                                          int idomain)
 498 {
 499   set<int> local_numbers;
 500   int local_index=0;
 501
 502   map<MED_EN::medGeometryElement,int>::const_iterator iter;
 503   for (iter = present_type_numbers.begin(); iter!=present_type_numbers.end();iter++)
 504   {
 505     int type=iter->first;
 506     int nodes_per_type= type%100;
 507
 508     if (!((type/100==m_mesh_dimension)
 509           ||(type==MED_EN::MED_POLYGON && m_mesh_dimension==2)
 510           ||(type==MED_EN::MED_POLYHEDRA && m_mesh_dimension==3))) continue;
 511
 512     if (type != MED_EN::MED_POLYGON && type != MED_EN::MED_POLYHEDRA)
 513     {
 514       for (int icell=0; icell<present_type_numbers[type]; icell++)
 515       {
 516         for (int inode=0; inode<nodes_per_type; inode++)
 517         {
 518           int global=type_connectivity[type][icell*nodes_per_type+inode];
 519           if(local_numbers.find(global)==local_numbers.end())
 520           {
 521             local_index++;
 522             local_numbers.insert(global);
 523             //m_node_loc_to_glob.insert(make_pair(make_pair(idomain,local_index),global));
 524             m_node_loc_to_glob[idomain].push_back(global);
 525             m_node_glob_to_loc.insert(make_pair(global,make_pair(idomain,local_index)));
 526             //          cout << "node : global ="<<global<<" local =("<<idomain<<","<<local_index<<")"<<endl;
 527           }
 528         }
 529
 530       }
 531     }
 532     else if (type== MED_EN::MED_POLYGON)
 533     {
 534       for ( unsigned i = 0; i < polygon_conn.size(); ++i )
 535       {
 536         int global=polygon_conn[i];
 537         if(local_numbers.find(global)==local_numbers.end())
 538         {
 539           local_index++;
 540           local_numbers.insert(global);
 541           m_node_loc_to_glob[idomain].push_back(global);
 542           m_node_glob_to_loc.insert(make_pair(global,make_pair(idomain,local_index)));
 543         }
 544       }
 545     }
 546     else if (type==MED_EN::MED_POLYHEDRA)
 547     {
 548       for ( unsigned i = 0; i < polyhedron_conn.size(); ++i )
 549       {
 550         int global=polyhedron_conn[i];
 551         if(local_numbers.find(global)==local_numbers.end())
 552         {
 553           local_index++;
 554           local_numbers.insert(global);
 555           m_node_loc_to_glob[idomain].push_back(global);
 556           m_node_glob_to_loc.insert(make_pair(global,make_pair(idomain,local_index)));
 557         }
 558       }
 559     }
 560
 561   }
 562   m_nb_nodes[idomain]=local_index;
 563 }
 564
 565 //================================================================================
 566 /*!
 567  * \brief Return true if the domain mesh contains a cell based on given global nodes
 568  */
 569 //================================================================================
 570
 571 bool ParallelTopology::hasCellWithNodes( const MESHCollection& new_collection,
 572                                          int                   domain,
 573                                          const set<int>&       globNodes)
 574 {
 575   // convert global nodes to local in the given domain
 576   set<int> nodes;
 577   set<int>::const_iterator n = globNodes.begin();
 578   for ( ; n != globNodes.end(); ++n )
 579     nodes.insert( convertGlobalNode( *n, domain ));
 580
 581   const MED_EN::medConnectivity connType = MED_EN::MED_NODAL;
 582   const MED_EN::medEntityMesh   entity   = MED_EN::MED_CELL;
 583
 584   // loop on all types of cells
 585   const MEDMEM::MESH* mesh = new_collection.getMesh( domain );
 586   int nbTypes = mesh->getNumberOfTypes( entity );
 587   const MED_EN::medGeometryElement * types = mesh->getTypes( entity );
 588   for ( int t = 0; t < nbTypes; ++t )
 589   {
 590     // get connectivity
 591     if ( !mesh->existConnectivity( connType, entity ))
 592       continue;
 593     int nbCell = mesh->getNumberOfElements( entity, types[t] );
 594     const int *conn, *index;
 595     conn  = mesh->getConnectivity(connType, entity, types[t]);
 596     index = mesh->getConnectivityIndex(connType, entity);
 597     // find a cell containing the first of given nodes,
 598     // then check if the found cell contains all the given nodes
 599     const int firstNode = *nodes.begin();
 600     for ( int i = 0; i < nbCell; ++i )
 601     {
 602       for ( int j = index[i]-1; j < index[i+1]-1; ++j )
 603         if ( conn[j] == firstNode )
 604         {
 605           unsigned nbSame = 0;
 606           for ( j = index[i]-1; j < index[i+1]-1; ++j )
 607             nbSame += nodes.count( conn[j] );
 608           if ( nbSame == nodes.size() )
 609             return true;
 610           break;
 611         }
 612     }
 613   }
 614   return false;
 615 }
 616
 617 ////creating face mapping
 618 void ParallelTopology::createFaceMapping(const MESHCollection& initial_collection,
 619                                          const MESHCollection& new_collection)
 620   //                     map<MED_EN::medGeometryElement,int*>& type_list,
 621   //                     map<MED_EN::medGeometryElement,int>& present_type_numbers,
 622   //                     int idomain
 623
 624 {
 625   // containers for the new topology
 626   vector<int> new_counts(m_nb_domain,0);
 627   vector<int> domain_counts(m_nb_domain,0);
 628   const Topology* old_topology=initial_collection.getTopology();
 629   int nb_domain_old=old_topology->nbDomain();
 630   int global_index=old_topology->getFaceNumber();
 631   //cout << "nb faces = " << global_index << endl;
 632   set <pair<int, pair<int,int> > > global_treated;
 633
 634   //definition of the d-1 constituent for the considered mesh dimension
 635   MED_EN::medEntityMesh constituent_entity;
 636   switch (m_mesh_dimension)
 637   {
 638   case 3:
 639     constituent_entity= MED_EN::MED_FACE;
 640     break;
 641   case 2:
 642     constituent_entity = MED_EN::MED_EDGE;
 643     break;
 644   }
 645
 646   for (int iold=0; iold<nb_domain_old;iold++)
 647   {
 648     if ( !initial_collection.getMesh(iold) ) continue;
 649     int nbtotalface = initial_collection.getMesh(iold)->getNumberOfElements(constituent_entity,MED_EN::MED_ALL_ELEMENTS);
 650     SCRUTE_MED(nbtotalface);
 651     const int* face_conn = 0;
 652     const int* face_offset = 0;
 653     if (nbtotalface >0)
 654     {
 655       face_conn = initial_collection.getMesh(iold)->getConnectivity(
 656                                                                     MED_EN::MED_NODAL,constituent_entity,MED_EN::MED_ALL_ELEMENTS);
 657       face_offset = initial_collection.getMesh(iold)->getConnectivityIndex(MED_EN::MED_NODAL,constituent_entity);
 658     }
 659     for (int iface=0;iface<nbtotalface; iface++)
 660     {
 661       int global_face_number = old_topology->convertFaceToGlobal(iold,iface+1);
 662
 663       //      int inode = face_offset[iface];
 664       for (int i=0; i<m_nb_domain; i++) domain_counts[i]=0;
 665       set <int> nodes;
 666       int nbnodes;
 667       {
 668         nbnodes=face_offset[iface+1]-face_offset[iface];
 669         for (int inode= face_offset[iface];inode < face_offset[iface+1]; inode++)
 670         {
 671           int node=face_conn[inode-1];
 672
 673           int global = old_topology->convertNodeToGlobal(iold,node);
 674           //        cout << "global node "<<global<<"ip "<<iold<< "noeud"<<node<<endl;
 675           nodes.insert(global);
 676           typedef INTERP_KERNEL::HashMultiMap<int,pair<int,int> >::iterator mmiter;
 677           pair<mmiter,mmiter> range=m_node_glob_to_loc.equal_range(global);
 678
 679           int ip;
 680           for (mmiter it=range.first; it !=range.second; it++)
 681           {
 682             ip=(it->second).first;
 683             domain_counts[ip]++;
 684           }
 685         }
 686       }
 687
 688       set<int>::const_iterator iter_node = nodes.begin();
 689       int numbers[3];
 690       numbers[2]=0; // for segments
 691       for (int i=0; i<nbnodes; i++)
 692       {
 693         numbers[i]=*iter_node;
 694         iter_node++;
 695       }
 696       set <pair<int, pair<int,int> > > ::iterator iter_triplets;
 697       pair<int, pair<int,int> > triplet = make_pair(numbers[0],make_pair(numbers[1],numbers[2]));
 698       iter_triplets=global_treated.find(triplet);
 699       if (iter_triplets==global_treated.end())
 700       {
 701         global_treated.insert(triplet);
 702         //  int nbnodes=face_offset[iface+1]-face_offset[iface];
 703         if (global_face_number == -1)
 704         {
 705           global_index++;
 706           global_face_number=global_index;
 707
 708         }
 709         //  SCRUTE_MED(nbnodes);
 710
 711         for (int inew=0;inew<m_nb_domain;inew++)
 712         {
 713           //     SCRUTE_MED(domain_counts[inew]);
 714           if(domain_counts[inew]==nbnodes)
 715           {
 716             if ( !hasCellWithNodes( new_collection, inew, nodes ))
 717               continue; // 0020861: EDF 1387 MED: Result of medsplitter gives standalone triangles
 718
 719             new_counts[inew]++;
 720             m_face_glob_to_loc.insert(make_pair(global_face_number,make_pair(inew,new_counts[inew])));
 721             //m_face_loc_to_glob.insert(make_pair(make_pair(inew,new_counts[inew]),global_face_number));
 722             m_face_loc_to_glob[inew].push_back(global_face_number);
 723           }
 724         }
 725         //cout << "global_face_number = " << global_face_number << endl;
 726       }
 727     }
 728   }
 729
 730   for (int inew=0;inew<m_nb_domain;inew++)
 731   {
 732     m_nb_faces[inew]=new_counts[inew];
 733     MESSAGE_MED(" Nb faces ["<<inew<<"]="<<m_nb_faces[inew]);
 734   }
 735   MESSAGE_MED(" total number of faces"<<getFaceNumber());
 736 }
 737
 738 ////creating node mapping
 739 void ParallelTopology::createFaceMapping2ndversion(const MESHCollection& initial_collection)
 740
 741 {
 742   // containers for the new topology
 743   vector<int> new_counts(m_nb_domain,0);
 744   vector<int> domain_counts(m_nb_domain,0);
 745   const Topology* old_topology=initial_collection.getTopology();
 746   int nb_domain_old=old_topology->nbDomain();
 747   //int global_index=old_topology->getFaceNumber();
 748   //  set <pair<int, pair<int,int> > > global_treated;
 749
 750   //definition of the d-1 constituent for the considered mesh dimension
 751   MED_EN::medEntityMesh constituent_entity;
 752   switch (m_mesh_dimension)
 753   {
 754   case 3:
 755     constituent_entity= MED_EN::MED_FACE;
 756     break;
 757   case 2:
 758     constituent_entity = MED_EN::MED_EDGE;
 759     break;
 760   }
 761
 762   for (int iold=0; iold<nb_domain_old;iold++)
 763   {
 764     int nbcell=old_topology->getCellNumber(iold);
 765
 766     const int* face_conn = initial_collection.getMesh(iold)->
 767       getConnectivity(MED_EN::MED_DESCENDING,MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
 768     const int* face_offset = initial_collection.getMesh(iold)->
 769       getConnectivityIndex(MED_EN::MED_DESCENDING,MED_EN::MED_CELL);
 770     MESSAGE_MED("end of connectivity calculation");
 771     set<int> global_treated;
 772     for (int icell=0; icell<nbcell; icell++)
 773     {
 774       int global_cell_number=old_topology->convertCellToGlobal(iold,icell+1);
 775       int inew = getCellDomainNumber(global_cell_number);
 776
 777       for (int iface = face_offset[icell]; iface < face_offset[icell+1]; iface++)
 778       {
 779         int global_face_number=old_topology->convertFaceToGlobal(iold,abs(face_conn[iface-1]));
 780         if (global_treated.find (global_face_number)==global_treated.end())
 781         {
 782           new_counts[inew]++;
 783           m_face_glob_to_loc.insert(make_pair(global_face_number,make_pair(inew,new_counts[inew])));
 784
 785           //m_face_loc_to_glob.insert(make_pair(make_pair(inew,new_counts[inew]),global_face_number));
 786           m_face_loc_to_glob[inew].push_back(global_face_number);
 787           global_treated.insert(global_face_number);
 788
 789         }
 790       }
 791     }
 792   }
 793
 794
 795   for (int inew=0;inew<m_nb_domain;inew++)
 796   {
 797     m_nb_faces[inew]=new_counts[inew];
 798     //  cout << " Nb faces ["<<inew<<"]="<<m_nb_faces[inew]<<endl;
 799   }
 800   MESSAGE_MED(" total number of faces"<<getFaceNumber());
 801 }
 802
 803
 804 //replacing a table of global numbering with a table with local numberings
 805 // type_connectivity contains global connectivity for each type in input
 806 // type_connectivity contains local connectivity for each type in output
 807 void ParallelTopology::convertToLocal(map<MED_EN::medGeometryElement,int*>& type_connectivity,
 808                                       map<MED_EN::medGeometryElement,int>& present_type_numbers,
 809                                       int idomain,
 810                                       MED_EN::medEntityMesh entity)
 811 {
 812   int dimension;
 813   switch (entity)
 814   {
 815   case MED_EN::MED_CELL:
 816     dimension=m_mesh_dimension;
 817     break;
 818   case MED_EN::MED_FACE:
 819     dimension=2;
 820     break;
 821   case MED_EN::MED_EDGE:
 822     dimension=1;
 823     break;
 824   }
 825
 826   MED_EN::MESH_ENTITIES::const_iterator currentEntity;
 827   list<MED_EN::medGeometryElement>::const_iterator iter;
 828   currentEntity  = MED_EN::meshEntities.find(MED_EN::MED_CELL);
 829
 830   for (iter = (*currentEntity).second.begin();iter != (*currentEntity).second.end(); iter++)
 831   {
 832     MED_EN::medGeometryElement type = (*iter);
 833     if (type/100 != dimension) continue;
 834     for (int inode=0; inode<present_type_numbers[type]*(type%100); inode++)
 835     {
 836       //      cout <<" inode :"<<inode<< " global = "<<type_connectivity[type][inode];
 837       int global = type_connectivity[type][inode];
 838       typedef INTERP_KERNEL::HashMultiMap<int,pair<int,int> >::iterator mmiter;
 839       pair<mmiter,mmiter> range=m_node_glob_to_loc.equal_range(global);
 840       for (mmiter it=range.first; it !=range.second; it++)
 841       {
 842         if ((it->second).first==idomain)
 843           type_connectivity[type][inode]=(it->second).second;
 844       }
 845       //      cout << " local = " <<type_connectivity[type][inode]<<endl;
 846     }
 847
 848   }
 849 }
 850
 851 //replacing a table of global numbering with a table with local numberings
 852 // type_connectivity contains global connectivity for each type in input
 853 // type_connectivity contains local connectivity for each type in output
 854 void ParallelTopology::convertToLocal2ndVersion(int* nodes, int nbnodes, int idomain)
 855 {
 856   for (int inode=0; inode<nbnodes; inode++)
 857   {
 858     //      cout <<" inode :"<<inode<< " global = "<<type_connectivity[type][inode];
 859     int global = nodes[inode];
 860     typedef INTERP_KERNEL::HashMultiMap<int,pair<int,int> >::iterator mmiter;
 861     pair<mmiter,mmiter> range=m_node_glob_to_loc.equal_range(global);
 862     for (mmiter it=range.first; it !=range.second; it++)
 863     {
 864       if ((it->second).first==idomain)
 865         nodes[inode]=(it->second).second;
 866     }
 867   }
 868 }
 869
 870
 871 //! computing arrays with node/node correspondencies
 872 void ParallelTopology::computeNodeNodeCorrespondencies(int idomain, vector<MEDMEM::MEDSKYLINEARRAY*>& corr) const
 873 {
 874   vector <int> sizes (m_nb_domain,0);
 875   vector <int*> values (m_nb_domain);
 876   for (int i=0; i<m_nb_domain; i++)
 877   {
 878     values[i]=0;
 879   }
 880
 881   for (int inode=0; inode<m_nb_nodes[idomain]; inode++)
 882   {
 883     //int global = (m_node_loc_to_glob.find(make_pair(idomain,inode+1)))->second;
 884     int global = m_node_loc_to_glob[idomain][inode];
 885     typedef INTERP_KERNEL::HashMultiMap<int,pair<int,int> >::const_iterator mm;
 886     pair<mm,mm> range = m_node_glob_to_loc.equal_range(global);
 887     for (mm it=range.first; it !=range.second; it++)
 888     {
 889       int id = (it->second).first;
 890       if (id !=idomain)
 891       {
 892         sizes[id]++;
 893       }
 894     }
 895   }
 896
 897   for (int ip=0; ip<m_nb_domain; ip++)
 898   {
 899     if (sizes[ip]>0)
 900       values[ip]=new int[2*sizes[ip]];
 901     sizes[ip]=0;
 902   }
 903
 904
 905   for (int inode=0; inode<m_nb_nodes[idomain]; inode++)
 906   {
 907     //int global = (m_node_loc_to_glob.find(make_pair(idomain,inode+1)))->second;
 908     int global = m_node_loc_to_glob[idomain][inode];
 909     typedef INTERP_KERNEL::HashMultiMap<int,pair<int,int> >::const_iterator mm;
 910     pair<mm,mm> range = m_node_glob_to_loc.equal_range(global);
 911     for (mm it=range.first; it !=range.second; it++)
 912     {
 913       int id = (it->second).first;
 914       if (id !=idomain)
 915       {
 916         values[id][sizes[id]*2]=inode+1;
 917         values[id][sizes[id]*2+1]=(it->second).second;
 918         sizes[id]++;
 919       }
 920     }
 921   }
 922
 923   for (int i=0; i< m_nb_domain; i++)
 924   {
 925     if (sizes[i]!=0)
 926     {
 927       int* index = new int[sizes[i]+1];
 928       for (int j=0; j<sizes[i]+1; j++)
 929         index[j]=j+1;
 930       corr[i]=new MEDMEM::MEDSKYLINEARRAY(sizes[i],2*sizes[i],index,values[i]);
 931       delete[] index;
 932       delete[] values[i];
 933     }
 934   }
 935 }
 936
 937 //================================================================================
 938 /*!
 939  * \brief computing arrays with cell/cell correspondencies
 940  *  \param idomain - domain for which to compute correspondencies
 941  *  \param corr - out correspondencies
 942  *  \param graph - graph containing new cell distribution among domains
 943  *  \param id_shift - shitf to get cell id global on this processor from id global
 944  *                    over all procs
 945  */
 946 //================================================================================
 947
 948 void ParallelTopology::computeCellCellCorrespondencies(int idomain, vector<MEDMEM::MEDSKYLINEARRAY*>& corr, const Graph* graph) const
 949 {
 950   vector <int> sizes (m_nb_domain,0);
 951   vector <int*> values (m_nb_domain);
 952   for (int i=0; i<m_nb_domain; i++)
 953   {
 954     values[i]=0;
 955   }
 956
 957   vector <INTERP_KERNEL::HashMultiMap<int,int> > cell_corresp;
 958   //TODO : remplacer int* par une map <int,int>
 959   //  vector<int*  > number_of_connections(m_nb_domain);
 960   //  vector<map<int,int> > number_of_connections;
 961   vector<map<int,int> > number_of_connections;
 962   cell_corresp.resize(m_nb_domain);
 963   number_of_connections.resize(m_nb_domain);
 964   //  for (int i=0; i<m_nb_domain; i++)
 965   //    {
 966   //      //    cell_corresp[i]=new multimap<int,int>;
 967   //      if (m_nb_cells[i] >0)
 968   //        {
 969   //          number_of_connections[i]=new int[m_nb_cells[idomain]];
 970   //          for (int j=0; j<m_nb_cells[idomain]; j++)
 971   //            number_of_connections[i][j]=0;
 972   //        }
 973   //    }
 974
 975   const MEDMEM::MEDSKYLINEARRAY* skylinegraph = graph->getGraph();
 976
 977   const int* index=skylinegraph->getIndex();
 978   const int* value=skylinegraph->getValue();
 979
 980   TGlob2DomainLoc::const_iterator gl_do_lo_end = m_glob_to_loc.end();
 981
 982   for (int icell=0; icell<m_nb_cells[idomain]; icell++)
 983   {
 984     int global= m_loc_to_glob[idomain][icell];
 985     for (int ii=index[global-1]-1; ii<index[global]-1; ii++)
 986     {
 987       int distant_global=value[ii];
 988
 989       const pair<int,int>& local = m_glob_to_loc.find(distant_global)->second;
 990
 991       if (local.first != idomain)
 992       {
 993         cell_corresp[local.first].insert(make_pair(icell+1,local.second));
 994         //              number_of_connections[local.first][icell]++;
 995         if (number_of_connections[local.first].find(icell)==number_of_connections[local.first].end())
 996           number_of_connections[local.first].insert(make_pair(icell,1));
 997         else
 998           number_of_connections[local.first][icell]++;
 999
1000       }
1001     }
1002   }
1003
1004   for (int inew=0; inew<m_nb_domain; inew++)
1005   {
1006     if (inew==idomain || number_of_connections[inew].empty()) continue;
1007
1008     int* new_index=new int[m_nb_cells[idomain]+1];
1009     new_index[0]=1;
1010     for (int i=0; i<m_nb_cells[idomain]; i++)
1011     {
1012
1013       if (number_of_connections[inew].find(i)!=number_of_connections[inew].end())
1014         new_index[i+1]=new_index[i]+number_of_connections[inew][i];
1015       else
1016         new_index[i+1]=new_index[i];
1017     }
1018     int* new_value;
1019     if (new_index[m_nb_cells[idomain]]-1 > 0)
1020       new_value=new int[new_index[m_nb_cells[idomain]]-1];
1021     else
1022       new_value=0;
1023
1024     int value_i=0;
1025
1026     //                      INTERP_KERNEL::HashMultiMap<int,int>::iterator iter=cell_corresp[inew].begin();
1027
1028     for (int i=0; i<m_nb_cells[idomain]; i++)
1029     {
1030       //          for (int j=new_index[i];j<new_index[i+1];j++,value_i++,iter++)
1031       //            new_value[value_i]=iter->second;
1032
1033       typedef INTERP_KERNEL::HashMultiMap<int,int>::iterator mmiter;
1034       pair<mmiter,mmiter> range=cell_corresp[inew].equal_range(i+1);
1035       for (mmiter it=range.first; it!=range.second; it++)
1036       {
1037         new_value[value_i]=it->second;
1038         value_i++;
1039       }
1040     }
1041     if (value_i>0)
1042     {
1043       corr[inew] = new MEDMEM::MEDSKYLINEARRAY(m_nb_cells[idomain], new_index[m_nb_cells[idomain]]-1, new_index,new_value,true);
1044     }
1045     else
1046     {
1047       corr[inew]=0;
1048       if (new_value !=0) delete[]new_value;
1049       delete[]new_index;
1050     }
1051
1052   }
1053
1054   //    for (int inew=0; inew<m_nb_domain; inew++)
1055   //    if (m_nb_cells[inew]>0)
1056   //      delete[] number_of_connections[inew];
1057
1058 }
1059
1060 //================================================================================
1061 /*!
1062  * \brief Return max global face number
1063  */
1064 //================================================================================
1065
1066 int ParallelTopology::getMaxGlobalFace() const
1067 {
1068   int max = 0;
1069   TGlob2LocsMap::const_iterator g_l_l = m_face_glob_to_loc.begin();
1070   for ( ; g_l_l != m_face_glob_to_loc.end(); ++g_l_l )
1071     if ( g_l_l->first > max )
1072       max = g_l_l->first;
1073   return max;
1074 }
1075
1076 void ParallelTopology::recreateFaceMapping(const TGeom2FacesByDomian& face_map)
1077 {
1078   m_face_glob_to_loc.clear();
1079   for (int i=0; i<m_nb_domain;i++)
1080     m_face_loc_to_glob[i].clear();
1081
1082   for (int idomain=0; idomain<m_nb_domain; idomain++)
1083   {
1084     int ilocal=1;
1085     TGeom2Faces::const_iterator iter = face_map[idomain].begin();
1086     for (; iter != face_map[idomain].end(); iter++)
1087     {
1088       for (unsigned i=0; i< (iter->second).size(); i++)
1089       {
1090         MEDSPLITTER_FaceModel* face = (iter->second)[i];
1091         //cout << "global :"<<face->getGlobal()<<" - "<<ilocal<<endl;
1092         m_face_glob_to_loc.insert(make_pair(face->getGlobal(),make_pair(idomain,ilocal)));
1093         m_face_loc_to_glob[idomain].push_back(face->getGlobal());
1094         ilocal++;
1095       }
1096     }
1097     m_nb_faces[idomain] =ilocal-1;
1098   }
1099 }
1100
1101 //================================================================================
1102 /*!
1103  * \brief Recreating cell and node mapping after send-reveive and fusion of domain meshes
1104  */
1105 //================================================================================
1106
1107 void ParallelTopology::recreateMappingAfterFusion(const vector<MEDMEM::MESH*>& meshes)
1108 {
1109   const char* LOC = "ParallelTopology::recreateMappingAfterFusion(): ";
1110
1111   m_glob_to_loc.clear();
1112   m_node_glob_to_loc.clear();
1113   m_face_glob_to_loc.clear();
1114
1115   for (int idomain=0; idomain<m_nb_domain; idomain++)
1116   {
1117     m_nb_cells[idomain] = m_nb_nodes[idomain] = m_nb_faces[idomain] = 0;
1118     m_loc_to_glob[idomain].clear();
1119     m_node_loc_to_glob[idomain].clear();
1120     m_face_loc_to_glob[idomain].clear();
1121
1122     if ( !meshes[idomain]->getCoordinateptr() ) continue; // empty domian
1123
1124     //creating cell maps
1125
1126     m_nb_cells[idomain]=meshes[idomain]->getNumberOfElements(MED_EN::MED_CELL,
1127                                                              MED_EN::MED_ALL_ELEMENTS);
1128     if ( m_cell_loc_to_glob_fuse[idomain].size() != m_nb_cells[idomain] )
1129       throw MED_EXCEPTION(MEDMEM::STRING(LOC)<<" invalid nb of fused cells");
1130
1131     m_loc_to_glob[idomain].swap( m_cell_loc_to_glob_fuse[idomain] );
1132
1133     for (int i=0; i<m_nb_cells[idomain]; i++)
1134     {
1135       int global=m_loc_to_glob[idomain][i];
1136       m_glob_to_loc[global]=make_pair(idomain,i+1);
1137     }
1138
1139     //creating node maps
1140
1141     m_nb_nodes[idomain]=meshes[idomain]->getNumberOfNodes();
1142     m_node_loc_to_glob[idomain] = ((MEDMEM::MeshFuse*)meshes[idomain])->getNodeNumbers();
1143     if ( m_node_loc_to_glob[idomain].size() != m_nb_nodes[idomain] )
1144       throw MED_EXCEPTION(MEDMEM::STRING(LOC)<<" invalid nb of fused nodes");
1145
1146     // setting mappings for all nodes
1147     for (int inode=0; inode<m_nb_nodes[idomain]; inode++)
1148     {
1149       int global_number=m_node_loc_to_glob[idomain][inode];
1150       m_node_glob_to_loc.insert(make_pair(global_number,make_pair(idomain,inode+1)));
1151     }
1152
1153
1154     //creating dimension d-1 component mappings
1155
1156     MED_EN::medEntityMesh constituent_entity =
1157       (m_mesh_dimension == 3 ? MED_EN::MED_FACE : MED_EN::MED_EDGE );
1158     m_nb_faces[idomain] = meshes[idomain]->getNumberOfElements(constituent_entity,
1159                                                                MED_EN::MED_ALL_ELEMENTS);
1160     if ( m_face_loc_to_glob_fuse[idomain].size() != m_nb_faces[idomain] )
1161       throw MED_EXCEPTION(MEDMEM::STRING(LOC)<<" invalid nb of fused faces of domain "<< idomain
1162                           << ": expect " << m_nb_faces[idomain]
1163                           << " but have " << m_face_loc_to_glob_fuse[idomain].size());
1164
1165     m_face_loc_to_glob[idomain].swap( m_face_loc_to_glob_fuse[idomain] );
1166
1167     for (int iface=0; iface<m_nb_faces[idomain]; iface++)
1168     {
1169       int global_number=m_face_loc_to_glob[idomain][iface];
1170       m_face_glob_to_loc.insert(make_pair(global_number,make_pair(idomain,iface+1)));
1171     }
1172   }
1173
1174 }