src/MEDMEM/MEDMEM_MeshFuse.cxx

   1 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
   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 // File      : MEDMEM_MeshFuse.cxx
  20 // Created   : Tue Jul  7 18:27:00 2009
  21 // Author    : Edward AGAPOV (eap)
  22
  23 #include "MEDMEM_MeshFuse.hxx"
  24
  25 #include "MEDMEM_Family.hxx"
  26 #include "MEDMEM_Group.hxx"
  27
  28 using namespace MEDMEM;
  29 using namespace MED_EN;
  30 using namespace std;
  31
  32 #define _NODE_TYPE_ MED_NONE
  33
  34 MeshFuse::MeshFuse():MESHING()
  35 {
  36 }
  37
  38 MeshFuse::~MeshFuse()
  39 {
  40 }
  41
  42 //================================================================================
  43 /*!
  44  * \brief set global numbers of nodes if MeshFuse has been filled via MESHING
  45  */
  46 //================================================================================
  47
  48 void MeshFuse::setNodeNumbers( const std::vector<int>& node_glob_numbers )
  49 {
  50   const char* LOC = "MeshFuse::setNodeNumbers(node_glob_numbers): ";
  51
  52   if ( !_node_glob_numbers.empty() )
  53     throw MEDEXCEPTION(STRING(LOC)<<"node numbers has been already set");
  54
  55   if ( node_glob_numbers.size() != getNumberOfNodes() &&
  56        node_glob_numbers.size() > 0 && getNumberOfNodes() > 0 )
  57     throw MEDEXCEPTION
  58       (STRING(LOC)<<"size of node_glob_numbers must be equal number of nodes in MeshFuse");
  59
  60   _node_glob_numbers = node_glob_numbers;
  61 }
  62
  63 //================================================================================
  64 /*!
  65  * \brief Fuse me and other mesh
  66  *  \param mesh - mesh to concatenate
  67  *  \param node_glob_numbers - node numbers used to merge coincident meshes entities
  68  */
  69 //================================================================================
  70
  71 void MeshFuse::concatenate( const MESH* mesh, const vector<int>& node_glob_numbers )
  72 {
  73   const char* LOC = "MeshFuse::concatenate( mesh, node_glob_numbers ): ";
  74   if ( !mesh || mesh->getNumberOfNodes() == 0) return;
  75
  76   _mesh = mesh;
  77
  78   if ( this->getNumberOfNodes() < 1 )
  79   {
  80     mesh->getCoordinates( MED_FULL_INTERLACE );// make fields filled in case of GRID
  81     mesh->getConnectivityptr();
  82
  83     // just copy mesh
  84     static_cast<MESH&>(*this) = *mesh;
  85
  86     _node_glob_numbers = node_glob_numbers;
  87     return;
  88   }
  89   // check feasibility
  90
  91   if ( mesh->getNumberOfNodes() > 0 && node_glob_numbers.empty() )
  92     throw MEDEXCEPTION(STRING(LOC)<<"merging without node global numbers not implemented yet");
  93
  94   if ( mesh->getNumberOfNodes() != node_glob_numbers.size() )
  95     throw MEDEXCEPTION(STRING(LOC)<<"invalid number of node global numbers");
  96
  97   if ( mesh->getSpaceDimension() != this->getSpaceDimension() ||
  98        mesh->getMeshDimension()  != this->getMeshDimension() )
  99     throw MEDEXCEPTION(STRING(LOC)<<"can't unite meshes with different dimensions so far, sorry");
 100
 101   // init
 102   _merged_of_type.clear();
 103   for ( int i = 0; i < NB_INDICES; ++i )
 104   {
 105     _nb_index[ i ].clear();
 106     _nb_index[ i ].resize( MED_ALL_ENTITIES );
 107   }
 108
 109   // concatenation
 110
 111   int final_nb_nodes = makeNewNodeIds( node_glob_numbers );
 112
 113   expandCoordinates(final_nb_nodes);
 114
 115   expandConnectivity(final_nb_nodes);
 116
 117   expandSupports();
 118
 119   // clear unnecessary data
 120   _new_elem_ids_of_type.clear();
 121 }
 122
 123 //================================================================================
 124 /*!
 125  * \brief Return number of nodes in the expanded mesh
 126  *  \param add_glob_numbers - global ids of nodes to append
 127  */
 128 //================================================================================
 129
 130 int MeshFuse::makeNewNodeIds(const vector<int>& add_glob_numbers)
 131 {
 132   // remember merged added nodes and make an array of new ids of added nodes
 133   vector<int>&       merged = _merged_of_type      [ _NODE_TYPE_ ];
 134   vector<int>& new_node_ids = _new_elem_ids_of_type[ _NODE_TYPE_ ];
 135   new_node_ids.resize( add_glob_numbers.size() );
 136
 137   // extand global node numbers
 138   _node_glob_numbers.reserve( _node_glob_numbers.size() + add_glob_numbers.size());
 139
 140   map<int,int> glob_ids; // to map global ids to new local ones
 141   for ( int i = 0; i < _node_glob_numbers.size(); ++i )
 142     glob_ids.insert( make_pair( _node_glob_numbers[i], i+1 ));
 143
 144   int next_loc_id = getNumberOfNodes() + 1;
 145
 146   for ( int i = 0; i < add_glob_numbers.size(); ++i )
 147   {
 148     map<int,int>::iterator  glob_loc =
 149       glob_ids.insert( make_pair( add_glob_numbers[i], next_loc_id )).first;
 150
 151     new_node_ids[i] = glob_loc->second;
 152
 153     if ( new_node_ids[i] == next_loc_id ) // unique global number
 154     {
 155       next_loc_id++;
 156       _node_glob_numbers.push_back( add_glob_numbers[i] );
 157     }
 158     else
 159     {
 160       merged.push_back( i );
 161     }
 162   }
 163   _nb_index[ INIT_OLD ][ MED_NODE ][ _NODE_TYPE_ ] = getNumberOfNodes();
 164   _nb_index[ INIT_ADD ][ MED_NODE ][ _NODE_TYPE_ ] = add_glob_numbers.size();
 165   _nb_index[ RSLT_ADD ][ MED_NODE ][ _NODE_TYPE_ ] = add_glob_numbers.size() - merged.size();
 166
 167   return next_loc_id - 1;
 168 }
 169
 170 //================================================================================
 171 /*!
 172  * \brief Update coordinates
 173  */
 174 //================================================================================
 175
 176 void MeshFuse::expandCoordinates(int final_nb_nodes)
 177 {
 178   const int dim = getSpaceDimension();
 179
 180   // create new coordinates
 181   double* coord = new double[ final_nb_nodes * dim ];
 182   MEDARRAY<double> medarray( coord, dim, final_nb_nodes, MED_FULL_INTERLACE,
 183                              /*shallowCopy=*/true,/*ownershipOfValues=*/false);
 184   // copy old coordinates
 185   int nb_old_coord = getNumberOfNodes() * dim;
 186   memcpy( coord, getCoordinates( MED_FULL_INTERLACE ), nb_old_coord * sizeof( double ));
 187
 188   // set new coords
 189   coord += nb_old_coord;
 190   const double* add_coord =_mesh->getCoordinates( MED_FULL_INTERLACE );
 191   if ( _merged_of_type[ _NODE_TYPE_ ].empty())
 192   {
 193     // no coincident nodes in the two meshes, just add coords
 194     memcpy( coord, add_coord, _mesh->getNumberOfNodes() * dim * sizeof( double ));
 195   }
 196   else
 197   {
 198     // copy coord of only unique nodes
 199     int first_added_node = getNumberOfNodes() + 1;
 200     const vector<int>& new_node_ids = _new_elem_ids_of_type[ _NODE_TYPE_ ];
 201     for ( int n = 0; n < new_node_ids.size(); ++n )
 202     {
 203       if ( new_node_ids[n] < first_added_node ) continue; // coincident node
 204       memcpy( coord, add_coord + n * dim, dim * sizeof( double ));
 205       coord += dim;
 206     }
 207   }
 208   _coordinate->setCoordinates( &medarray, /*shallowCopy=*/true );
 209
 210   _numberOfNodes = final_nb_nodes;
 211 }
 212
 213 //================================================================================
 214 /*!
 215  * \brief Concatenate connectivity of meshes
 216  *
 217  * Current implementation impies that cells can't coincide in the meshes
 218  */
 219 //================================================================================
 220
 221 void MeshFuse::expandConnectivity(int final_nb_nodes)
 222 {
 223   const medConnectivity nodal   = MED_NODAL;
 224
 225   // fill in _nb_index[ INIT_OLD ]
 226   for ( medEntityMesh entity = MED_CELL; entity < MED_ALL_ENTITIES; ++entity )
 227   {
 228     if ( existConnectivity( nodal, entity ))
 229     {
 230       const int *               index = this->getGlobalNumberingIndex(entity);
 231       const medGeometryElement* types = this->getTypes(entity);
 232       int                    nb_types = this->getNumberOfTypes(entity);
 233       for ( int t = 0; t < nb_types; ++t )
 234         _nb_index[ INIT_OLD ][ entity ][ types[t] ] = index[t+1]-index[0];
 235     }
 236   }
 237
 238   CONNECTIVITY *prev_connectivity = 0, *cell_connectivity = 0;
 239
 240   // loop on all entities
 241   for ( medEntityMesh entity = MED_CELL; entity < MED_ALL_ENTITIES; ++entity )
 242   {
 243     if ( entity == MED_FACE && getMeshDimension() == 2 )
 244       continue; // there can be EDGE
 245
 246     if ( !_mesh->existConnectivity( nodal, entity ))
 247     {
 248       // no entity in added mesh
 249       if ( existConnectivity( nodal, entity ) && prev_connectivity )
 250         prev_connectivity->setConstituent
 251           ( new CONNECTIVITY( *getConnectivityptr()->getConstituent( entity )));
 252       break;
 253     }
 254     if ( !existConnectivity( nodal, entity ))
 255     {
 256       // no entity in the old mesh - fully copy the added connectivity
 257       CONNECTIVITY* connectivity =
 258         new CONNECTIVITY( *_mesh->getConnectivityptr()->getConstituent( entity ));
 259       connectivity->setNumberOfNodes( final_nb_nodes );
 260       updateNodeIds( connectivity );
 261       if ( entity == MED_CELL )
 262         cell_connectivity = connectivity;
 263       else
 264         cell_connectivity->setConstituent( connectivity );
 265
 266       break;
 267     }
 268     // Fill in _nb_index[ INIT_ADD ]
 269     const int *                   index = _mesh->getGlobalNumberingIndex(entity);
 270     const medGeometryElement* add_types = _mesh->getTypes(entity);
 271     int                    nb_add_types = _mesh->getNumberOfTypes(entity);
 272     for ( int t = 0; t < nb_add_types; ++t )
 273       _nb_index[ INIT_ADD ][ entity ][ add_types[t] ] = index[t+1]-index[0];
 274
 275     // Unite connectivities of std types
 276
 277     // count types
 278     set<medGeometryElement> types;
 279     types.insert( this->getTypes(entity), this->getTypes(entity) + this->getNumberOfTypes(entity));
 280     types.insert(_mesh->getTypes(entity),_mesh->getTypes(entity) +_mesh->getNumberOfTypes(entity));
 281
 282     int sum_old = 0, sum_add = 0;
 283
 284     // make all data
 285     vector< TConnData >          conn_of_type( types.size() );
 286     vector< medGeometryElement > type_vec( types.size() );
 287     vector< int  >               count(1,1);
 288     set<medGeometryElement>::iterator type = types.begin();
 289     for ( int t = 0; type != types.end(); ++type, ++t )
 290     {
 291       int max_conn_len = 0;
 292       if ( this->getNumberOfElements( entity, *type ))
 293         max_conn_len += this->getConnectivityLength( nodal, entity, *type);
 294       if ( _mesh->getNumberOfElements( entity, *type ))
 295         max_conn_len += _mesh->getConnectivityLength( nodal, entity, *type);
 296       conn_of_type[t]._connectivity.reserve( max_conn_len );
 297
 298       int nb_old = appendConnectivity( conn_of_type[t], this, entity, *type );
 299       sum_old += nb_old;
 300       _nb_index[ INIT_OLD ][ entity ][ *type ] = sum_old;
 301
 302       int nb_add = appendConnectivity( conn_of_type[t],_mesh, entity, *type );
 303       sum_add += nb_add;
 304       _nb_index[ RSLT_ADD ][ entity ][ *type ] = sum_add;
 305
 306       count.push_back( count.back() + conn_of_type[t]._nb_elems );
 307       type_vec[t] = *type;
 308     }
 309     // make new connectivity
 310     CONNECTIVITY* connectivity = new CONNECTIVITY( types.size(), entity );
 311     connectivity->setNumberOfNodes ( final_nb_nodes );
 312     connectivity->setGeometricTypes( &type_vec[0], entity );
 313     connectivity->setCount         ( &count   [0], entity );
 314     for ( int t = 0; t < types.size(); ++t )
 315       connectivity->setNodal( &conn_of_type[t]._connectivity[0],
 316                               entity, type_vec[t],
 317                               &conn_of_type[t]._index[0] );
 318
 319     // store connectivity of an entity
 320     if ( !prev_connectivity )
 321     {
 322       cell_connectivity = connectivity;
 323       prev_connectivity = cell_connectivity;
 324     }
 325     else
 326     {
 327       prev_connectivity->setConstituent( connectivity );
 328       prev_connectivity = connectivity;
 329     }
 330
 331   } // loop on entities
 332
 333   if ( cell_connectivity )
 334   {
 335     delete _connectivity;
 336     _connectivity = cell_connectivity;
 337   }
 338 }
 339
 340 //================================================================================
 341 /*!
 342  * \brief Update node ids in the copied connectivity of theadded mesh
 343  */
 344 //================================================================================
 345
 346 void MeshFuse::updateNodeIds( CONNECTIVITY* connectivity )
 347 {
 348   const medConnectivity   nodal = MED_NODAL;
 349   const medGeometryElement type = MED_ALL_ELEMENTS;
 350
 351   const vector<int>& new_node_ids = _new_elem_ids_of_type[ _NODE_TYPE_ ];
 352
 353   medEntityMesh entity = connectivity->getEntity();
 354
 355   for ( ; entity < MED_ALL_ENTITIES; entity++ )
 356   {
 357     // Collect all connectivities with their lengths
 358
 359     list< pair< const int*, int > > conn_len_list;
 360
 361     if ( connectivity->existConnectivity( nodal, entity ))
 362       conn_len_list.push_back
 363         ( make_pair( connectivity->getConnectivity(nodal,entity,type),
 364                      connectivity->getConnectivityLength(nodal,entity,type)));
 365
 366     // Convert them
 367
 368     list< pair< const int*, int > >::iterator conn_len = conn_len_list.begin();
 369     for ( ; conn_len != conn_len_list.end(); ++conn_len )
 370     {
 371       int* conn = (int* ) conn_len->first;
 372       int* conn_end = conn + conn_len->second;
 373       for ( ; conn < conn_end; ++conn )
 374         *conn = new_node_ids[ *conn-1 ];
 375     }
 376   }
 377 }
 378
 379 //================================================================================
 380 /*!
 381  * \brief Concatenate connectivities
 382  *  \param data - storage of resulting connectivities
 383  *  \param mesh - mesh to get connectivity from
 384  *  \param entity - mesh entity to process
 385  *  \param type - geom type to process
 386  *  \retval int - nb of appended elements
 387  */
 388 //================================================================================
 389
 390 int MeshFuse::appendConnectivity( MeshFuse::TConnData& data,
 391                                   const MESH*          mesh,
 392                                   medEntityMesh        entity,
 393                                   medGeometryElement   type)
 394 {
 395   // get connectivity of type
 396
 397   const int* conn, *index = 0;
 398   int nb_elem, conn_len;
 399
 400   nb_elem = mesh->getNumberOfElements ( entity, type );
 401   if ( !nb_elem ) return 0;
 402   conn    = mesh->getConnectivity     ( MED_NODAL, entity, type );
 403   index   = mesh->getConnectivityIndex( MED_NODAL, entity );
 404   int shift = getElemNbShift( entity, type, (mesh == this ? INIT_OLD : INIT_ADD), /*prev=*/true);
 405   index += shift;
 406   conn_len = index[ nb_elem ] - index[ 0 ];
 407
 408   bool need_index = ( type == MED_POLYGON || type == MED_POLYHEDRA );
 409   if ( !need_index )
 410     data._index.resize( 1, 0 ); // for safe access to pointer even if no real index exists
 411
 412   // Append
 413
 414   int old_size = data._nb_elems;
 415   data._nb_elems += nb_elem;
 416
 417   if ( mesh == this )
 418   {
 419     // append connectivity to data as is
 420     data._connectivity.insert( data._connectivity.end(), conn, conn + conn_len );
 421     if ( need_index )
 422     {
 423       if ( data._index.empty() )
 424         data._index.insert( data._index.end(), index, index + nb_elem + 1 );
 425       else
 426         for ( int i = 0; i < nb_elem; ++i )
 427           data._index.push_back( data._index.back() + index[i+1] - index[i] );
 428     }
 429   }
 430   else
 431   {
 432     // convert connectivity of other mesh
 433
 434     vector<int> other_conn( conn_len );
 435     const vector<int>& new_node_ids = _new_elem_ids_of_type[ _NODE_TYPE_ ];
 436     if ( type == MED_POLYHEDRA )
 437       {
 438         for ( int n = 0; n < conn_len; ++n )
 439           if ( conn[ n ] > 0 )
 440             other_conn[ n ] = new_node_ids[ conn[ n ]-1 ];
 441       }
 442     else
 443       {
 444         for ( int n = 0; n < conn_len; ++n )
 445           other_conn[ n ] = new_node_ids[ conn[ n ]-1 ];
 446       }
 447     if ( entity == MED_CELL || _merged_of_type[ _NODE_TYPE_ ].empty() )
 448     {
 449       // store converted connectivity in data
 450       data._connectivity.insert( data._connectivity.end(), other_conn.begin(), other_conn.end());
 451       if ( need_index )
 452       {
 453         if ( data._index.empty() && index[0] == 1 )
 454           data._index.insert( data._index.end(), index, index + nb_elem + 1 );
 455         else
 456         {
 457           data._index.push_back( 1 );
 458           for ( int i = 0; i < nb_elem; ++i )
 459             data._index.push_back( data._index.back() + index[i+1] - index[i] );
 460         }
 461       }
 462     }
 463     else
 464     {
 465       // exclude coincident elements from connectivity
 466
 467       if ( need_index && data._index.empty() )
 468         data._index.push_back( 1 );
 469
 470       // find added elements possibly(!) coincident with old ones
 471       vector<int>& merged = _merged_of_type[ type ]; // to fill in
 472       int first_added_node = _nb_index[ INIT_OLD ][ MED_NODE ][ _NODE_TYPE_ ] + 1;
 473       for ( int i = 0; i < nb_elem; ++i )
 474       {
 475         // count coincident nodes
 476         int nb_coincident_nodes = 0;
 477         for ( int n = index[i]; n < index[i+1]; ++n )
 478           nb_coincident_nodes += int( other_conn[ n-1 ] < first_added_node );
 479
 480         if ( nb_coincident_nodes == index[i+1] - index[i] )
 481           merged.push_back( i );
 482       }
 483       // find old elements equal to merged, if no equal exist there is zero in array
 484       vector<int>& equalo = _equalo_of_type[ type ];
 485       findEqualOldElements( entity, type, equalo );
 486       if ( equalo.size() < merged.size() )
 487         equalo.resize( merged.size(), 0 );
 488
 489       // fill connectivity
 490       int rm_i = 0, nb_rm = 0;
 491       for ( int i = 0; i < nb_elem; ++i )
 492       {
 493         bool is_merged = ( rm_i < merged.size() && i == merged[rm_i] && equalo[rm_i] );
 494         if ( is_merged )
 495         {
 496           data._nb_elems--;
 497           rm_i++;
 498           nb_rm++;
 499         }
 500         else
 501         {
 502           for ( int n = index[i]; n < index[i+1]; ++n )
 503             data._connectivity.push_back( other_conn[ n-1 ] );
 504           if ( need_index )
 505             data._index.push_back( data._index.back() + index[i+1] - index[i] );
 506         }
 507       }
 508       if ( nb_rm == 0 )
 509         merged.clear(), equalo.clear();
 510     }
 511   }
 512   return data._nb_elems - old_size;
 513 }
 514
 515 //================================================================================
 516 /*!
 517  * \brief Return updated old support
 518  */
 519 //================================================================================
 520
 521 template< class TSUPPORT >
 522 TSUPPORT* MeshFuse::updateOldSupport(TSUPPORT* support) const
 523 {
 524   if ( support->isOnAllElements() )
 525   {
 526     //support->update(); -- this changes old nb elems to nb elems after fuse
 527   }
 528   else if ( support->getEntity() != MED_NODE )
 529   {
 530     // element numbers of partial support change if some elements had been added
 531     // before a type of the entity
 532     const medGeometryElement* types = support->getTypes();
 533     int                    nb_types = support->getNumberOfTypes();
 534     MEDSKYLINEARRAY *        number = support->getnumber();
 535     for ( int t = 0; t < nb_types; ++t )
 536     {
 537       int shift = getElemNbShift( support->getEntity(), types[t], RSLT_ADD, /*prev=*/true);
 538       if ( shift == 0 ) continue;
 539       int nb_elems = support->getNumberOfElements( types[t] );
 540       for ( int j = 0; j < nb_elems; ++j )
 541         number->setIJ( t+1, j+1, number->getIJ( t+1, j+1 ) + shift );
 542     }
 543   }
 544   return support;
 545 }
 546
 547
 548 //================================================================================
 549 /*!
 550  * \brief Creates a copy of support being added or extands the present one
 551  *  \param add_support - support of the added mesh
 552  *  \param same_name_support - the present support with the same name
 553  *  \retval SUPPORT* - resulting support
 554  */
 555 //================================================================================
 556
 557 template< class TSUPPORT >
 558 TSUPPORT* MeshFuse::makeSupport(const TSUPPORT* add_support,
 559                                 TSUPPORT*       same_name_support)
 560 {
 561   if ( same_name_support && same_name_support->getEntity() != add_support->getEntity() )
 562     throw MEDEXCEPTION("MeshFuse: Supports with equal names and different entity!");
 563
 564   // make resulting support
 565
 566   TSUPPORT* res_support = same_name_support;
 567   if ( !same_name_support )
 568   {
 569     res_support = new TSUPPORT;
 570     res_support->setName       ( add_support->getName() );
 571     res_support->setDescription( add_support->getDescription() );
 572     res_support->setEntity     ( add_support->getEntity() );
 573     res_support->setMesh       ( this );
 574   }
 575   else if ( same_name_support->isOnAllElements() && add_support->isOnAllElements() )
 576   {
 577     same_name_support->update();
 578     return same_name_support;
 579   }
 580
 581   // count nb of types
 582
 583   set<medGeometryElement> add_types, old_types, all_types;
 584   if ( res_support->getEntity() == MED_NODE )
 585     add_types.insert( _NODE_TYPE_ );
 586   else
 587     add_types.insert( add_support->getTypes(),
 588                       add_support->getTypes() + add_support->getNumberOfTypes() );
 589   all_types = add_types;
 590   if ( same_name_support )
 591   {
 592     if ( same_name_support->getEntity() == MED_NODE )
 593       old_types.insert( _NODE_TYPE_ );
 594     else
 595       old_types.insert(same_name_support->getTypes(),
 596                        same_name_support->getTypes()+same_name_support->getNumberOfTypes() );
 597     all_types.insert( old_types.begin(), old_types.end() );
 598   }
 599
 600   // make all data
 601
 602   int nb_geom_types = all_types.size();
 603   vector< medGeometryElement > geom_types( nb_geom_types );
 604   vector< vector< int > >        elements( nb_geom_types );
 605   vector< int >               nb_elements( nb_geom_types );
 606   vector< int >                     index( 1,1 );
 607   set<medGeometryElement>::iterator type = all_types.begin();
 608   for ( int t = 0; type != all_types.end(); ++type, ++t )
 609   {
 610     uniteSupportElements( add_types.count( *type ) ? add_support       : 0,
 611                           old_types.count( *type ) ? same_name_support : 0,
 612                           *type,
 613                           elements[ t ]);
 614     nb_elements[ t ] = elements[ t ].size();
 615     index.push_back( index.back() + nb_elements[ t ] );
 616   }
 617
 618   // set data to support
 619
 620   res_support->setAll( false );
 621   res_support->setNumberOfGeometricType( nb_geom_types );
 622   res_support->setGeometricType( &geom_types[0] );
 623   res_support->setNumberOfElements( &nb_elements[0] );
 624   res_support->setNumber (new MEDSKYLINEARRAY( nb_geom_types, index[nb_geom_types]-1));
 625   res_support->getnumber()->setIndex( & index[0] );
 626   for ( int t = 0; t < nb_geom_types; ++t )
 627     res_support->getnumber()->setI( t+1, & elements[t][0] );
 628
 629   return res_support;
 630 }
 631
 632 //================================================================================
 633 /*!
 634  * \brief Concatenate families and groups
 635  */
 636 //================================================================================
 637
 638 void MeshFuse::expandSupports()
 639 {
 640   // we unite supports hiving same names
 641
 642   // make maps of updated old supports
 643   map< string, FAMILY* > name_to_family;
 644   map< string, GROUP* >  name_to_group;
 645   set<int> family_ids;
 646
 647   vector<FAMILY*>* fams[4] = { &_familyNode, &_familyCell, &_familyFace, &_familyEdge};
 648   for ( int i = 0; i < 4; ++i )
 649     for ( int f = 0; f < fams[i]->size(); ++f )
 650     {
 651       name_to_family.insert
 652         ( make_pair( fams[i]->at(f)->getName(), updateOldSupport( fams[i]->at(f) )));
 653       family_ids.insert( fams[i]->at(f)->getIdentifier() );
 654     }
 655
 656   vector<GROUP*>* groups[4] = { &_groupNode, &_groupCell, &_groupFace, &_groupEdge };
 657   for ( int i = 0; i < 4; ++i )
 658     for ( int g = 0; g < groups[i]->size(); ++g )
 659     {
 660       name_to_group.insert
 661         ( make_pair( groups[i]->at(g)->getName(), updateOldSupport( groups[i]->at(g) )));
 662     }
 663
 664   // unite supports
 665   vector<FAMILY*> add_fams[4]={ _mesh->getFamilies(MED_NODE),
 666                                 _mesh->getFamilies(MED_CELL),
 667                                 _mesh->getFamilies(MED_FACE),
 668                                 _mesh->getFamilies(MED_EDGE) };
 669   for ( int i = 0; i < 4; ++i )
 670     for ( int f = 0; f < add_fams[i].size(); ++f )
 671     {
 672       FAMILY* add_fam = add_fams[i][f];
 673       map<string,FAMILY*>::iterator name_fam = name_to_family.find( add_fam->getName());
 674       FAMILY* res_fam = makeSupport( add_fam,
 675                                      name_fam == name_to_family.end() ? 0 : name_fam->second);
 676       if ( name_to_family.insert( make_pair( res_fam->getName(), res_fam )).second )
 677       {
 678         fams[ i ]->push_back( res_fam );
 679         int id = add_fam->getIdentifier();
 680         if ( family_ids.count ( id ))
 681           id = ( id < 0 ) ? *family_ids.begin() - 1 : *family_ids.rbegin() + 1;
 682         res_fam->setIdentifier( id );
 683         family_ids.insert( id );
 684       }
 685       // update group names in a family
 686       vector<string> res_gr_names( res_fam->getGroupsNames(),
 687                                    res_fam->getGroupsNames()+res_fam->getNumberOfGroups());
 688       for ( int g = 0; g < add_fam->getNumberOfGroups(); ++g )
 689         if ( find( res_gr_names.begin(), res_gr_names.end(), add_fam->getGroupName(g+1))
 690              == res_gr_names.end())
 691           res_gr_names.push_back( add_fam->getGroupName( g+1 ));
 692       if ( res_fam->getNumberOfGroups() < res_gr_names.size() )
 693       {
 694         res_fam->setNumberOfGroups( res_gr_names.size() );
 695         res_fam->setGroupsNames( &res_gr_names[0] );
 696       }
 697     }
 698
 699   vector<GROUP*> add_grp[4] ={ _mesh->getGroups(MED_NODE),
 700                                _mesh->getGroups(MED_CELL),
 701                                _mesh->getGroups(MED_FACE),
 702                                _mesh->getGroups(MED_EDGE) };
 703   for ( int i = 0; i < 4; ++i )
 704   {
 705     for ( int g = 0; g < add_grp[i].size(); ++g )
 706     {
 707       map< string, GROUP* >::iterator name_grp = name_to_group.find( add_grp[i][g]->getName());
 708 //       STRING out("OLD GRP: ");
 709 //       if ( name_grp == name_to_group.end() ) out << "none";
 710 //       else out << *(name_grp->second);
 711 //       out._s << endl;
 712 //       out << "ADD GRP: " << *(add_grp[i][g]); out._s << endl; out << "";
 713       GROUP* res_grp = makeSupport( add_grp[i][g],
 714                                     name_grp == name_to_group.end() ? 0 : name_grp->second);
 715       if ( !name_to_group.count( res_grp->getName() ))
 716         groups[ i ]->push_back( res_grp );
 717       if ( res_grp->getFamilies().size() < add_grp[i][g]->getFamilies().size() )
 718       {
 719         res_grp->setFamilies( add_grp[i][g]->getFamilies() );
 720         res_grp->setNumberOfFamilies( add_grp[i][g]->getNumberOfFamilies() );
 721       }
 722 //       out << "RES GRP: " << *res_grp;
 723 //       cout << out << endl;
 724     }
 725     // update pointers to families in groups
 726     for ( int g = 0; g < groups[i]->size(); ++g )
 727     {
 728       GROUP* grp = groups[i]->at(g);
 729       if ( grp->getNumberOfFamilies() < 1 ) continue;
 730       vector<FAMILY*> fams = grp->getFamilies();
 731       for ( int f = 0; f < fams.size(); ++f )
 732         fams[f] = name_to_family[ fams[f]->getName() ];
 733       grp->setFamilies( fams );
 734       grp->setNumberOfFamilies( fams.size() );
 735     }
 736   }
 737 }
 738
 739 //================================================================================
 740 /*!
 741  * \brief Return shift for conversion of element numbers
 742  *  \param which - to select a global numbering index
 743  *  \param prev - true means "shift of type previous to given type"
 744  */
 745 //================================================================================
 746
 747 int MeshFuse::getElemNbShift( const medEntityMesh& entity,
 748                               medGeometryElement   type,
 749                               int which, bool prev ) const
 750 {
 751   //if ( type == MED_NONE ) type = MED_POINT1;
 752
 753   const TNbOfGeom & shift_of_type = _nb_index[ which ][ int(entity) ];
 754
 755   TNbOfGeom::const_iterator type_shift = shift_of_type.lower_bound( type );
 756   if ( type_shift == shift_of_type.end() )
 757     return shift_of_type.empty() ? 0 : shift_of_type.rbegin()->second;
 758
 759   // get shift of exactly the type or of the previos type
 760   if ( ( prev && type_shift->first == type ) || type_shift->first > type )
 761   {
 762     if ( type_shift == shift_of_type.begin() )
 763       return 0;
 764     else
 765       type_shift--;
 766   }
 767
 768   return type_shift->second;
 769 }
 770
 771 //================================================================================
 772 /*!
 773  * \brief Fills in elements of support of given type
 774  *  \param add_support - support of the added mesh
 775  *  \param old_support - the present support with the same name
 776  *  \param type - geometric type to process
 777  *  \param elements - output array of element
 778  */
 779 //================================================================================
 780
 781 void MeshFuse::uniteSupportElements(const SUPPORT*     add_support,
 782                                     SUPPORT*           old_support,
 783                                     medGeometryElement type,
 784                                     vector<int> &      elements)
 785 {
 786   int sup_type = ( type/100 == 0 ? MED_ALL_ELEMENTS : type );
 787
 788   const medEntityMesh entity = (add_support ? add_support : old_support )->getEntity();
 789
 790   // Add old elements
 791
 792   int nb_add_elems = add_support ? add_support->getNumberOfElements( sup_type ) : 0;
 793   int nb_old_elems = 0;
 794   if ( old_support )
 795   {
 796     nb_old_elems = old_support->getNumberOfElements( sup_type );
 797     elements.reserve( nb_old_elems + nb_add_elems );
 798     int add_shift = getElemNbShift( entity, type, RSLT_ADD, /*prev=*/1);
 799     int old_shift = getElemNbShift( entity, type, INIT_OLD, /*prev=*/1);
 800     int shift = 1 + add_shift + old_shift;
 801     // convertion of numbers is already done in updateOldSupport()
 802     if ( old_support->isOnAllElements() )
 803       for ( int i = 0; i < nb_old_elems; ++i )
 804         elements.push_back( i + shift );
 805     else
 806       elements.insert( elements.end(), old_support->getNumber( sup_type ),
 807                        old_support->getNumber( sup_type ) + nb_old_elems );
 808     if ( nb_add_elems == 0 )
 809       return;
 810   }
 811   else
 812   {
 813     elements.reserve( nb_add_elems );
 814   }
 815
 816   // Convert and add elements of the added support
 817
 818   const int * add_elems = add_support->isOnAllElements() ? 0 : add_support->getNumber( sup_type );
 819
 820   int add_shift = getElemNbShift( entity, type, RSLT_ADD, /*prev=*/1);
 821   int old_shift = getElemNbShift( entity, type, INIT_OLD, /*prev=*/0);
 822
 823   if ( _merged_of_type[ type ].empty() )
 824   {
 825     // no coicident elements in the old and the added meshes - just unite two lists
 826     int shift = add_support->isOnAllElements() ? 1 + add_shift + old_shift : old_shift;
 827     if ( add_support->isOnAllElements() )
 828       for ( int i = 0; i < nb_add_elems; ++i )
 829         elements.push_back( i + shift );
 830     else
 831       for ( int i = 0; i < nb_add_elems; ++i )
 832         elements.push_back( add_elems[i] + shift );
 833   }
 834   else
 835   {
 836     // some elements of the added mesh coincide with old ones,
 837     // so conversion is not so easy, and some support members can
 838     // be twice in it - skip them
 839     vector<int> & new_elem_ids = _new_elem_ids_of_type[ type ];
 840     if ( new_elem_ids.empty() )
 841       makeNewElemIds( entity, type, new_elem_ids );
 842
 843     set< int > old_elems_in_support( elements.begin(), elements.end() );
 844     int last_old_nb = old_shift + add_shift;
 845     int shift = 1 + getElemNbShift( entity, type, INIT_ADD, /*prev=*/1);
 846     for ( int i = 0; i < nb_add_elems; ++i )
 847     {
 848       int new_id = new_elem_ids[ add_elems ? add_elems[i]-shift : i ];
 849       if ( new_id > last_old_nb || old_elems_in_support.count( new_id ) == 0 )
 850         elements.push_back( new_id );
 851     }
 852   }
 853 }
 854
 855 //================================================================================
 856 /*!
 857  * \brief Fills in ids of elements of added mesh in the resulting mesh
 858  *  \param type - element type to treat
 859  *  \param new_ids - output array
 860  */
 861 //================================================================================
 862
 863 void MeshFuse::makeNewElemIds(medEntityMesh      entity,
 864                               medGeometryElement type,
 865                               vector< int > &    new_ids)
 866 {
 867   const char* LOC = "MeshFuse::makeNewElemIds(): ";
 868   if ( entity == MED_NODE )
 869     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"we must not be here"));
 870
 871   vector<int>& merged_i = _merged_of_type[ type ];
 872   vector<int>::iterator rm_i = merged_i.begin();
 873
 874   // find ids for merged added elements
 875   vector< int >& old_ids = _equalo_of_type[ type ];
 876 //   if ( old_ids.empty() && !merged_i.empty() )
 877 //     findEqualOldElements( entity, type, old_ids );
 878   vector< int >::iterator old_id = old_ids.begin();
 879
 880   // nb of added elements
 881   int add_nb_elems = _mesh->getNumberOfElements( entity, type );
 882   new_ids.resize( add_nb_elems );
 883
 884   // new id of the first added element
 885   int old_shift = getElemNbShift( entity, type, INIT_OLD, /*prev=*/0);
 886   int add_shift = getElemNbShift( entity, type, RSLT_ADD, /*prev=*/1);
 887   int elem_id = old_shift + add_shift + 1;
 888
 889   // all new ids
 890   // (this works implying that numbers in SUPPORT are in inceasing order!)
 891   for (int i_elem = 0; i_elem < add_nb_elems; )
 892   {
 893     int i_until = ( rm_i == merged_i.end() ? add_nb_elems : *rm_i++ );
 894     // increase elem_id until the next merged element
 895     while ( i_elem < i_until )
 896       new_ids[ i_elem++ ] = elem_id++;
 897     // insert id of old element equal to merged one, if any
 898     if ( i_elem < add_nb_elems )
 899     {
 900       if ( *old_id )
 901         new_ids[ i_elem++ ] = *old_id++;
 902       else
 903         new_ids[ i_elem++ ] = elem_id++, old_id++; // no equal old elem exist
 904     }
 905   }
 906 }
 907
 908 //================================================================================
 909 /*!
 910  * \brief Finds ids of elements of the old mesh equal to merged elements of the added one
 911  */
 912 //================================================================================
 913
 914 void MeshFuse::findEqualOldElements(medEntityMesh      entity,
 915                                     medGeometryElement type,
 916                                     vector< int > &    old_ids)
 917 {
 918   // poly element can coincide with any type of the same entity
 919   const bool isPoly = ( type == MED_POLYGON || type == MED_POLYHEDRA );
 920   const medGeometryElement checkType = isPoly ? MED_ALL_ELEMENTS : type;
 921
 922   if ( !_mesh->getNumberOfElements(entity, type) ||
 923        ! this->getNumberOfElements(entity, checkType) )
 924     return;
 925
 926   int old_nb_elems_end, old_nb_elems_start;
 927   if ( isPoly )
 928     {
 929       old_nb_elems_start = 0;
 930       old_nb_elems_end   = this->getNumberOfElements( entity, MED_ALL_ELEMENTS );
 931     }
 932   else
 933     {
 934       // if this method is called when connectivity of the entity is not yet concatenated:
 935       old_nb_elems_start = getElemNbShift( entity, type, INIT_OLD, /*prev=*/true);
 936       old_nb_elems_end   = getElemNbShift( entity, type, INIT_OLD, /*prev=*/false);
 937       // if this method is called after expandConnectivity() and this mesh contains all elements
 938       //   int old_nb_elems =
 939     }
 940   const int *old_conn, *old_index, *add_conn, *add_index;
 941   add_conn  = _mesh->getConnectivity( MED_NODAL, entity, type );
 942   old_conn  =  this->getConnectivity( MED_NODAL, entity, checkType );
 943   add_index = _mesh->getConnectivityIndex( MED_NODAL, entity );
 944   old_index =  this->getConnectivityIndex( MED_NODAL, entity );
 945
 946   const vector<int>& new_node_ids = _new_elem_ids_of_type[ _NODE_TYPE_ ];
 947
 948   const vector<int>& merged_i = _merged_of_type[ type ];
 949   vector<int>::const_iterator rm_i = merged_i.begin();
 950
 951   old_ids.reserve( merged_i.size() );
 952
 953   for ( ; rm_i != merged_i.end(); ++rm_i ) // rm_i counts from 0
 954   {
 955     // get nodes of rm_i-th added face
 956     const int* add_elem_conn = add_conn + add_index[ *rm_i ]-1;
 957     int    nb_add_elem_nodes = add_index[ *rm_i+1 ] - add_index[ *rm_i ];
 958     set<int> add_elem_nodes;
 959     for ( int n = 0; n < nb_add_elem_nodes; ++n )
 960       add_elem_nodes.insert( new_node_ids[ add_elem_conn[n]-1 ]);
 961
 962     // look for equal old face among all old faces
 963     const int* old_node = old_conn;
 964     int        new_node = *add_elem_nodes.begin();
 965     int  found_old_elem = 0;
 966     for ( int i = old_nb_elems_start; i < old_nb_elems_end && !found_old_elem; ++i )
 967     {
 968       int nb_old_elem_nodes = old_index[ i+1 ] - old_index[ i ];
 969       for ( int j = 0; j < nb_old_elem_nodes; ++j, ++old_node )
 970       {
 971         if ( new_node != *old_node ) continue;
 972         // found common node, compare all nodes
 973         const int* old_elem_conn = old_conn + old_index[ i ]-1;
 974         set<int> old_elem_nodes( old_elem_conn, old_elem_conn + nb_old_elem_nodes);
 975         if ( add_elem_nodes == old_elem_nodes )
 976         {
 977           found_old_elem = i + 1;
 978           break;
 979         }
 980       }
 981     }
 982     // Issue 0020809. Its possible that no old element exists
 983 //     if ( !found_old_elem )
 984 //       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"all nodes of added elements are merged "
 985 //                                    "but equal old element not found. "
 986 //                                    "Please check global nodes numbers."));
 987     old_ids.push_back( found_old_elem );
 988   }
 989 }
 990
 991 //================================================================================
 992 /*!
 993  * \brief Unite numbers of elements taking into account their types
 994  */
 995 //================================================================================
 996
 997 void MeshFuse::append( medEntityMesh      entity,
 998                        vector<int>&       numbers,
 999                        const vector<int>& add_numbers )
