Back to index

salome-med  6.5.0
MEDMEM_MeshFuse.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
00002 //
00003 // This library is free software; you can redistribute it and/or
00004 // modify it under the terms of the GNU Lesser General Public
00005 // License as published by the Free Software Foundation; either
00006 // version 2.1 of the License.
00007 //
00008 // This library is distributed in the hope that it will be useful,
00009 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00010 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00011 // Lesser General Public License for more details.
00012 //
00013 // You should have received a copy of the GNU Lesser General Public
00014 // License along with this library; if not, write to the Free Software
00015 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00016 //
00017 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00018 //
00019 // File      : MEDMEM_MeshFuse.cxx
00020 // Created   : Tue Jul  7 18:27:00 2009
00021 // Author    : Edward AGAPOV (eap)
00022 
00023 #include "MEDMEM_MeshFuse.hxx"
00024 
00025 #include "MEDMEM_Family.hxx"
00026 #include "MEDMEM_Group.hxx"
00027 
00028 using namespace MEDMEM;
00029 using namespace MED_EN;
00030 using namespace std;
00031 
00032 #define _NODE_TYPE_ MED_NONE
00033 
00034 MeshFuse::MeshFuse():MESHING()
00035 {
00036 }
00037 
00038 MeshFuse::~MeshFuse()
00039 {
00040 }
00041 
00042 //================================================================================
00046 //================================================================================
00047 
00048 void MeshFuse::setNodeNumbers( const std::vector<int>& node_glob_numbers )
00049 {
00050   const char* LOC = "MeshFuse::setNodeNumbers(node_glob_numbers): ";
00051 
00052   if ( !_node_glob_numbers.empty() )
00053     throw MEDEXCEPTION(STRING(LOC)<<"node numbers has been already set");
00054 
00055   if ( node_glob_numbers.size() != getNumberOfNodes() &&
00056        node_glob_numbers.size() > 0 && getNumberOfNodes() > 0 )
00057     throw MEDEXCEPTION
00058       (STRING(LOC)<<"size of node_glob_numbers must be equal number of nodes in MeshFuse");
00059 
00060   _node_glob_numbers = node_glob_numbers;
00061 }
00062 
00063 //================================================================================
00069 //================================================================================
00070 
00071 void MeshFuse::concatenate( const MESH* mesh, const vector<int>& node_glob_numbers )
00072 {
00073   const char* LOC = "MeshFuse::concatenate( mesh, node_glob_numbers ): ";
00074   if ( !mesh || mesh->getNumberOfNodes() == 0) return;
00075     
00076   _mesh = mesh;
00077 
00078   if ( this->getNumberOfNodes() < 1 )
00079   {
00080     mesh->getCoordinates( MED_FULL_INTERLACE );// make fields filled in case of GRID
00081     mesh->getConnectivityptr();
00082 
00083     // just copy mesh
00084     static_cast<MESH&>(*this) = *mesh;
00085 
00086     _node_glob_numbers = node_glob_numbers;
00087     return;
00088   }
00089   // check feasibility
00090 
00091   if ( mesh->getNumberOfNodes() > 0 && node_glob_numbers.empty() )
00092     throw MEDEXCEPTION(STRING(LOC)<<"merging without node global numbers not implemented yet");
00093 
00094   if ( mesh->getNumberOfNodes() != node_glob_numbers.size() )
00095     throw MEDEXCEPTION(STRING(LOC)<<"invalid number of node global numbers");
00096 
00097   if ( mesh->getSpaceDimension() != this->getSpaceDimension() ||
00098        mesh->getMeshDimension()  != this->getMeshDimension() )
00099     throw MEDEXCEPTION(STRING(LOC)<<"can't unite meshes with different dimensions so far, sorry");
00100   
00101   // init
00102   _merged_of_type.clear();
00103   for ( int i = 0; i < NB_INDICES; ++i )
00104   {
00105     _nb_index[ i ].clear();
00106     _nb_index[ i ].resize( MED_ALL_ENTITIES );
00107   }
00108 
00109   // concatenation
00110 
00111   int final_nb_nodes = makeNewNodeIds( node_glob_numbers );
00112 
00113   expandCoordinates(final_nb_nodes);
00114 
00115   expandConnectivity(final_nb_nodes);
00116 
00117   expandSupports();
00118 
00119   // clear unnecessary data
00120   _new_elem_ids_of_type.clear();
00121 }
00122 
00123 //================================================================================
00128 //================================================================================
00129 
00130 int MeshFuse::makeNewNodeIds(const vector<int>& add_glob_numbers)
00131 {
00132   // remember merged added nodes and make an array of new ids of added nodes
00133   vector<int>&       merged = _merged_of_type      [ _NODE_TYPE_ ];
00134   vector<int>& new_node_ids = _new_elem_ids_of_type[ _NODE_TYPE_ ];
00135   new_node_ids.resize( add_glob_numbers.size() );
00136 
00137   // extand global node numbers
00138   _node_glob_numbers.reserve( _node_glob_numbers.size() + add_glob_numbers.size());
00139 
00140   map<int,int> glob_ids; // to map global ids to new local ones
00141   for ( int i = 0; i < _node_glob_numbers.size(); ++i )
00142     glob_ids.insert( make_pair( _node_glob_numbers[i], i+1 ));
00143 
00144   int next_loc_id = getNumberOfNodes() + 1;
00145 
00146   for ( int i = 0; i < add_glob_numbers.size(); ++i )
00147   {
00148     map<int,int>::iterator  glob_loc =
00149       glob_ids.insert( make_pair( add_glob_numbers[i], next_loc_id )).first;
00150 
00151     new_node_ids[i] = glob_loc->second;
00152 
00153     if ( new_node_ids[i] == next_loc_id ) // unique global number
00154     {
00155       next_loc_id++;
00156       _node_glob_numbers.push_back( add_glob_numbers[i] );
00157     }
00158     else
00159     {
00160       merged.push_back( i );
00161     }
00162   }
00163   _nb_index[ INIT_OLD ][ MED_NODE ][ _NODE_TYPE_ ] = getNumberOfNodes();
00164   _nb_index[ INIT_ADD ][ MED_NODE ][ _NODE_TYPE_ ] = add_glob_numbers.size();
00165   _nb_index[ RSLT_ADD ][ MED_NODE ][ _NODE_TYPE_ ] = add_glob_numbers.size() - merged.size();
00166 
00167   return next_loc_id - 1;
00168 }
00169 
00170 //================================================================================
00174 //================================================================================
00175 
00176 void MeshFuse::expandCoordinates(int final_nb_nodes)
00177 {
00178   const int dim = getSpaceDimension();
00179 
00180   // create new coordinates
00181   double* coord = new double[ final_nb_nodes * dim ];
00182   MEDARRAY<double> medarray( coord, dim, final_nb_nodes, MED_FULL_INTERLACE,
00183                              /*shallowCopy=*/true,/*ownershipOfValues=*/false);
00184   // copy old coordinates
00185   int nb_old_coord = getNumberOfNodes() * dim;
00186   memcpy( coord, getCoordinates( MED_FULL_INTERLACE ), nb_old_coord * sizeof( double ));
00187 
00188   // set new coords
00189   coord += nb_old_coord;
00190   const double* add_coord =_mesh->getCoordinates( MED_FULL_INTERLACE );
00191   if ( _merged_of_type[ _NODE_TYPE_ ].empty())
00192   {
00193     // no coincident nodes in the two meshes, just add coords
00194     memcpy( coord, add_coord, _mesh->getNumberOfNodes() * dim * sizeof( double ));
00195   }
00196   else
00197   {
00198     // copy coord of only unique nodes
00199     int first_added_node = getNumberOfNodes() + 1;
00200     const vector<int>& new_node_ids = _new_elem_ids_of_type[ _NODE_TYPE_ ];
00201     for ( int n = 0; n < new_node_ids.size(); ++n )
00202     {
00203       if ( new_node_ids[n] < first_added_node ) continue; // coincident node
00204       memcpy( coord, add_coord + n * dim, dim * sizeof( double ));
00205       coord += dim;
00206     }
00207   }
00208   _coordinate->setCoordinates( &medarray, /*shallowCopy=*/true );
00209 
00210   _numberOfNodes = final_nb_nodes;
00211 }
00212 
00213 //================================================================================
00219 //================================================================================
00220 
00221 void MeshFuse::expandConnectivity(int final_nb_nodes)
00222 {
00223   const medConnectivity nodal   = MED_NODAL;
00224 
00225   // fill in _nb_index[ INIT_OLD ]
00226   for ( medEntityMesh entity = MED_CELL; entity < MED_ALL_ENTITIES; ++entity )
00227   {
00228     if ( existConnectivity( nodal, entity ))
00229     {
00230       const int *               index = this->getGlobalNumberingIndex(entity);
00231       const medGeometryElement* types = this->getTypes(entity);
00232       int                    nb_types = this->getNumberOfTypes(entity);
00233       for ( int t = 0; t < nb_types; ++t )
00234         _nb_index[ INIT_OLD ][ entity ][ types[t] ] = index[t+1]-index[0];
00235     }
00236   }
00237 
00238   CONNECTIVITY *prev_connectivity = 0, *cell_connectivity = 0;
00239 
00240   // loop on all entities
00241   for ( medEntityMesh entity = MED_CELL; entity < MED_ALL_ENTITIES; ++entity )
00242   {
00243     if ( entity == MED_FACE && getMeshDimension() == 2 )
00244       continue; // there can be EDGE
00245 
00246     if ( !_mesh->existConnectivity( nodal, entity ))
00247     {
00248       // no entity in added mesh
00249       if ( existConnectivity( nodal, entity ) && prev_connectivity )
00250         prev_connectivity->setConstituent
00251           ( new CONNECTIVITY( *getConnectivityptr()->getConstituent( entity )));
00252       break;
00253     }
00254     if ( !existConnectivity( nodal, entity ))
00255     {
00256       // no entity in the old mesh - fully copy the added connectivity
00257       CONNECTIVITY* connectivity =
00258         new CONNECTIVITY( *_mesh->getConnectivityptr()->getConstituent( entity ));
00259       connectivity->setNumberOfNodes( final_nb_nodes );
00260       updateNodeIds( connectivity );
00261       if ( entity == MED_CELL )
00262         cell_connectivity = connectivity;
00263       else
00264         cell_connectivity->setConstituent( connectivity );
00265 
00266       break;
00267     }
00268     // Fill in _nb_index[ INIT_ADD ]
00269     const int *                   index = _mesh->getGlobalNumberingIndex(entity);
00270     const medGeometryElement* add_types = _mesh->getTypes(entity);
00271     int                    nb_add_types = _mesh->getNumberOfTypes(entity);
00272     for ( int t = 0; t < nb_add_types; ++t )
00273       _nb_index[ INIT_ADD ][ entity ][ add_types[t] ] = index[t+1]-index[0];
00274 
00275     // Unite connectivities of std types
00276 
00277     // count types
00278     set<medGeometryElement> types;
00279     types.insert( this->getTypes(entity), this->getTypes(entity) + this->getNumberOfTypes(entity));
00280     types.insert(_mesh->getTypes(entity),_mesh->getTypes(entity) +_mesh->getNumberOfTypes(entity));
00281 
00282     int sum_old = 0, sum_add = 0;
00283 
00284     // make all data
00285     vector< TConnData >          conn_of_type( types.size() );
00286     vector< medGeometryElement > type_vec( types.size() );
00287     vector< int  >               count(1,1);
00288     set<medGeometryElement>::iterator type = types.begin();
00289     for ( int t = 0; type != types.end(); ++type, ++t )
00290     {
00291       int max_conn_len = 0;
00292       if ( this->getNumberOfElements( entity, *type ))
00293         max_conn_len += this->getConnectivityLength( nodal, entity, *type);
00294       if ( _mesh->getNumberOfElements( entity, *type ))
00295         max_conn_len += _mesh->getConnectivityLength( nodal, entity, *type);
00296       conn_of_type[t]._connectivity.reserve( max_conn_len );
00297 
00298       int nb_old = appendConnectivity( conn_of_type[t], this, entity, *type );
00299       sum_old += nb_old;
00300       _nb_index[ INIT_OLD ][ entity ][ *type ] = sum_old;
00301 
00302       int nb_add = appendConnectivity( conn_of_type[t],_mesh, entity, *type );
00303       sum_add += nb_add;
00304       _nb_index[ RSLT_ADD ][ entity ][ *type ] = sum_add;
00305 
00306       count.push_back( count.back() + conn_of_type[t]._nb_elems );
00307       type_vec[t] = *type;
00308     }
00309     // make new connectivity
00310     CONNECTIVITY* connectivity = new CONNECTIVITY( types.size(), entity );
00311     connectivity->setNumberOfNodes ( final_nb_nodes );
00312     connectivity->setGeometricTypes( &type_vec[0], entity );
00313     connectivity->setCount         ( &count   [0], entity );
00314     for ( int t = 0; t < types.size(); ++t )
00315       connectivity->setNodal( &conn_of_type[t]._connectivity[0],
00316                               entity, type_vec[t],
00317                               &conn_of_type[t]._index[0] );
00318 
00319     // store connectivity of an entity
00320     if ( !prev_connectivity )
00321     {
00322       cell_connectivity = connectivity;
00323       prev_connectivity = cell_connectivity;
00324     }
00325     else
00326     {
00327       prev_connectivity->setConstituent( connectivity );
00328       prev_connectivity = connectivity;
00329     }
00330 
00331   } // loop on entities
00332 
00333   if ( cell_connectivity )
00334   {
00335     delete _connectivity;
00336     _connectivity = cell_connectivity;
00337   }
00338 }
00339 
00340 //================================================================================
00344 //================================================================================
00345 
00346 void MeshFuse::updateNodeIds( CONNECTIVITY* connectivity )
00347 {
00348   const medConnectivity   nodal = MED_NODAL;
00349   const medGeometryElement type = MED_ALL_ELEMENTS;
00350 
00351   const vector<int>& new_node_ids = _new_elem_ids_of_type[ _NODE_TYPE_ ];
00352 
00353   medEntityMesh entity = connectivity->getEntity();
00354 
00355   for ( ; entity < MED_ALL_ENTITIES; entity++ )
00356   {
00357     // Collect all connectivities with their lengths
00358 
00359     list< pair< const int*, int > > conn_len_list;
00360 
00361     if ( connectivity->existConnectivity( nodal, entity ))
00362       conn_len_list.push_back
00363         ( make_pair( connectivity->getConnectivity(nodal,entity,type),
00364                      connectivity->getConnectivityLength(nodal,entity,type)));
00365 
00366     // Convert them
00367 
00368     list< pair< const int*, int > >::iterator conn_len = conn_len_list.begin();
00369     for ( ; conn_len != conn_len_list.end(); ++conn_len )
00370     {
00371       int* conn = (int* ) conn_len->first;
00372       int* conn_end = conn + conn_len->second;
00373       for ( ; conn < conn_end; ++conn )
00374         *conn = new_node_ids[ *conn-1 ];
00375     }
00376   }
00377 }
00378 
00379 //================================================================================
00388 //================================================================================
00389 
00390 int MeshFuse::appendConnectivity( MeshFuse::TConnData& data,
00391                                   const MESH*          mesh,
00392                                   medEntityMesh        entity,
00393                                   medGeometryElement   type)
00394 {
00395   // get connectivity of type
00396 
00397   const int* conn, *index = 0;
00398   int nb_elem, conn_len;
00399 
00400   nb_elem = mesh->getNumberOfElements ( entity, type );
00401   if ( !nb_elem ) return 0;
00402   conn    = mesh->getConnectivity     ( MED_NODAL, entity, type );
00403   index   = mesh->getConnectivityIndex( MED_NODAL, entity );
00404   int shift = getElemNbShift( entity, type, (mesh == this ? INIT_OLD : INIT_ADD), /*prev=*/true);
00405   index += shift;
00406   conn_len = index[ nb_elem ] - index[ 0 ];
00407 
00408   bool need_index = ( type == MED_POLYGON || type == MED_POLYHEDRA );
00409   if ( !need_index )
00410     data._index.resize( 1, 0 ); // for safe access to pointer even if no real index exists
00411 
00412   // Append
00413 
00414   int old_size = data._nb_elems;
00415   data._nb_elems += nb_elem;
00416 
00417   if ( mesh == this )
00418   {
00419     // append connectivity to data as is
00420     data._connectivity.insert( data._connectivity.end(), conn, conn + conn_len );
00421     if ( need_index )
00422     {
00423       if ( data._index.empty() )
00424         data._index.insert( data._index.end(), index, index + nb_elem + 1 );
00425       else
00426         for ( int i = 0; i < nb_elem; ++i )
00427           data._index.push_back( data._index.back() + index[i+1] - index[i] );
00428     }
00429   }
00430   else
00431   {
00432     // convert connectivity of other mesh
00433 
00434     vector<int> other_conn( conn_len );
00435     const vector<int>& new_node_ids = _new_elem_ids_of_type[ _NODE_TYPE_ ];
00436     if ( type == MED_POLYHEDRA )
00437       {
00438         for ( int n = 0; n < conn_len; ++n )
00439           if ( conn[ n ] > 0 )
00440             other_conn[ n ] = new_node_ids[ conn[ n ]-1 ];
00441       }
00442     else
00443       {
00444         for ( int n = 0; n < conn_len; ++n )
00445           other_conn[ n ] = new_node_ids[ conn[ n ]-1 ];
00446       }
00447     if ( entity == MED_CELL || _merged_of_type[ _NODE_TYPE_ ].empty() )
00448     {
00449       // store converted connectivity in data
00450       data._connectivity.insert( data._connectivity.end(), other_conn.begin(), other_conn.end());
00451       if ( need_index )
00452       {
00453         if ( data._index.empty() && index[0] == 1 )
00454           data._index.insert( data._index.end(), index, index + nb_elem + 1 );
00455         else
00456         {
00457           data._index.push_back( 1 );
00458           for ( int i = 0; i < nb_elem; ++i )
00459             data._index.push_back( data._index.back() + index[i+1] - index[i] );
00460         }
00461       }
00462     }
00463     else
00464     {
00465       // exclude coincident elements from connectivity
00466 
00467       if ( need_index && data._index.empty() )
00468         data._index.push_back( 1 );
00469 
00470       // find added elements possibly(!) coincident with old ones
00471       vector<int>& merged = _merged_of_type[ type ]; // to fill in
00472       int first_added_node = _nb_index[ INIT_OLD ][ MED_NODE ][ _NODE_TYPE_ ] + 1;
00473       for ( int i = 0; i < nb_elem; ++i )
00474       {
00475         // count coincident nodes
00476         int nb_coincident_nodes = 0;
00477         for ( int n = index[i]; n < index[i+1]; ++n )
00478           nb_coincident_nodes += int( other_conn[ n-1 ] < first_added_node );
00479 
00480         if ( nb_coincident_nodes == index[i+1] - index[i] )
00481           merged.push_back( i );
00482       }
00483       // find old elements equal to merged, if no equal exist there is zero in array
00484       vector<int>& equalo = _equalo_of_type[ type ];
00485       findEqualOldElements( entity, type, equalo );
00486       if ( equalo.size() < merged.size() )
00487         equalo.resize( merged.size(), 0 );
00488 
00489       // fill connectivity
00490       int rm_i = 0, nb_rm = 0;
00491       for ( int i = 0; i < nb_elem; ++i )
00492       {
00493         bool is_merged = ( rm_i < merged.size() && i == merged[rm_i] && equalo[rm_i] );
00494         if ( is_merged )
00495         {
00496           data._nb_elems--;
00497           rm_i++;
00498           nb_rm++;
00499         }
00500         else
00501         {
00502           for ( int n = index[i]; n < index[i+1]; ++n )
00503             data._connectivity.push_back( other_conn[ n-1 ] );
00504           if ( need_index )
00505             data._index.push_back( data._index.back() + index[i+1] - index[i] );
00506         }
00507       }
00508       if ( nb_rm == 0 )
00509         merged.clear(), equalo.clear();
00510     }
00511   }
00512   return data._nb_elems - old_size;
00513 }
00514 
00515 //================================================================================
00519 //================================================================================
00520 
00521 template< class TSUPPORT >
00522 TSUPPORT* MeshFuse::updateOldSupport(TSUPPORT* support) const 
00523 {
00524   if ( support->isOnAllElements() )
00525   {
00526     //support->update(); -- this changes old nb elems to nb elems after fuse
00527   }
00528   else if ( support->getEntity() != MED_NODE )
00529   {
00530     // element numbers of partial support change if some elements had been added
00531     // before a type of the entity
00532     const medGeometryElement* types = support->getTypes();
00533     int                    nb_types = support->getNumberOfTypes();
00534     MEDSKYLINEARRAY *        number = support->getnumber();
00535     for ( int t = 0; t < nb_types; ++t )
00536     {
00537       int shift = getElemNbShift( support->getEntity(), types[t], RSLT_ADD, /*prev=*/true);
00538       if ( shift == 0 ) continue;
00539       int nb_elems = support->getNumberOfElements( types[t] );
00540       for ( int j = 0; j < nb_elems; ++j )
00541         number->setIJ( t+1, j+1, number->getIJ( t+1, j+1 ) + shift );
00542     }
00543   }
00544   return support;
00545 }
00546 
00547 
00548 //================================================================================
00555 //================================================================================
00556 
00557 template< class TSUPPORT >
00558 TSUPPORT* MeshFuse::makeSupport(const TSUPPORT* add_support,
00559                                 TSUPPORT*       same_name_support)
00560 {
00561   if ( same_name_support && same_name_support->getEntity() != add_support->getEntity() )
00562     throw MEDEXCEPTION("MeshFuse: Supports with equal names and different entity!");
00563 
00564   // make resulting support
00565 
00566   TSUPPORT* res_support = same_name_support;
00567   if ( !same_name_support )
00568   {
00569     res_support = new TSUPPORT;
00570     res_support->setName       ( add_support->getName() );
00571     res_support->setDescription( add_support->getDescription() );
00572     res_support->setEntity     ( add_support->getEntity() );
00573     res_support->setMesh       ( this );
00574   }
00575   else if ( same_name_support->isOnAllElements() && add_support->isOnAllElements() )
00576   {
00577     same_name_support->update();
00578     return same_name_support;
00579   }
00580 
00581   // count nb of types
00582 
00583   set<medGeometryElement> add_types, old_types, all_types;
00584   if ( res_support->getEntity() == MED_NODE )
00585     add_types.insert( _NODE_TYPE_ );
00586   else
00587     add_types.insert( add_support->getTypes(),
00588                       add_support->getTypes() + add_support->getNumberOfTypes() );
00589   all_types = add_types;
00590   if ( same_name_support )
00591   {
00592     if ( same_name_support->getEntity() == MED_NODE )
00593       old_types.insert( _NODE_TYPE_ );
00594     else
00595       old_types.insert(same_name_support->getTypes(),
00596                        same_name_support->getTypes()+same_name_support->getNumberOfTypes() );
00597     all_types.insert( old_types.begin(), old_types.end() );
00598   }
00599 
00600   // make all data
00601 
00602   int nb_geom_types = all_types.size();
00603   vector< medGeometryElement > geom_types( nb_geom_types );
00604   vector< vector< int > >        elements( nb_geom_types );
00605   vector< int >               nb_elements( nb_geom_types );
00606   vector< int >                     index( 1,1 );
00607   set<medGeometryElement>::iterator type = all_types.begin();
00608   for ( int t = 0; type != all_types.end(); ++type, ++t )
00609   {
00610     uniteSupportElements( add_types.count( *type ) ? add_support       : 0,
00611                           old_types.count( *type ) ? same_name_support : 0,
00612                           *type,
00613                           elements[ t ]);
00614     nb_elements[ t ] = elements[ t ].size();
00615     index.push_back( index.back() + nb_elements[ t ] );
00616   }
00617 
00618   // set data to support
00619 
00620   res_support->setAll( false );
00621   res_support->setNumberOfGeometricType( nb_geom_types );
00622   res_support->setGeometricType( &geom_types[0] );
00623   res_support->setNumberOfElements( &nb_elements[0] );
00624   res_support->setNumber (new MEDSKYLINEARRAY( nb_geom_types, index[nb_geom_types]-1));
00625   res_support->getnumber()->setIndex( & index[0] );
00626   for ( int t = 0; t < nb_geom_types; ++t )
00627     res_support->getnumber()->setI( t+1, & elements[t][0] );
00628 
00629   return res_support;
00630 }
00631 
00632 //================================================================================
00636 //================================================================================
00637 
00638 void MeshFuse::expandSupports()
00639 {
00640   // we unite supports hiving same names
00641 
00642   // make maps of updated old supports
00643   map< string, FAMILY* > name_to_family;
00644   map< string, GROUP* >  name_to_group;
00645   set<int> family_ids;
00646 
00647   vector<FAMILY*>* fams[4] = { &_familyNode, &_familyCell, &_familyFace, &_familyEdge};
00648   for ( int i = 0; i < 4; ++i )
00649     for ( int f = 0; f < fams[i]->size(); ++f )
00650     {
00651       name_to_family.insert
00652         ( make_pair( fams[i]->at(f)->getName(), updateOldSupport( fams[i]->at(f) )));
00653       family_ids.insert( fams[i]->at(f)->getIdentifier() );
00654     }
00655 
00656   vector<GROUP*>* groups[4] = { &_groupNode, &_groupCell, &_groupFace, &_groupEdge };
00657   for ( int i = 0; i < 4; ++i )
00658     for ( int g = 0; g < groups[i]->size(); ++g )
00659     {
00660       name_to_group.insert
00661         ( make_pair( groups[i]->at(g)->getName(), updateOldSupport( groups[i]->at(g) )));
00662     }
00663 
00664   // unite supports
00665   vector<FAMILY*> add_fams[4]={ _mesh->getFamilies(MED_NODE),
00666                                 _mesh->getFamilies(MED_CELL),
00667                                 _mesh->getFamilies(MED_FACE),
00668                                 _mesh->getFamilies(MED_EDGE) };
00669   for ( int i = 0; i < 4; ++i )
00670     for ( int f = 0; f < add_fams[i].size(); ++f )
00671     {
00672       FAMILY* add_fam = add_fams[i][f];
00673       map<string,FAMILY*>::iterator name_fam = name_to_family.find( add_fam->getName());
00674       FAMILY* res_fam = makeSupport( add_fam,
00675                                      name_fam == name_to_family.end() ? 0 : name_fam->second);
00676       if ( name_to_family.insert( make_pair( res_fam->getName(), res_fam )).second )
00677       {
00678         fams[ i ]->push_back( res_fam );
00679         int id = add_fam->getIdentifier();
00680         if ( family_ids.count ( id ))
00681           id = ( id < 0 ) ? *family_ids.begin() - 1 : *family_ids.rbegin() + 1;
00682         res_fam->setIdentifier( id );
00683         family_ids.insert( id );
00684       }
00685       // update group names in a family
00686       vector<string> res_gr_names( res_fam->getGroupsNames(),
00687                                    res_fam->getGroupsNames()+res_fam->getNumberOfGroups());
00688       for ( int g = 0; g < add_fam->getNumberOfGroups(); ++g )
00689         if ( find( res_gr_names.begin(), res_gr_names.end(), add_fam->getGroupName(g+1))
00690              == res_gr_names.end())
00691           res_gr_names.push_back( add_fam->getGroupName( g+1 ));
00692       if ( res_fam->getNumberOfGroups() < res_gr_names.size() )
00693       {
00694         res_fam->setNumberOfGroups( res_gr_names.size() );
00695         res_fam->setGroupsNames( &res_gr_names[0] );
00696       }
00697     }
00698 
00699   vector<GROUP*> add_grp[4] ={ _mesh->getGroups(MED_NODE),
00700                                _mesh->getGroups(MED_CELL),
00701                                _mesh->getGroups(MED_FACE),
00702                                _mesh->getGroups(MED_EDGE) };
00703   for ( int i = 0; i < 4; ++i )
00704   {
00705     for ( int g = 0; g < add_grp[i].size(); ++g )
00706     {
00707       map< string, GROUP* >::iterator name_grp = name_to_group.find( add_grp[i][g]->getName());
00708 //       STRING out("OLD GRP: ");
00709 //       if ( name_grp == name_to_group.end() ) out << "none";
00710 //       else out << *(name_grp->second);
00711 //       out._s << endl;
00712 //       out << "ADD GRP: " << *(add_grp[i][g]); out._s << endl; out << "";
00713       GROUP* res_grp = makeSupport( add_grp[i][g],
00714                                     name_grp == name_to_group.end() ? 0 : name_grp->second);
00715       if ( !name_to_group.count( res_grp->getName() ))
00716         groups[ i ]->push_back( res_grp );
00717       if ( res_grp->getFamilies().size() < add_grp[i][g]->getFamilies().size() )
00718       {
00719         res_grp->setFamilies( add_grp[i][g]->getFamilies() );
00720         res_grp->setNumberOfFamilies( add_grp[i][g]->getNumberOfFamilies() );
00721       }
00722 //       out << "RES GRP: " << *res_grp;
00723 //       cout << out << endl;
00724     }
00725     // update pointers to families in groups
00726     for ( int g = 0; g < groups[i]->size(); ++g )
00727     {
00728       GROUP* grp = groups[i]->at(g);
00729       if ( grp->getNumberOfFamilies() < 1 ) continue;
00730       vector<FAMILY*> fams = grp->getFamilies();
00731       for ( int f = 0; f < fams.size(); ++f )
00732         fams[f] = name_to_family[ fams[f]->getName() ];
00733       grp->setFamilies( fams );
00734       grp->setNumberOfFamilies( fams.size() );
00735     }
00736   }
00737 }
00738 
00739 //================================================================================
00745 //================================================================================
00746 
00747 int MeshFuse::getElemNbShift( const medEntityMesh& entity,
00748                               medGeometryElement   type,
00749                               int which, bool prev ) const
00750 {
00751   //if ( type == MED_NONE ) type = MED_POINT1;
00752 
00753   const TNbOfGeom & shift_of_type = _nb_index[ which ][ int(entity) ];
00754 
00755   TNbOfGeom::const_iterator type_shift = shift_of_type.lower_bound( type );
00756   if ( type_shift == shift_of_type.end() )
00757     return shift_of_type.empty() ? 0 : shift_of_type.rbegin()->second;
00758 
00759   // get shift of exactly the type or of the previos type
00760   if ( ( prev && type_shift->first == type ) || type_shift->first > type )
00761   {
00762     if ( type_shift == shift_of_type.begin() )
00763       return 0;
00764     else
00765       type_shift--;
00766   }
00767 
00768   return type_shift->second;
00769 }
00770 
00771 //================================================================================
00779 //================================================================================
00780 
00781 void MeshFuse::uniteSupportElements(const SUPPORT*     add_support,
00782                                     SUPPORT*           old_support,
00783                                     medGeometryElement type,
00784                                     vector<int> &      elements)
00785 {
00786   int sup_type = ( type/100 == 0 ? MED_ALL_ELEMENTS : type );
00787 
00788   const medEntityMesh entity = (add_support ? add_support : old_support )->getEntity();
00789 
00790   // Add old elements
00791 
00792   int nb_add_elems = add_support ? add_support->getNumberOfElements( sup_type ) : 0;
00793   int nb_old_elems = 0;
00794   if ( old_support )
00795   {
00796     nb_old_elems = old_support->getNumberOfElements( sup_type );
00797     elements.reserve( nb_old_elems + nb_add_elems );
00798     int add_shift = getElemNbShift( entity, type, RSLT_ADD, /*prev=*/1);
00799     int old_shift = getElemNbShift( entity, type, INIT_OLD, /*prev=*/1);
00800     int shift = 1 + add_shift + old_shift;
00801     // convertion of numbers is already done in updateOldSupport()
00802     if ( old_support->isOnAllElements() )
00803       for ( int i = 0; i < nb_old_elems; ++i )
00804         elements.push_back( i + shift );
00805     else
00806       elements.insert( elements.end(), old_support->getNumber( sup_type ),
00807                        old_support->getNumber( sup_type ) + nb_old_elems );
00808     if ( nb_add_elems == 0 )
00809       return;
00810   }
00811   else
00812   {
00813     elements.reserve( nb_add_elems );
00814   }
00815 
00816   // Convert and add elements of the added support
00817 
00818   const int * add_elems = add_support->isOnAllElements() ? 0 : add_support->getNumber( sup_type );
00819 
00820   int add_shift = getElemNbShift( entity, type, RSLT_ADD, /*prev=*/1);
00821   int old_shift = getElemNbShift( entity, type, INIT_OLD, /*prev=*/0);
00822 
00823   if ( _merged_of_type[ type ].empty() )
00824   {
00825     // no coicident elements in the old and the added meshes - just unite two lists
00826     int shift = add_support->isOnAllElements() ? 1 + add_shift + old_shift : old_shift;
00827     if ( add_support->isOnAllElements() )
00828       for ( int i = 0; i < nb_add_elems; ++i )
00829         elements.push_back( i + shift );
00830     else
00831       for ( int i = 0; i < nb_add_elems; ++i )
00832         elements.push_back( add_elems[i] + shift );
00833   }
00834   else
00835   {
00836     // some elements of the added mesh coincide with old ones,
00837     // so conversion is not so easy, and some support members can
00838     // be twice in it - skip them
00839     vector<int> & new_elem_ids = _new_elem_ids_of_type[ type ];
00840     if ( new_elem_ids.empty() )
00841       makeNewElemIds( entity, type, new_elem_ids );
00842 
00843     set< int > old_elems_in_support( elements.begin(), elements.end() );
00844     int last_old_nb = old_shift + add_shift;
00845     int shift = 1 + getElemNbShift( entity, type, INIT_ADD, /*prev=*/1);
00846     for ( int i = 0; i < nb_add_elems; ++i )
00847     {
00848       int new_id = new_elem_ids[ add_elems ? add_elems[i]-shift : i ];
00849       if ( new_id > last_old_nb || old_elems_in_support.count( new_id ) == 0 )
00850         elements.push_back( new_id );
00851     }
00852   }
00853 }
00854 
00855 //================================================================================
00861 //================================================================================
00862 
00863 void MeshFuse::makeNewElemIds(medEntityMesh      entity,
00864                               medGeometryElement type,
00865                               vector< int > &    new_ids)
00866 {
00867   const char* LOC = "MeshFuse::makeNewElemIds(): ";
00868   if ( entity == MED_NODE )
00869     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"we must not be here"));
00870 
00871   vector<int>& merged_i = _merged_of_type[ type ];
00872   vector<int>::iterator rm_i = merged_i.begin();
00873 
00874   // find ids for merged added elements
00875   vector< int >& old_ids = _equalo_of_type[ type ];
00876 //   if ( old_ids.empty() && !merged_i.empty() )
00877 //     findEqualOldElements( entity, type, old_ids );
00878   vector< int >::iterator old_id = old_ids.begin();
00879 
00880   // nb of added elements
00881   int add_nb_elems = _mesh->getNumberOfElements( entity, type );
00882   new_ids.resize( add_nb_elems );
00883 
00884   // new id of the first added element
00885   int old_shift = getElemNbShift( entity, type, INIT_OLD, /*prev=*/0);
00886   int add_shift = getElemNbShift( entity, type, RSLT_ADD, /*prev=*/1);
00887   int elem_id = old_shift + add_shift + 1;
00888 
00889   // all new ids
00890   // (this works implying that numbers in SUPPORT are in inceasing order!)
00891   for (int i_elem = 0; i_elem < add_nb_elems; )
00892   {
00893     int i_until = ( rm_i == merged_i.end() ? add_nb_elems : *rm_i++ );
00894     // increase elem_id until the next merged element
00895     while ( i_elem < i_until )
00896       new_ids[ i_elem++ ] = elem_id++;
00897     // insert id of old element equal to merged one, if any
00898     if ( i_elem < add_nb_elems )
00899     {
00900       if ( *old_id )
00901         new_ids[ i_elem++ ] = *old_id++;
00902       else
00903         new_ids[ i_elem++ ] = elem_id++, old_id++; // no equal old elem exist
00904     }
00905   }
00906 }
00907 
00908 //================================================================================
00912 //================================================================================
00913 
00914 void MeshFuse::findEqualOldElements(medEntityMesh      entity,
00915                                     medGeometryElement type,
00916                                     vector< int > &    old_ids)
00917 {
00918   // poly element can coincide with any type of the same entity
00919   const bool isPoly = ( type == MED_POLYGON || type == MED_POLYHEDRA );
00920   const medGeometryElement checkType = isPoly ? MED_ALL_ELEMENTS : type;
00921 
00922   if ( !_mesh->getNumberOfElements(entity, type) ||
00923        ! this->getNumberOfElements(entity, checkType) )
00924     return;
00925 
00926   int old_nb_elems_end, old_nb_elems_start;
00927   if ( isPoly )
00928     {
00929       old_nb_elems_start = 0;
00930       old_nb_elems_end   = this->getNumberOfElements( entity, MED_ALL_ELEMENTS );
00931     }
00932   else
00933     {
00934       // if this method is called when connectivity of the entity is not yet concatenated:
00935       old_nb_elems_start = getElemNbShift( entity, type, INIT_OLD, /*prev=*/true);
00936       old_nb_elems_end   = getElemNbShift( entity, type, INIT_OLD, /*prev=*/false);
00937       // if this method is called after expandConnectivity() and this mesh contains all elements
00938       //   int old_nb_elems = 
00939     }
00940   const int *old_conn, *old_index, *add_conn, *add_index;
00941   add_conn  = _mesh->getConnectivity( MED_NODAL, entity, type );
00942   old_conn  =  this->getConnectivity( MED_NODAL, entity, checkType );
00943   add_index = _mesh->getConnectivityIndex( MED_NODAL, entity );
00944   old_index =  this->getConnectivityIndex( MED_NODAL, entity );
00945 
00946   const vector<int>& new_node_ids = _new_elem_ids_of_type[ _NODE_TYPE_ ];
00947 
00948   const vector<int>& merged_i = _merged_of_type[ type ];
00949   vector<int>::const_iterator rm_i = merged_i.begin();
00950 
00951   old_ids.reserve( merged_i.size() );
00952 
00953   for ( ; rm_i != merged_i.end(); ++rm_i ) // rm_i counts from 0
00954   {
00955     // get nodes of rm_i-th added face
00956     const int* add_elem_conn = add_conn + add_index[ *rm_i ]-1;
00957     int    nb_add_elem_nodes = add_index[ *rm_i+1 ] - add_index[ *rm_i ];
00958     set<int> add_elem_nodes;
00959     for ( int n = 0; n < nb_add_elem_nodes; ++n )
00960       add_elem_nodes.insert( new_node_ids[ add_elem_conn[n]-1 ]);
00961 
00962     // look for equal old face among all old faces
00963     const int* old_node = old_conn;
00964     int        new_node = *add_elem_nodes.begin();
00965     int  found_old_elem = 0;
00966     for ( int i = old_nb_elems_start; i < old_nb_elems_end && !found_old_elem; ++i )
00967     {
00968       int nb_old_elem_nodes = old_index[ i+1 ] - old_index[ i ];
00969       for ( int j = 0; j < nb_old_elem_nodes; ++j, ++old_node )
00970       {
00971         if ( new_node != *old_node ) continue;
00972         // found common node, compare all nodes
00973         const int* old_elem_conn = old_conn + old_index[ i ]-1;
00974         set<int> old_elem_nodes( old_elem_conn, old_elem_conn + nb_old_elem_nodes);
00975         if ( add_elem_nodes == old_elem_nodes )
00976         {
00977           found_old_elem = i + 1;
00978           break;
00979         }
00980       }
00981     }
00982     // Issue 0020809. Its possible that no old element exists
00983 //     if ( !found_old_elem )
00984 //       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"all nodes of added elements are merged "
00985 //                                    "but equal old element not found. "
00986 //                                    "Please check global nodes numbers."));
00987     old_ids.push_back( found_old_elem );
00988   }
00989 }
00990 
00991 //================================================================================
00995 //================================================================================
00996 
00997 void MeshFuse::append( medEntityMesh      entity,
00998                        vector<int>&       numbers,
00999                        const vector<int>& add_numbers )
01000 {
01001   const char* LOC="MeshFuse::append(): ";
01002 
01003   int nb_types = getNumberOfTypes( entity );
01004   if ( numbers.empty() || add_numbers.empty() ||
01005        ( nb_types < 2 && _merged_of_type[ getElementType( entity, 1 )].empty() ))
01006   {
01007     numbers.insert( numbers.end(), add_numbers.begin(), add_numbers.end() );
01008     return;
01009   }
01010 
01011   vector<int> result;
01012   result.reserve( numbers.size() + add_numbers.size() );
01013   const int* old_nums = & numbers[0];
01014   const int* add_nums = & add_numbers[0];
01015 
01016   MEDMEM::PointerOf<medGeometryElement> types;
01017   types.setShallowAndOwnership( getTypes(entity));
01018 
01019   for ( int t = 0; t < nb_types; ++t )
01020   {
01021     int nb_old =
01022       getElemNbShift( entity, types[t], INIT_OLD, /*prev=*/false) -
01023       getElemNbShift( entity, types[t], INIT_OLD, /*prev=*/true);
01024     int nb_add =
01025       getElemNbShift( entity, types[t], INIT_ADD, /*prev=*/false) -
01026       getElemNbShift( entity, types[t], INIT_ADD, /*prev=*/true);
01027 
01028     result.insert( result.end(), old_nums, old_nums + nb_old );
01029     old_nums += nb_old;
01030 
01031     const vector<int>& loc_merged = _merged_of_type[ types[t] ];
01032     if ( loc_merged.empty() )
01033     {
01034       result.insert( result.end(), add_nums, add_nums + nb_add );
01035     }
01036     else
01037     {
01038       vector<int>::const_iterator imerged = loc_merged.begin();
01039       vector<int>::const_iterator equalid = _equalo_of_type[ types[t] ].begin();
01040 
01041       int from = 0, to = -1;
01042       while ( from < nb_add )
01043       {
01044         while ( imerged != loc_merged.end() && *equalid == 0 )
01045           imerged++, equalid++;
01046         to = ( imerged == loc_merged.end() ? nb_add : ( equalid++, *imerged++ ));
01047         if ( to > from )
01048           result.insert( result.end(), add_nums + from, add_nums + to );
01049         from = to + 1;
01050       }
01051     }
01052     add_nums += nb_add;
01053   }
01054   if ( result.size() != getNumberOfElements( entity, MED_ALL_ELEMENTS ))
01055     throw MED_EXCEPTION(MEDMEM::STRING(LOC) << "invalid nb of numbers of entity " << entity
01056                         << ": expect " << getNumberOfElements( entity, MED_ALL_ELEMENTS)
01057                         << " but get " << result.size());
01058 
01059   numbers.swap(result);
01060 }