1 // Copyright (C) 2007-2012 CEA/DEN, EDF R&D, OPEN CASCADE
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.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 // File : MEDMEM_MeshFuse.cxx
20 // Created : Tue Jul 7 18:27:00 2009
21 // Author : Edward AGAPOV (eap)
23 #include "MEDMEM_MeshFuse.hxx"
25 #include "MEDMEM_Family.hxx"
26 #include "MEDMEM_Group.hxx"
28 using namespace MEDMEM;
29 using namespace MED_EN;
32 #define _NODE_TYPE_ MED_NONE
34 MeshFuse::MeshFuse():MESHING()
42 //================================================================================
44 * \brief set global numbers of nodes if MeshFuse has been filled via MESHING
46 //================================================================================
48 void MeshFuse::setNodeNumbers( const std::vector<int>& node_glob_numbers )
50 const char* LOC = "MeshFuse::setNodeNumbers(node_glob_numbers): ";
52 if ( !_node_glob_numbers.empty() )
53 throw MEDEXCEPTION(STRING(LOC)<<"node numbers has been already set");
55 if ( node_glob_numbers.size() != getNumberOfNodes() &&
56 node_glob_numbers.size() > 0 && getNumberOfNodes() > 0 )
58 (STRING(LOC)<<"size of node_glob_numbers must be equal number of nodes in MeshFuse");
60 _node_glob_numbers = node_glob_numbers;
63 //================================================================================
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
69 //================================================================================
71 void MeshFuse::concatenate( const MESH* mesh, const vector<int>& node_glob_numbers )
73 const char* LOC = "MeshFuse::concatenate( mesh, node_glob_numbers ): ";
74 if ( !mesh || mesh->getNumberOfNodes() == 0) return;
78 if ( this->getNumberOfNodes() < 1 )
80 mesh->getCoordinates( MED_FULL_INTERLACE );// make fields filled in case of GRID
81 mesh->getConnectivityptr();
84 static_cast<MESH&>(*this) = *mesh;
86 _node_glob_numbers = node_glob_numbers;
91 if ( mesh->getNumberOfNodes() > 0 && node_glob_numbers.empty() )
92 throw MEDEXCEPTION(STRING(LOC)<<"merging without node global numbers not implemented yet");
94 if ( mesh->getNumberOfNodes() != node_glob_numbers.size() )
95 throw MEDEXCEPTION(STRING(LOC)<<"invalid number of node global numbers");
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");
102 _merged_of_type.clear();
103 for ( int i = 0; i < NB_INDICES; ++i )
105 _nb_index[ i ].clear();
106 _nb_index[ i ].resize( MED_ALL_ENTITIES );
111 int final_nb_nodes = makeNewNodeIds( node_glob_numbers );
113 expandCoordinates(final_nb_nodes);
115 expandConnectivity(final_nb_nodes);
119 // clear unnecessary data
120 _new_elem_ids_of_type.clear();
123 //================================================================================
125 * \brief Return number of nodes in the expanded mesh
126 * \param add_glob_numbers - global ids of nodes to append
128 //================================================================================
130 int MeshFuse::makeNewNodeIds(const vector<int>& add_glob_numbers)
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() );
137 // extand global node numbers
138 _node_glob_numbers.reserve( _node_glob_numbers.size() + add_glob_numbers.size());
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 ));
144 int next_loc_id = getNumberOfNodes() + 1;
146 for ( int i = 0; i < add_glob_numbers.size(); ++i )
148 map<int,int>::iterator glob_loc =
149 glob_ids.insert( make_pair( add_glob_numbers[i], next_loc_id )).first;
151 new_node_ids[i] = glob_loc->second;
153 if ( new_node_ids[i] == next_loc_id ) // unique global number
156 _node_glob_numbers.push_back( add_glob_numbers[i] );
160 merged.push_back( i );
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();
167 return next_loc_id - 1;
170 //================================================================================
172 * \brief Update coordinates
174 //================================================================================
176 void MeshFuse::expandCoordinates(int final_nb_nodes)
178 const int dim = getSpaceDimension();
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 ));
189 coord += nb_old_coord;
190 const double* add_coord =_mesh->getCoordinates( MED_FULL_INTERLACE );
191 if ( _merged_of_type[ _NODE_TYPE_ ].empty())
193 // no coincident nodes in the two meshes, just add coords
194 memcpy( coord, add_coord, _mesh->getNumberOfNodes() * dim * sizeof( double ));
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 )
203 if ( new_node_ids[n] < first_added_node ) continue; // coincident node
204 memcpy( coord, add_coord + n * dim, dim * sizeof( double ));
208 _coordinate->setCoordinates( &medarray, /*shallowCopy=*/true );
210 _numberOfNodes = final_nb_nodes;
213 //================================================================================
215 * \brief Concatenate connectivity of meshes
217 * Current implementation impies that cells can't coincide in the meshes
219 //================================================================================
221 void MeshFuse::expandConnectivity(int final_nb_nodes)
223 const medConnectivity nodal = MED_NODAL;
225 // fill in _nb_index[ INIT_OLD ]
226 for ( medEntityMesh entity = MED_CELL; entity < MED_ALL_ENTITIES; ++entity )
228 if ( existConnectivity( nodal, entity ))
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];
238 CONNECTIVITY *prev_connectivity = 0, *cell_connectivity = 0;
240 // loop on all entities
241 for ( medEntityMesh entity = MED_CELL; entity < MED_ALL_ENTITIES; ++entity )
243 if ( entity == MED_FACE && getMeshDimension() == 2 )
244 continue; // there can be EDGE
246 if ( !_mesh->existConnectivity( nodal, entity ))
248 // no entity in added mesh
249 if ( existConnectivity( nodal, entity ) && prev_connectivity )
250 prev_connectivity->setConstituent
251 ( new CONNECTIVITY( *getConnectivityptr()->getConstituent( entity )));
254 if ( !existConnectivity( nodal, entity ))
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;
264 cell_connectivity->setConstituent( connectivity );
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];
275 // Unite connectivities of std 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));
282 int sum_old = 0, sum_add = 0;
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 )
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 );
298 int nb_old = appendConnectivity( conn_of_type[t], this, entity, *type );
300 _nb_index[ INIT_OLD ][ entity ][ *type ] = sum_old;
302 int nb_add = appendConnectivity( conn_of_type[t],_mesh, entity, *type );
304 _nb_index[ RSLT_ADD ][ entity ][ *type ] = sum_add;
306 count.push_back( count.back() + conn_of_type[t]._nb_elems );
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],
317 &conn_of_type[t]._index[0] );
319 // store connectivity of an entity
320 if ( !prev_connectivity )
322 cell_connectivity = connectivity;
323 prev_connectivity = cell_connectivity;
327 prev_connectivity->setConstituent( connectivity );
328 prev_connectivity = connectivity;
331 } // loop on entities
333 if ( cell_connectivity )
335 delete _connectivity;
336 _connectivity = cell_connectivity;
340 //================================================================================
342 * \brief Update node ids in the copied connectivity of theadded mesh
344 //================================================================================
346 void MeshFuse::updateNodeIds( CONNECTIVITY* connectivity )
348 const medConnectivity nodal = MED_NODAL;
349 const medGeometryElement type = MED_ALL_ELEMENTS;
351 const vector<int>& new_node_ids = _new_elem_ids_of_type[ _NODE_TYPE_ ];
353 medEntityMesh entity = connectivity->getEntity();
355 for ( ; entity < MED_ALL_ENTITIES; entity++ )
357 // Collect all connectivities with their lengths
359 list< pair< const int*, int > > conn_len_list;
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)));
368 list< pair< const int*, int > >::iterator conn_len = conn_len_list.begin();
369 for ( ; conn_len != conn_len_list.end(); ++conn_len )
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 ];
379 //================================================================================
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
388 //================================================================================
390 int MeshFuse::appendConnectivity( MeshFuse::TConnData& data,
392 medEntityMesh entity,
393 medGeometryElement type)
395 // get connectivity of type
397 const int* conn, *index = 0;
398 int nb_elem, conn_len;
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);
406 conn_len = index[ nb_elem ] - index[ 0 ];
408 bool need_index = ( type == MED_POLYGON || type == MED_POLYHEDRA );
410 data._index.resize( 1, 0 ); // for safe access to pointer even if no real index exists
414 int old_size = data._nb_elems;
415 data._nb_elems += nb_elem;
419 // append connectivity to data as is
420 data._connectivity.insert( data._connectivity.end(), conn, conn + conn_len );
423 if ( data._index.empty() )
424 data._index.insert( data._index.end(), index, index + nb_elem + 1 );
426 for ( int i = 0; i < nb_elem; ++i )
427 data._index.push_back( data._index.back() + index[i+1] - index[i] );
432 // convert connectivity of other mesh
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 )
438 for ( int n = 0; n < conn_len; ++n )
440 other_conn[ n ] = new_node_ids[ conn[ n ]-1 ];
444 for ( int n = 0; n < conn_len; ++n )
445 other_conn[ n ] = new_node_ids[ conn[ n ]-1 ];
447 if ( entity == MED_CELL || _merged_of_type[ _NODE_TYPE_ ].empty() )
449 // store converted connectivity in data
450 data._connectivity.insert( data._connectivity.end(), other_conn.begin(), other_conn.end());
453 if ( data._index.empty() && index[0] == 1 )
454 data._index.insert( data._index.end(), index, index + nb_elem + 1 );
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] );
465 // exclude coincident elements from connectivity
467 if ( need_index && data._index.empty() )
468 data._index.push_back( 1 );
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 )
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 );
480 if ( nb_coincident_nodes == index[i+1] - index[i] )
481 merged.push_back( i );
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 );
490 int rm_i = 0, nb_rm = 0;
491 for ( int i = 0; i < nb_elem; ++i )
493 bool is_merged = ( rm_i < merged.size() && i == merged[rm_i] && equalo[rm_i] );
502 for ( int n = index[i]; n < index[i+1]; ++n )
503 data._connectivity.push_back( other_conn[ n-1 ] );
505 data._index.push_back( data._index.back() + index[i+1] - index[i] );
509 merged.clear(), equalo.clear();
512 return data._nb_elems - old_size;
515 //================================================================================
517 * \brief Return updated old support
519 //================================================================================
521 template< class TSUPPORT >
522 TSUPPORT* MeshFuse::updateOldSupport(TSUPPORT* support) const
524 if ( support->isOnAllElements() )
526 //support->update(); -- this changes old nb elems to nb elems after fuse
528 else if ( support->getEntity() != MED_NODE )
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 )
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 );
548 //================================================================================
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
555 //================================================================================
557 template< class TSUPPORT >
558 TSUPPORT* MeshFuse::makeSupport(const TSUPPORT* add_support,
559 TSUPPORT* same_name_support)
561 if ( same_name_support && same_name_support->getEntity() != add_support->getEntity() )
562 throw MEDEXCEPTION("MeshFuse: Supports with equal names and different entity!");
564 // make resulting support
566 TSUPPORT* res_support = same_name_support;
567 if ( !same_name_support )
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 );
575 else if ( same_name_support->isOnAllElements() && add_support->isOnAllElements() )
577 same_name_support->update();
578 return same_name_support;
583 set<medGeometryElement> add_types, old_types, all_types;
584 if ( res_support->getEntity() == MED_NODE )
585 add_types.insert( _NODE_TYPE_ );
587 add_types.insert( add_support->getTypes(),
588 add_support->getTypes() + add_support->getNumberOfTypes() );
589 all_types = add_types;
590 if ( same_name_support )
592 if ( same_name_support->getEntity() == MED_NODE )
593 old_types.insert( _NODE_TYPE_ );
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() );
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 )
610 uniteSupportElements( add_types.count( *type ) ? add_support : 0,
611 old_types.count( *type ) ? same_name_support : 0,
614 nb_elements[ t ] = elements[ t ].size();
615 index.push_back( index.back() + nb_elements[ t ] );
618 // set data to support
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] );
632 //================================================================================
634 * \brief Concatenate families and groups
636 //================================================================================
638 void MeshFuse::expandSupports()
640 // we unite supports hiving same names
642 // make maps of updated old supports
643 map< string, FAMILY* > name_to_family;
644 map< string, GROUP* > name_to_group;
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 )
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() );
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 )
661 ( make_pair( groups[i]->at(g)->getName(), updateOldSupport( groups[i]->at(g) )));
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 )
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 )
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 );
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() )
694 res_fam->setNumberOfGroups( res_gr_names.size() );
695 res_fam->setGroupsNames( &res_gr_names[0] );
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 )
705 for ( int g = 0; g < add_grp[i].size(); ++g )
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);
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() )
719 res_grp->setFamilies( add_grp[i][g]->getFamilies() );
720 res_grp->setNumberOfFamilies( add_grp[i][g]->getNumberOfFamilies() );
722 // out << "RES GRP: " << *res_grp;
723 // cout << out << endl;
725 // update pointers to families in groups
726 for ( int g = 0; g < groups[i]->size(); ++g )
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() );
739 //================================================================================
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"
745 //================================================================================
747 int MeshFuse::getElemNbShift( const medEntityMesh& entity,
748 medGeometryElement type,
749 int which, bool prev ) const
751 //if ( type == MED_NONE ) type = MED_POINT1;
753 const TNbOfGeom & shift_of_type = _nb_index[ which ][ int(entity) ];
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;
759 // get shift of exactly the type or of the previos type
760 if ( ( prev && type_shift->first == type ) || type_shift->first > type )
762 if ( type_shift == shift_of_type.begin() )
768 return type_shift->second;
771 //================================================================================
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
779 //================================================================================
781 void MeshFuse::uniteSupportElements(const SUPPORT* add_support,
782 SUPPORT* old_support,
783 medGeometryElement type,
784 vector<int> & elements)
786 int sup_type = ( type/100 == 0 ? MED_ALL_ELEMENTS : type );
788 const medEntityMesh entity = (add_support ? add_support : old_support )->getEntity();
792 int nb_add_elems = add_support ? add_support->getNumberOfElements( sup_type ) : 0;
793 int nb_old_elems = 0;
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 );
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 )
813 elements.reserve( nb_add_elems );
816 // Convert and add elements of the added support
818 const int * add_elems = add_support->isOnAllElements() ? 0 : add_support->getNumber( sup_type );
820 int add_shift = getElemNbShift( entity, type, RSLT_ADD, /*prev=*/1);
821 int old_shift = getElemNbShift( entity, type, INIT_OLD, /*prev=*/0);
823 if ( _merged_of_type[ type ].empty() )
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 );
831 for ( int i = 0; i < nb_add_elems; ++i )
832 elements.push_back( add_elems[i] + shift );
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 );
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 )
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 );
855 //================================================================================
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
861 //================================================================================
863 void MeshFuse::makeNewElemIds(medEntityMesh entity,
864 medGeometryElement type,
865 vector< int > & new_ids)
867 const char* LOC = "MeshFuse::makeNewElemIds(): ";
868 if ( entity == MED_NODE )
869 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"we must not be here"));
871 vector<int>& merged_i = _merged_of_type[ type ];
872 vector<int>::iterator rm_i = merged_i.begin();
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();
880 // nb of added elements
881 int add_nb_elems = _mesh->getNumberOfElements( entity, type );
882 new_ids.resize( add_nb_elems );
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;
890 // (this works implying that numbers in SUPPORT are in inceasing order!)
891 for (int i_elem = 0; i_elem < add_nb_elems; )
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 )
901 new_ids[ i_elem++ ] = *old_id++;
903 new_ids[ i_elem++ ] = elem_id++, old_id++; // no equal old elem exist
908 //================================================================================
910 * \brief Finds ids of elements of the old mesh equal to merged elements of the added one
912 //================================================================================
914 void MeshFuse::findEqualOldElements(medEntityMesh entity,
915 medGeometryElement type,
916 vector< int > & old_ids)
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;
922 if ( !_mesh->getNumberOfElements(entity, type) ||
923 ! this->getNumberOfElements(entity, checkType) )
926 int old_nb_elems_end, old_nb_elems_start;
929 old_nb_elems_start = 0;
930 old_nb_elems_end = this->getNumberOfElements( entity, MED_ALL_ELEMENTS );
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 =
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 );
946 const vector<int>& new_node_ids = _new_elem_ids_of_type[ _NODE_TYPE_ ];
948 const vector<int>& merged_i = _merged_of_type[ type ];
949 vector<int>::const_iterator rm_i = merged_i.begin();
951 old_ids.reserve( merged_i.size() );
953 for ( ; rm_i != merged_i.end(); ++rm_i ) // rm_i counts from 0
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 ]);
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 )
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 )
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 )
977 found_old_elem = i + 1;
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 );
991 //================================================================================
993 * \brief Unite numbers of elements taking into account their types
995 //================================================================================
997 void MeshFuse::append( medEntityMesh entity,
998 vector<int>& numbers,
999 const vector<int>& add_numbers )
1001 const char* LOC="MeshFuse::append(): ";
1003 int nb_types = getNumberOfTypes( entity );
1004 if ( numbers.empty() || add_numbers.empty() ||
1005 ( nb_types < 2 && _merged_of_type[ getElementType( entity, 1 )].empty() ))
1007 numbers.insert( numbers.end(), add_numbers.begin(), add_numbers.end() );
1012 result.reserve( numbers.size() + add_numbers.size() );
1013 const int* old_nums = & numbers[0];
1014 const int* add_nums = & add_numbers[0];
1016 MEDMEM::PointerOf<medGeometryElement> types;
1017 types.setShallowAndOwnership( getTypes(entity));
1019 for ( int t = 0; t < nb_types; ++t )
1022 getElemNbShift( entity, types[t], INIT_OLD, /*prev=*/false) -
1023 getElemNbShift( entity, types[t], INIT_OLD, /*prev=*/true);
1025 getElemNbShift( entity, types[t], INIT_ADD, /*prev=*/false) -
1026 getElemNbShift( entity, types[t], INIT_ADD, /*prev=*/true);
1028 result.insert( result.end(), old_nums, old_nums + nb_old );
1031 const vector<int>& loc_merged = _merged_of_type[ types[t] ];
1032 if ( loc_merged.empty() )
1034 result.insert( result.end(), add_nums, add_nums + nb_add );
1038 vector<int>::const_iterator imerged = loc_merged.begin();
1039 vector<int>::const_iterator equalid = _equalo_of_type[ types[t] ].begin();
1041 int from = 0, to = -1;
1042 while ( from < nb_add )
1044 while ( imerged != loc_merged.end() && *equalid == 0 )
1045 imerged++, equalid++;
1046 to = ( imerged == loc_merged.end() ? nb_add : ( equalid++, *imerged++ ));
1048 result.insert( result.end(), add_nums + from, add_nums + to );
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());
1059 numbers.swap(result);