1000 {
1001   const char* LOC="MeshFuse::append(): ";
1002
1003   int nb_types = getNumberOfTypes( entity );
1004   if ( numbers.empty() || add_numbers.empty() ||
1005        ( nb_types < 2 && _merged_of_type[ getElementType( entity, 1 )].empty() ))
1006   {
1007     numbers.insert( numbers.end(), add_numbers.begin(), add_numbers.end() );
1008     return;
1009   }
1010
1011   vector<int> result;
1012   result.reserve( numbers.size() + add_numbers.size() );
1013   const int* old_nums = & numbers[0];
1014   const int* add_nums = & add_numbers[0];
1015
1016   MEDMEM::PointerOf<medGeometryElement> types;
1017   types.setShallowAndOwnership( getTypes(entity));
1018
1019   for ( int t = 0; t < nb_types; ++t )
1020   {
1021     int nb_old =
1022       getElemNbShift( entity, types[t], INIT_OLD, /*prev=*/false) -
1023       getElemNbShift( entity, types[t], INIT_OLD, /*prev=*/true);
1024     int nb_add =
1025       getElemNbShift( entity, types[t], INIT_ADD, /*prev=*/false) -
1026       getElemNbShift( entity, types[t], INIT_ADD, /*prev=*/true);
1027
1028     result.insert( result.end(), old_nums, old_nums + nb_old );
1029     old_nums += nb_old;
1030
1031     const vector<int>& loc_merged = _merged_of_type[ types[t] ];
1032     if ( loc_merged.empty() )
1033     {
1034       result.insert( result.end(), add_nums, add_nums + nb_add );
1035     }
1036     else
1037     {
1038       vector<int>::const_iterator imerged = loc_merged.begin();
1039       vector<int>::const_iterator equalid = _equalo_of_type[ types[t] ].begin();
1040
1041       int from = 0, to = -1;
1042       while ( from < nb_add )
1043       {
1044         while ( imerged != loc_merged.end() && *equalid == 0 )
1045           imerged++, equalid++;
1046         to = ( imerged == loc_merged.end() ? nb_add : ( equalid++, *imerged++ ));
1047         if ( to > from )
1048           result.insert( result.end(), add_nums + from, add_nums + to );
1049         from = to + 1;
1050       }
1051     }
1052     add_nums += nb_add;
1053   }
1054   if ( result.size() != getNumberOfElements( entity, MED_ALL_ELEMENTS ))
1055     throw MED_EXCEPTION(MEDMEM::STRING(LOC) << "invalid nb of numbers of entity " << entity
1056                         << ": expect " << getNumberOfElements( entity, MED_ALL_ELEMENTS)
1057                         << " but get " << result.size());
1058
1059   numbers.swap(result);
1060 }