src/MEDMEM/MEDMEM_DriverTools.cxx

   1 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 // This library is free software; you can redistribute it and/or
   7 // modify it under the terms of the GNU Lesser General Public
   8 // License as published by the Free Software Foundation; either
   9 // version 2.1 of the License.
  10 //
  11 // This library is distributed in the hope that it will be useful,
  12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 // Lesser General Public License for more details.
  15 //
  16 // You should have received a copy of the GNU Lesser General Public
  17 // License along with this library; if not, write to the Free Software
  18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22
  23 #include "MEDMEM_DriverTools.hxx"
  24 #include "MEDMEM_STRING.hxx"
  25 #include "MEDMEM_Exception.hxx"
  26 #include "MEDMEM_Mesh.hxx"
  27 #include "MEDMEM_Group.hxx"
  28 #include "MEDMEM_Field.hxx"
  29
  30 #include <iomanip>
  31 #include <algorithm>
  32
  33 using namespace std;
  34 using namespace MED_EN;
  35
  36 #define DUMP_LINES_LIMIT 20
  37
  38 namespace MEDMEM {
  39
  40 // avoid coping sortedNodeIDs
  41 _maille::_maille(const _maille& ma)
  42   : sommets(ma.sommets), geometricType(ma.geometricType),
  43     reverse(ma.reverse), sortedNodeIDs(0), _ordre(ma._ordre)
  44 {
  45 }
  46
  47 // Cet opérateur permet d'ordonner les mailles dans un set suivant l'ordre requis par MED
  48 bool _maille::operator < (const _maille& ma) const
  49 {
  50   // si le type géométrique differe, la comparaison est basée dessus
  51   // sinon on se base sur une comparaison des numéros de sommets
  52   // we compare _mailles of only same geometricType due to maillageByType usage
  53   size_t l=sommets.size();
  54   if ( dimension() > 3 ) { // poly
  55     size_t l2 = ma.sommets.size();
  56     if ( l != l2 )
  57       return l < l2;
  58   }
  59   const int* v1 = getSortedNodes();
  60   const int* v2 = ma.getSortedNodes();
  61   for ( const int* vEnd = v1 + l; v1 < vEnd; ++v1, ++v2 )
  62     if(*v1 != *v2)
  63       return *v1 < *v2;
  64   return false; // cas d'égalité
  65
  66 //   if(geometricType==ma.geometricType)
  67 //   {
  68 //     // construction de deux vecteur temporaire contenant les numeros de sommets
  69 //     // pour faire le tri et les comparaisons
  70 //     size_t l=sommets.size();
  71 //     if ( dimension() > 3 ) { // poly
  72 //       size_t l2 = ma.sommets.size();
  73 //       if ( l != l2 )
  74 //         return l < l2;
  75 //     }
  76 //     std::vector<int> v1(l);
  77 //     std::vector<int> v2(l);
  78 //     for (unsigned int i=0; i!=l; ++i)
  79 //     {
  80 //       v1[i]=nodeNum(i);
  81 //       v2[i]=ma.nodeNum(i);
  82 //     }
  83 //     std::sort(v1.begin(), v1.end());
  84 //     std::sort(v2.begin(), v2.end());
  85 //     for(std::vector<int>::const_iterator i1=v1.begin(), i2=v2.begin(); i1!=v1.end(); ++i1, ++i2)
  86 //       if(*i1 != *i2)
  87 //         return *i1 < *i2;
  88 //     return false; // cas d'égalité
  89 //   }
  90 //   else
  91 //     return geometricType<ma.geometricType;
  92 }
  93
  94 // creates if needed and return sortedNodeIDs
  95 const int* _maille::getSortedNodes() const
  96 {
  97   if ( !sortedNodeIDs )
  98   {
  99     size_t l=sommets.size();
 100     sortedNodeIDs = new int[ l ];
 101
 102     for (size_t i=0; i!=l; ++i)
 103       sortedNodeIDs[i]=nodeID(i);
 104     std::sort( sortedNodeIDs, sortedNodeIDs + l );
 105   }
 106   return sortedNodeIDs;
 107 }
 108
 109 _link _maille::link(int i) const
 110 {
 111   ASSERT_MED ( i >= 0 && i < (int)sommets.size() );
 112   int i2 = ( i + 1 == (int)sommets.size() ) ? 0 : i + 1;
 113   if ( reverse )
 114     return make_pair( sommets[i2]->first, sommets[i]->first );
 115   else
 116     return make_pair( sommets[i]->first, sommets[i2]->first );
 117 }
 118
 119 // retourne l'entité d'une maille en fonction de la dimension du maillage.
 120 MED_EN::medEntityMesh _maille::getEntity(const int meshDimension) const throw (MEDEXCEPTION)
 121 {
 122   const char * LOC = "_maille::getEntity(const int meshDimension)" ;
 123 //   BEGIN_OF_MED(LOC);
 124
 125   int mailleDimension = this->dimensionWithPoly();
 126   medEntityMesh entity;
 127   if (mailleDimension == meshDimension)
 128     entity = MED_CELL;
 129   else
 130     switch (mailleDimension)
 131       {
 132       case 0 :
 133         entity = MED_NODE;
 134         break;
 135       case 1 :
 136         entity = MED_EDGE;
 137         break;
 138       case 2 :
 139         entity = MED_FACE;
 140         break;
 141       default :
 142         throw MEDEXCEPTION(LOCALIZED(STRING(LOC) << "Impossible de determiner l'entite de la maille."));
 143       }
 144 return entity;
 145
 146 //END_OF_MED(LOC);
 147 }
 148
 149 void _maillageByDimIterator::init(const int dim, const bool convertPoly )
 150 {
 151   myIt = myImed->maillageByType.begin();
 152   myEnd = myImed->maillageByType.end();
 153   myDim = dim;
 154   myConvertPoly = convertPoly;
 155   nbRemovedByType = & myImed->nbRemovedByType;
 156 }
 157
 158 std::ostream& operator << (std::ostream& os, const _maille& ma)
 159 {
 160     os << "maille " << ma.ordre() << " (" << ma.geometricType << ") : < ";
 161     os << ma.nodeNum(0);
 162     for( unsigned i=1; i!=ma.sommets.size(); ++i)
 163         os << ", " << ma.nodeNum( i );
 164     os << " > sortedNodeIDs: ";
 165     if ( ma.sortedNodeIDs ) {
 166       os << "< ";
 167       for( unsigned i=0; i!=ma.sommets.size(); ++i)
 168         os << ( i ? ", " : "" ) << ma.sortedNodeIDs[ i ];
 169       os << " >";
 170     }
 171     else {
 172       os << "NULL";
 173     }
 174     if ( ma.isMerged() )
 175       os << " MERGED ";
 176     return os;
 177 }
 178
 179 std::ostream& operator << (std::ostream& os, const _groupe& gr)
 180 {
 181     os << "--- Groupe " << gr.nom << " --- " << std::endl ;
 182     os << " -> liste des sous-groupes : ";
 183     for( std::vector<int>::const_iterator i=gr.groupes.begin(); i!=gr.groupes.end(); ++i)
 184             os << *i << " ";
 185
 186     os << std::endl << " -> liste des "<< gr.mailles.size() << " mailles : " << std::endl;
 187
 188     _groupe::TMailleIter i1=gr.mailles.begin();
 189     int l;
 190     for(l = 0; l < DUMP_LINES_LIMIT && i1!=gr.mailles.end(); i1++, l++)
 191             os << setw(3) << l+1 << " " << *(*i1) << std::endl;
 192
 193     if ( l == DUMP_LINES_LIMIT )
 194       os << "   ... skip " << gr.mailles.size() - l << " mailles" << endl;
 195
 196     os << " relocMap, size=" << gr.relocMap.size() << endl;
 197     map<unsigned,int>::const_iterator it = gr.relocMap.begin();
 198     for ( l = 0; l < DUMP_LINES_LIMIT && it != gr.relocMap.end(); ++it, ++l )
 199       os << " (" << it->first << "," << it->second << ")";
 200     if ( gr.relocMap.size() > 0 )
 201       os << endl;
 202     return os;
 203 }
 204
 205 std::ostream& operator << (std::ostream& os, const _noeud& no)
 206 {
 207     os << "noeud " << no.number << " : < ";
 208     std::vector<double>::const_iterator i=no.coord.begin();
 209     os << *i++ ;
 210     for( ; i!=no.coord.end(); ++i)
 211         os << ", " << *i;
 212     os << " >";
 213     return os;
 214 }
 215
 216 void MEDMEM::_fieldBase::dump(std::ostream& os) const
 217 {
 218   os << "field " << "<" << _name << ">" << endl <<
 219     "  nb sub: " << _sub.size() << endl <<
 220     "  group index: " << _group_id << endl <<
 221     "  type: " << _type << endl;
 222   os << "  subcomponents:" << endl;
 223   vector< _sub_data >::const_iterator sub_data = _sub.begin();
 224   for ( ; sub_data != _sub.end(); ++sub_data ) {
 225     os << "    group index: " << sub_data->_supp_id <<
 226       ", " << sub_data->nbComponents() << " comp names: ";
 227     for ( int i_comp = 0; i_comp < sub_data->nbComponents(); ++i_comp )
 228       os << " |" << sub_data->_comp_names[ i_comp ] << "|";
 229     os << endl;
 230   }
 231 }
 232
 233 std::ostream& operator << (std::ostream& os, const _fieldBase * f)
 234 {
 235   f->dump( os );
 236   return os;
 237 }
 238
 239 std::ostream& operator << (std::ostream& os, const _intermediateMED& mi)
 240 {
 241   int l;
 242   _maillageByDimIterator maIt( mi );
 243   while ( const set<_maille >* maillage = maIt.nextType() )
 244   {
 245     os << "Set des " << maillage->size()
 246        << " mailles of type " << maIt.type() << ": "<< std::endl;
 247     std::set<_maille>::const_iterator i=maillage->begin();
 248     for( l = 0; l < DUMP_LINES_LIMIT && i!=maillage->end(); ++i, ++l)
 249       os << setw(3) << l+1 << " " <<*i << std::endl;
 250     if ( l == DUMP_LINES_LIMIT )
 251       os << "   ... skip " << maillage->size() - l << " mailles" << endl;
 252   }
 253   os << std::endl << "Vector des " << mi.groupes.size() << " groupes : " << std::endl;
 254   for (unsigned int k=0; k!=mi.groupes.size(); k++)
 255     os << std::setw(3) << k << " " << mi.groupes[k] << std::endl;
 256
 257   os << std::endl << "map des " << mi.points.size() << " noeuds : " << std::endl;
 258   std::map<int,_noeud>::const_iterator j=mi.points.begin();
 259   for( l = 0; l < DUMP_LINES_LIMIT && j!=mi.points.end(); ++j, ++l)
 260     os << j->second << std::endl;
 261   if ( l == DUMP_LINES_LIMIT )
 262     os << "   ... skip " << mi.points.size() - l << " noeuds" << endl;
 263
 264   os << endl << mi.fields.size() << " fields:" << endl;
 265   std::list<_fieldBase* >::const_iterator fIt = mi.fields.begin();
 266   for ( l = 0; fIt != mi.fields.end(); ++fIt, ++l )
 267     os << " - " << l+1 << " " << *fIt << endl;
 268
 269   return os;
 270 }
 271
 272
 273 //=======================================================================
 274 //function : treatGroupes
 275 //purpose  : detect groupes of mixed dimension and erase groupes that
 276 //           won't be converted
 277 //=======================================================================
 278
 279 void _intermediateMED::treatGroupes()
 280 {
 281   const char * LOC = "_intermediateMED::treatGroupes() : ";
 282   BEGIN_OF_MED(LOC);
 283
 284   if ( myGroupsTreated )
 285     return;
 286   myGroupsTreated = true;
 287
 288   // --------------------
 289   // erase useless group
 290   // --------------------
 291
 292   // decrease hierarchical depth of subgroups
 293   vector<int>::iterator j;
 294   for (unsigned int i=0; i!=this->groupes.size(); ++i)
 295   {
 296     _groupe& grp = groupes[i];
 297     //INFOS_MED( i << " " << grp.nom );
 298     j = grp.groupes.begin();
 299     while( j!=grp.groupes.end() ) {
 300       int grpInd = *j-1;
 301       if ( grpInd < 0 || grpInd >= (int)groupes.size() ) {
 302         throw MEDEXCEPTION(LOCALIZED(STRING(LOC) << "Bad subgroup index: " << grpInd <<
 303                                      ", in " << i << " groupe.nom=" << grp.nom));
 304       }
 305       _groupe & sub_grp = groupes[ grpInd ];
 306       if ( !sub_grp.groupes.empty() ) {
 307         MESSAGE_MED("High hierarchical depth of subgroups in group " << i );
 308         *j = sub_grp.groupes[0]; // replace j with its 1st subgroup
 309         // push back the rest subs
 310         for ( int k = 1; k < (int)sub_grp.groupes.size(); ++k )
 311           grp.groupes.push_back( sub_grp.groupes[ k ]);
 312         // vector maybe is reallocated: restart iterator
 313         j = grp.groupes.begin();
 314       }
 315       else
 316         j++;
 317     }
 318     // remove empty sub-groupes
 319     j = grp.groupes.begin();
 320     while ( j!=grp.groupes.end() ) {
 321       if ( groupes[*j-1].empty() ) {
 322         grp.groupes.erase( j );
 323         j = grp.groupes.begin();
 324       }
 325       else
 326         j++;
 327     }
 328   }
 329   // get indices of groups that are field support -
 330   // do not erase them and their subgroups
 331   std::set<int> groups2convert;
 332   std::list< _fieldBase* >::const_iterator fIt = fields.begin();
 333   for ( ; fIt != fields.end(); ++fIt )
 334     (*fIt)->getGroupIds( groups2convert, true );
 335
 336   // keep named groups and their subgroups
 337   for (unsigned int i=0; i!=this->groupes.size(); ++i)
 338   {
 339     _groupe& grp = groupes[i];
 340     if ( grp.empty() || grp.nom.empty() )
 341       continue;
 342     groups2convert.insert( i );
 343     for( j = grp.groupes.begin(); j!=grp.groupes.end(); ++j)
 344       groups2convert.insert( *j-1 );
 345   }
 346   // erase groups that are not in groups2convert
 347   for (unsigned int i=0; i!=this->groupes.size(); ++i)
 348   {
 349     if ( groups2convert.find( i ) == groups2convert.end() ) {
 350       _groupe& grp = groupes[i];
 351       grp.mailles.clear();
 352       grp.groupes.clear();
 353       MESSAGE_MED( "Erase " << i << "-th group " << grp.nom );
 354     }
 355   }
 356
 357   // ---------------------------------------------------
 358   // define if there are groups with mixed entity types
 359   // ---------------------------------------------------
 360
 361   hasMixedCells = false;
 362   for (unsigned int i=0; i!=this->groupes.size(); ++i)
 363   {
 364     _groupe& grp = groupes[i];
 365     if ( grp.groupes.empty() )
 366       continue;
 367
 368     // check if sub-groups have different dimension
 369     j = grp.groupes.begin();
 370     int dim = groupes[*j-1].mailles[0]->dimension();
 371     for( j++; j!=grp.groupes.end(); ++j) {
 372       int dim2 = groupes[*j-1].mailles[0]->dimension();
 373       if ( dim != dim2 ) {
 374         if ( dim == 0 || dim2 == 0 || dim + dim2 == 9 ) {
 375           // cant create a group of nodes plus anything else,
 376           // neither a group of polygones + polyhedrons
 377           grp.mailles.clear();
 378           grp.groupes.clear();
 379           MESSAGE_MED( "Erase mixed dim group with nodes:" << i << "-th group " << grp.nom );
 380           break;
 381         }
 382         int d1 = dim  < 4 ? dim  : dim  - 2; // care of poly
 383         int d2 = dim2 < 4 ? dim2 : dim2 - 2;
 384         if ( d1 != d2 ) {
 385           hasMixedCells = true;
 386           MESSAGE_MED( "Mixed dim group: " << i << "-th " << grp.nom <<
 387                    "  dim1 = " << dim << " dim2 = " << dim2 );
 388         }
 389       }
 390     }
 391   }
 392
 393 //   if ( hasMixedCells )
 394 //     INFOS_MED( "There will be groups of mixed dimention" );
 395   END_OF_MED(LOC);
 396 }
 397
 398 void _intermediateMED::numerotationMaillage()
 399 {
 400   const char* LOC = "_intermediateMED::numerotationMaillage() : ";
 401   BEGIN_OF_MED(LOC);
 402   if ( myMaillesNumerated )
 403     return;
 404   myMaillesNumerated = true;
 405
 406   treatGroupes();
 407
 408   set<_maille>::const_iterator i, iEnd;
 409
 410   // numerotation mailles of type MED_POINT1 by node number
 411   bool hasPointMailles = false;
 412   if ( const set<_maille> * points = _maillageByDimIterator( *this, 0 ).nextType() ) {
 413     hasPointMailles = true;
 414     numerotationPoints();
 415     for ( i = points->begin(), iEnd = points->end(); i != iEnd; ++i )
 416       i->setOrdre( i->sommets[0]->second.number ); // is merged if point is merged
 417   }
 418
 419   // loop on entities
 420   for ( int dim = 1; dim <= 3; ++dim )
 421   {
 422     int iterDim = hasMixedCells ? -1 : dim;
 423     const bool skipFirstType = ( hasPointMailles && hasMixedCells );
 424
 425     // check if any numeration is present
 426     bool hasNumeration = true;
 427     _maillageByDimIterator entityMailles( *this, iterDim, true );
 428     if ( skipFirstType ) entityMailles.nextType();
 429     while ( const set<_maille> * typeMailles = entityMailles.nextType() ) {
 430       if ( typeMailles->begin()->ordre() == 0 || typeMailles->rbegin()->ordre() == 0 ) {
 431         hasNumeration = false;
 432         break;
 433       }
 434     }
 435     // check if re-numeration is needed
 436     bool ok = false, renumEntity = false;
 437     if ( hasNumeration )
 438     {
 439       ok = true;
 440       _maillageByDimIterator entityMailles( *this, iterDim, true );
 441       if ( skipFirstType ) entityMailles.nextType();
 442
 443       int prevNbElems = 0;
 444       while ( const set<_maille> * typeMailles = entityMailles.nextType() )
 445       {
 446         unsigned minOrdre = INT_MAX, maxOrdre = 0;
 447         for ( i = typeMailles->begin(), iEnd = typeMailles->end(); i!=iEnd; ++i) {
 448           if ( i->ordre() < minOrdre ) minOrdre = i->ordre();
 449           if ( i->ordre() > maxOrdre ) maxOrdre = i->ordre();
 450         }
 451         unsigned typeSize = entityMailles.sizeWithoutMerged();
 452         if ( typeSize != maxOrdre - minOrdre + 1 )
 453           ok = false;
 454         if ( prevNbElems != 0 ) {
 455           if ( minOrdre == 1 )
 456             renumEntity = true;
 457           else if ( prevNbElems+1 != (int)minOrdre )
 458             ok = false;
 459         }
 460         prevNbElems += typeSize;
 461       }
 462
 463       if ( ok && renumEntity ) // each geom type was numerated separately
 464       {
 465         entityMailles.init( iterDim, true );
 466         if ( skipFirstType ) entityMailles.nextType();
 467         prevNbElems = entityMailles.nextType()->size(); // no need to renumber the first type
 468         while ( const set<_maille> * typeMailles = entityMailles.nextType() ) {
 469           for ( i = typeMailles->begin(), iEnd = typeMailles->end(); i!=iEnd; ++i)
 470             i->setOrdre( i->ordre() + prevNbElems );
 471           prevNbElems += typeMailles->size();
 472         }
 473       }
 474     }
 475     if ( !ok )
 476     {
 477       int i_maille=0;
 478       entityMailles.init( iterDim, true );
 479       if ( skipFirstType ) entityMailles.nextType();
 480       while ( const set<_maille> * typeMailles = entityMailles.nextType() )
 481         for ( i = typeMailles->begin(), iEnd = typeMailles->end(); i!=iEnd; ++i)
 482           i->setOrdre( ++i_maille );
 483     }
 484   }
 485   END_OF_MED(LOC);
 486 }
 487
 488 bool _intermediateMED::numerotationPoints()
 489 {
 490   if ( !myNodesNumerated ) // is negative if numerated by merge
 491   {
 492     int i_noeud=0;
 493     for( map<int,_noeud>::iterator i=points.begin(); i!=points.end(); ++i)
 494       i->second.number = ++i_noeud ;
 495     myNodesNumerated = true;
 496     return true;
 497   }
 498   return false;
 499 }
 500
 501 int _intermediateMED::nbMerged(int type) const //!< nb nodes removed by merge
 502 {
 503   TNbByType::const_iterator typeNb = nbRemovedByType.find( type );
 504   return ( typeNb == nbRemovedByType.end() ? 0 : typeNb->second );
 505 }
 506
 507
 508 /*!
 509  * \if developper
 510  * create a MED COORDINATE from the intermediate structure.
 511  * \endif
 512  */
 513 COORDINATE * _intermediateMED::getCoordinate(const string & coordinateSystem)
 514 {
 515     const medModeSwitch mode=MED_FULL_INTERLACE;
 516     int spaceDimension=points.begin()->second.coord.size();
 517     int numberOfNodes=points.size() - nbMerged( MED_POINT1 );
 518
 519     // creation du tableau des coordonnees en mode MED_FULL_INTERLACE
 520     double * coord = new double[spaceDimension*numberOfNodes];
 521     double * xyz = coord;
 522     for( std::map<int,_noeud>::const_iterator i=points.begin(); i!=points.end(); ++i )
 523       if ( i->second.number > 0 ) {
 524         std::copy(i->second.coord.begin(), i->second.coord.end(), xyz );
 525         xyz += spaceDimension;
 526       }
 527
 528     // creation de l'objet COORDINATE
 529     COORDINATE * coordinate = new COORDINATE(spaceDimension, numberOfNodes, mode);
 530     coordinate->setCoordinates(mode,coord);
 531     delete [] coord;
 532     coordinate->setCoordinatesSystem(coordinateSystem);
 533     return coordinate;
 534 }
 535
 536
 537 /*!
 538  * \if developper
 539  * create a MED CONNECTIVITY from the intermediate structure.
 540  * \endif
 541  */
 542 CONNECTIVITY * _intermediateMED::getConnectivity()
 543 {
 544   const char * LOC = "_intermediateMED::getConnectivity() : ";
 545   BEGIN_OF_MED(LOC);
 546
 547   int numberOfNodes=points.size() - nbMerged( MED_POINT1 ); // number of nodes in the mesh
 548   medEntityMesh entity;
 549   CONNECTIVITY *Connectivity = NULL, *Constituent = NULL;
 550
 551   set<_maille>::const_iterator i, iEnd; // iterateurs sur les mailles
 552
 553   // find out mesh dimension
 554   medGeometryElement meshDim;
 555   _maillageByDimIterator allMailles( *this, -1, true );
 556   while ( allMailles.nextType() )
 557     meshDim = allMailles.dim();
 558
 559   // renumerote les points de 1 a n (pour le cas ou certains points ne sont pas presents dans le maillage d'origine)
 560   numerotationPoints();
 561
 562   // loop on entities
 563   for ( int dim = 0; dim <= 3; ++dim )
 564   {
 565     // skip nodes and elements of <dimension_maillage - 2> or less dimension
 566     // Unfortunately, it is impossible because of MESH::createFamilies() that requires
 567     // presence of connectivity even for nodes!
 568     //int dimension_maillage_moin_2=maillage.rbegin()->dimension() - 2;
 569
 570     // tableau de travail : nombre d'elements pour chaque type geometrique
 571     vector<int> count;
 572     count.reserve( maillageByType.size() );
 573     count.push_back( 1 );
 574
 575     // tableau de travail : stockage des types geometriques pour UNE entite
 576     vector<medGeometryElement> types;
 577     types.reserve( maillageByType.size() );
 578
 579     // iterator returning mailles of each type of an entity,
 580     // but if hasMixedCells, we iterate on all types at every dim, since
 581     // in this case we store POINT1 elems as MED_NODE and
 582     // elems of all the rest types as MED_CELL
 583     int iterDim = hasMixedCells ? -1 : dim;
 584     _maillageByDimIterator entityMailles( *this, iterDim, /*convertPoly=*/true );
 585
 586     // count nb of types and nb mailles of each type
 587     int dimension=0;
 588     if ( dim == 0 ) {
 589       if ( entityMailles.nextType() && entityMailles.dim() == 0 )
 590       {
 591         count.push_back( count.back() + numberOfNodes );
 592         types.push_back( entityMailles.type() );
 593       }
 594     }
 595     else {
 596       while ( entityMailles.nextType() )
 597       {
 598         //if ( entityMailles.dim() > 3 ) break; // ignore poly
 599
 600         dimension = entityMailles.dim();
 601         if ( dimension == 0 ) continue; // if hasMixedCells, iterator returns all types
 602
 603         count.push_back( count.back() + entityMailles.sizeWithoutMerged() );
 604         types.push_back( entityMailles.type() );
 605       }
 606     }
 607     int numberOfTypes = types.size(); // nombre de types de l'entite
 608     if ( numberOfTypes == 0 )
 609       continue;
 610
 611     if ( dimension == meshDim ) entity=MED_CELL;
 612     else if (dimension==2 )     entity=MED_FACE;
 613     else if (dimension==1 )     entity=MED_EDGE;
 614     else if (dimension==0 )     entity=MED_NODE;
 615
 616     Connectivity = new CONNECTIVITY ( numberOfTypes, entity );
 617     Connectivity->setEntityDimension( dimension );
 618     Connectivity->setNumberOfNodes  ( numberOfNodes );
 619     Connectivity->setGeometricTypes ( &types[0], entity);
 620     Connectivity->setCount          ( &count[0], entity );
 621
 622     int prevNbElems = 1; // in previous type
 623     for (int k=0; k!=numberOfTypes; ++k )
 624     {
 625       set<_maille> & typeMailles = maillageByType[ types[k] ];
 626       i = typeMailles.begin(), iEnd = typeMailles.end();
 627       int nbMailles = count[k+1]-count[k];
 628       int* connectivity = 0, *index = 0;
 629
 630       switch ( types[k] )
 631       {
 632       case MED_POLYGON:
 633         {
 634           // put polygones in order of increasing number
 635           vector<const _maille*> orderedPoly( nbMailles );
 636           for ( ; i != iEnd; ++i )
 637             if ( !i->isMerged() )
 638               orderedPoly[ i->ordre() - prevNbElems ] = &(*i);
 639
 640           // make index
 641           int* polyIndex = index = new int[ nbMailles + 1 ];
 642           vector<const _maille*>::iterator poly = orderedPoly.begin(), polyEnd = orderedPoly.end();
 643           for ( *polyIndex++ = 1; polyIndex < index+nbMailles+1; ++poly, ++polyIndex)
 644             *polyIndex = polyIndex[-1] + (*poly)->sommets.size();
 645
 646           // make connectivity
 647           int nbNodes = polyIndex[-1];
 648           int* conn = connectivity = new int[ nbNodes ];
 649           for ( poly = orderedPoly.begin(); poly != polyEnd; ++poly) {
 650             for ( int j = 0, nbNodes = (*poly)->sommets.size(); j < nbNodes; ++j )
 651               *conn++ = (*poly)->nodeNum( j );
 652           }
 653           break;
 654         }
 655
 656       case MED_POLYHEDRA:
 657         {
 658           if ( typeMailles.size() != polyherdalNbFaceNodes.size() )
 659             throw MEDEXCEPTION (LOCALIZED(STRING(LOC) << "Missing info on polyhedron faces"));
 660
 661           typedef TPolyherdalNbFaceNodes::iterator TPolyFaNoIter;
 662           TPolyFaNoIter polyFaNo, polyFaNoEnd = polyherdalNbFaceNodes.end();
 663
 664           // put poly's in order of increasing number and count size of connectivity
 665           vector<TPolyFaNoIter> orderedPolyFaNo( nbMailles );
 666           int connSize = 0;
 667           for ( polyFaNo = polyherdalNbFaceNodes.begin(); polyFaNo != polyFaNoEnd; ++polyFaNo )
 668             if ( !polyFaNo->first->isMerged() )
 669             {
 670               orderedPolyFaNo[ polyFaNo->first->ordre() - prevNbElems ] = polyFaNo;
 671               connSize += polyFaNo->first->sommets.size() + polyFaNo->second.size() - 1;
 672             }
 673           vector<TPolyFaNoIter>::iterator pfnIter, pfnEnd = orderedPolyFaNo.end();
 674
 675           // make index and connectivity
 676           int* conn = connectivity = new int[ connSize ];
 677           int* ind  = index        = new int[ nbMailles+1 ];
 678           *ind++ = 1;
 679           for ( pfnIter = orderedPolyFaNo.begin(); pfnIter != pfnEnd; ++pfnIter)
 680           {
 681             const _maille * poly = (*pfnIter)->first;
 682             const vector<int> & nbFaceNodes = (*pfnIter)->second;
 683             int nbNodes = 0;
 684             for ( unsigned iFace = 0; iFace < nbFaceNodes.size(); ++iFace )
 685             {
 686               for ( int j = 0, nbFNodes = nbFaceNodes[iFace]; j < nbFNodes; ++j )
 687                 *conn++ = poly->nodeNum( nbNodes++ );
 688               *conn++ = -1;
 689             }
 690             conn--;
 691             *ind = ind[-1] + nbNodes;
 692             ++ind;
 693           }
 694           break;
 695         }
 696
 697       default: // CLASSIC TYPES
 698
 699         // copie des sommets dans connectivity et set dans Connectivity
 700         int nbSommetsParMaille = i->sommets.size();
 701         int nbSommets = nbMailles * nbSommetsParMaille;
 702         connectivity = new int[ nbSommets ];
 703         if ( entity==MED_NODE )
 704         {
 705           for (int l=0; l!=nbSommets; ++l)
 706             connectivity[l] = l+1;
 707         }
 708         else
 709         {
 710           for ( ; i != iEnd; ++i ) { // loop on elements of geom type
 711             if ( i->isMerged() )
 712               continue;
 713             int* mailleConn = connectivity + nbSommetsParMaille * ( i->ordre() - prevNbElems );
 714             if ( i->reverse )
 715               for ( int n=nbSommetsParMaille-1; n!=-1; --n)
 716                 *mailleConn++ = i->nodeNum( n );
 717             else
 718               for ( int n=0; n != nbSommetsParMaille; ++n)
 719                 *mailleConn++ = i->nodeNum( n );
 720           }
 721           // DO NOT ERASE, maillage will be used while fields construction
 722           //maillage.erase(j);    ; // dangereux, mais optimise la memoire consommee!
 723         }
 724       }
 725
 726       Connectivity->setNodal (connectivity, entity, types[k], index);
 727       delete [] connectivity;
 728       delete [] index; index = 0;
 729
 730       prevNbElems += nbMailles;
 731     }
 732
 733     if ( Constituent )
 734       Connectivity->setConstituent (Constituent);
 735     // stocke Connectivity pour utilisation dans setConstituent lors de la boucle suivante
 736     Constituent = Connectivity;
 737
 738     if ( entity == MED_CELL )
 739       break; // necessary if hasMixedCells
 740   }
 741
 742   END_OF_MED(LOC);
 743   return Connectivity;
 744 }
 745
 746
 747 /*!
 748  * \if developper
 749  * fill the arguments vector of groups from the intermediate structure.
 750  * This function must be called before getConnectivity()
 751  * WARNING: new GROUP on all elements are invalid: numbers are not set!
 752  * to make them valid it is necessary to update() them after setting
 753  * connectivity to mesh
 754  * \endif
 755  */
 756 void _intermediateMED::getGroups(vector<GROUP *> & _groupCell,
 757                                  vector<GROUP *> & _groupFace,
 758                                  vector<GROUP *> & _groupEdge,
 759                                  vector<GROUP *> & _groupNode, MESH * _ptrMesh)
 760 {
 761   const char* LOC = "_intermediateMED::getGroups() : ";
 762   BEGIN_OF_MED(LOC);
 763
 764   //medGroupes.resize( groupes.size(), 0 );
 765   if (maillageByType.size() == 0) {
 766     INFOS_MED( "Erreur : il n'y a plus de mailles");
 767     return;
 768   }
 769
 770   // get indices of groups that are field support - do not skip them
 771   std::set<int> support_groups;
 772   std::list< _fieldBase* >::const_iterator fIt = fields.begin();
 773   for ( ; fIt != fields.end(); ++fIt )
 774     (*fIt)->getGroupIds( support_groups, false );
 775
 776   numerotationMaillage(); // Renumerotation des mailles par entite
 777
 778   int dimension_maillage=getMeshDimension();
 779
 780   for (size_t i=0; i!=this->groupes.size(); ++i)
 781   {
 782     _groupe& grp = groupes[i];
 783     if ( grp.medGroup /*medGroupes[ i ]*/ )
 784       continue; // grp already converted into a MED group
 785
 786     bool isFieldSupport = ( support_groups.find( i ) != support_groups.end() );
 787     // si le groupe est vide, ou si le groupe n'est pas nommé : on passe au suivant
 788     if ( grp.empty() || ( grp.nom.empty() && !isFieldSupport )) {
 789       if ( !grp.nom.empty() ) {
 790         INFOS_MED("Skip group " << i << " " << grp.nom);
 791 //       } else {
 792 //         MESSAGE_MED("Skip group " << i << " <" << grp.nom << "> " << ( grp.empty() ? ": empty" : ""));
 793       }
 794       continue;
 795     }
 796
 797     // Build a set of mailles: sort mailles by type and exclude maille doubling
 798     bool isSelfIntersect = false;
 799     typedef set< set<_maille>::iterator, _mailleIteratorCompare > TMailleSet;
 800     TMailleSet mailleSet;
 801     if( grp.groupes.size() ) {// le groupe i contient des sous-maillages
 802       int nb_elem = 0;
 803       for( vector<int>::iterator j=grp.groupes.begin(); j!=grp.groupes.end(); ++j)
 804       {
 805         nb_elem += groupes[*j-1].mailles.size();
 806         _groupe::TMailleIter maIt=groupes[*j-1].mailles.begin();
 807         for( ; maIt!=groupes[*j-1].mailles.end(); ++maIt) {
 808 //           TMailleSet::const_iterator ma_it = mailleSet.find( *maIt );
 809 //           if ( ma_it != mailleSet.end() ) {
 810 //             MESSAGE_MED("EQUAL ELEMS: " << *ma_it << " AND " << *maIt);
 811 //           }
 812 //           else
 813             mailleSet.insert( *maIt );
 814         }
 815       }
 816       if ( nb_elem != (int)mailleSet.size() ) { // Self intersecting compound group
 817         isSelfIntersect = true;
 818         INFOS_MED("Self intersecting group: " << i << " <" << grp.nom << ">"
 819               << ", mailleSet.size = " << mailleSet.size() << ", sum nb elems = " << nb_elem);
 820         for( vector<int>::iterator j=grp.groupes.begin(); j!=grp.groupes.end(); ++j)
 821           INFOS_MED(" in sub-group "<<  *j << " <" << groupes[*j-1].nom << "> "
 822                 << groupes[*j-1].mailles.size() << " mailles of type "
 823                 << groupes[*j-1].mailles[0]->geometricType);
 824         // if a group serve as a support to a field, then the _field is to be converted
 825         // into several MEDMEM::FIELDs, one per sub-group; if we already skipped some sub-groups,
 826         // we roll back the loop on groups to create MED groups for skipped ones
 827         if ( isFieldSupport ) {
 828           if ( grp.nom.empty() ) // clear group existing for the sake of field only
 829             grp.groupes.clear();
 830           std::set<int> sub_grps;
 831           for ( fIt = fields.begin(); fIt != fields.end(); ++fIt ) {
 832             _fieldBase * field = (*fIt);
 833             if ( field->_group_id == (int)i ) {
 834               field->_group_id = -1; // -> a field by support
 835               field->getGroupIds( sub_grps, false );
 836             }
 837           }
 838           if ( (int)i > *sub_grps.begin() ) { // roll back
 839             support_groups.erase( i );
 840             support_groups.insert( sub_grps.begin(), sub_grps.end() );
 841             i = *sub_grps.begin() - 1;
 842             continue;
 843           }
 844         }
 845       }
 846     }
 847     else {
 848       _groupe::TMailleIter maIt=grp.mailles.begin();
 849       for(; maIt!=grp.mailles.end(); ++maIt)
 850         mailleSet.insert( *maIt );
 851       if ( grp.mailles.size() != mailleSet.size() )
 852         INFOS_MED( "Self intersecting group: " << i << " <" << grp.nom << ">"
 853                << ", mailleSet.size = " << mailleSet.size()
 854                << ", nb elems = " << grp.mailles.size());
 855     }
 856     // MEDMEM does not allow constituents of <mesh_dim>-2 and less dimension
 857     // but do not skip node group
 858     int group_min_dim = (**mailleSet.begin()).dimensionWithPoly();
 859     int group_max_dim = (**(--mailleSet.end())).dimensionWithPoly();
 860     if ( group_max_dim != 0 && group_min_dim <= dimension_maillage - 2  ) {
 861       MESSAGE_MED("Skip group: " << i << " <" << grp.nom << "> - too small dimension: "
 862               << group_min_dim);
 863       continue;
 864     }
 865
 866     // initialise groupe_entity a l'entite de la premiere maille du groupe
 867     medEntityMesh groupe_entity = (**mailleSet.rbegin()).getEntity(dimension_maillage);
 868     if ( hasMixedCells && group_min_dim > 0 )
 869       groupe_entity = MED_CELL;
 870
 871     // total nb of elements in mesh of groupe_entity
 872     int totalNbElements = 0;
 873     if ( groupe_entity == MED_NODE )
 874       totalNbElements = points.size() - nbMerged( MED_POINT1 );
 875     else {
 876       int entityDim = hasMixedCells ? -1 : group_min_dim;
 877       _maillageByDimIterator allMailles( *this, entityDim, true );
 878       while ( allMailles.nextType() )
 879         if ( allMailles.dim() > 0 )
 880           totalNbElements += allMailles.sizeWithoutMerged();
 881     }
 882     const bool isOnAll = ((int) mailleSet.size() == totalNbElements );
 883
 884     // if !isOnAll, build a map _maille::ordre() -> index in GROUP.getNumber(MED_ALL_ELEMENTS).
 885     // It is used while fields building.
 886     if ( !isOnAll || isSelfIntersect || isFieldSupport ) {
 887       TMailleSet::iterator maIt = mailleSet.begin();
 888       for ( int iMa = 0; maIt != mailleSet.end(); maIt++ )
 889         grp.relocMap.insert( make_pair( (*maIt)->ordre(), ++iMa ));
 890     }
 891     //Parcours des mailles (a partir de la deuxieme) pour compter les types geometriques
 892     int nb_geometric_types=1;
 893     TMailleSet::iterator j=mailleSet.begin();
 894     medGeometryElement geometrictype=(**j).geometricType;
 895     for ( ++j ; j!=mailleSet.end(); ++j )
 896     {
 897       //Compte nombre de types geometriques
 898       if ( (**j).geometricType != geometrictype ) // si on change de type geometrique
 899       {
 900         nb_geometric_types++;
 901         geometrictype=(**j).geometricType;
 902       }
 903     }
 904
 905     MED_EN::medGeometryElement * tab_types_geometriques = new MED_EN::medGeometryElement[nb_geometric_types];
 906     int * tab_index_types_geometriques = new int[nb_geometric_types+1];
 907     int * tab_numeros_elements = new int[mailleSet.size()];
 908     int * tab_nombres_elements = new int[nb_geometric_types];
 909
 910     //Remplit tableaux entree des methodes set
 911     int indice_mailles=0/*, maxOrdre = -1*/;
 912     j=mailleSet.begin();
 913     geometrictype=(**j).geometricType;
 914     tab_index_types_geometriques[0]=1;
 915     int indice_types_geometriques=1;
 916     tab_types_geometriques[0]=geometrictype;
 917     //parcours des mailles du groupe
 918     for (  ; j!=mailleSet.end(); ++j , ++indice_mailles)
 919     {
 920       const _maille& ma = **j;
 921       tab_numeros_elements[indice_mailles]= ma.ordre();
 922 //       if ( maxOrdre < tab_numeros_elements[indice_mailles] )
 923 //         maxOrdre = tab_numeros_elements[indice_mailles];
 924       if (ma.geometricType != geometrictype)
 925       {
 926         tab_index_types_geometriques[indice_types_geometriques]=indice_mailles+1;
 927         geometrictype=ma.geometricType;
 928         tab_types_geometriques[indice_types_geometriques]=geometrictype;
 929         ++indice_types_geometriques;
 930       }
 931     }
 932     tab_index_types_geometriques[indice_types_geometriques]=indice_mailles+1;
 933     for (int k=0; k != nb_geometric_types; ++k)
 934     {
 935       tab_nombres_elements[k] = tab_index_types_geometriques[k+1]-tab_index_types_geometriques[k];
 936     }
 937     //INFOS_MED( "MAX ORDRE in grp " << grp.nom << " entity " << groupe_entity << " : " << maxOrdre);
 938
 939     //Determination type entite du groupe
 940     vector <GROUP *> * vect_group;
 941     switch ( groupe_entity )
 942     {
 943     case MED_CELL :
 944       vect_group= & _groupCell;
 945       break;
 946     case MED_FACE :
 947       vect_group= & _groupFace;
 948       break;
 949     case MED_EDGE :
 950       vect_group= & _groupEdge;
 951       break;
 952     case MED_NODE :
 953       vect_group= & _groupNode;
 954       break;
 955     }
 956
 957     //Creation nouveau groupe MED
 958     GROUP * new_group = new GROUP();
 959     grp.medGroup = new_group;
 960     //Appel methodes set
 961     //new_group->setTotalNumberOfElements(mailleSet.size());
 962     new_group->setName(grp.nom);
 963     new_group->setMesh(_ptrMesh);
 964     if ( _ptrMesh )
 965       _ptrMesh->removeReference();
 966     new_group->setNumberOfGeometricType(nb_geometric_types);
 967     new_group->setGeometricType(tab_types_geometriques);
 968     new_group->setNumberOfElements(tab_nombres_elements);
 969     new_group->setNumber(tab_index_types_geometriques,tab_numeros_elements);
 970     new_group->setEntity(groupe_entity);
 971     new_group->setAll( isOnAll );
 972
 973     vect_group->push_back(new_group);
 974
 975     // Issue 0021311. Use case: a gibi group has references (recorded in pile 1)
 976     // and several names (pile 27) refer (pile 10) to this group.
 977     // We create a copy of this group per each named reference
 978     for ( unsigned iRef = 0 ; iRef < grp.refNames.size(); ++iRef )
 979       if ( !grp.refNames[ iRef ].empty() )
 980       {
 981         vect_group->push_back( new GROUP( *new_group ));
 982         vect_group->back()->setName( grp.refNames[ iRef ] );
 983       }
 984
 985     delete [] tab_types_geometriques;
 986     delete [] tab_index_types_geometriques;
 987     delete [] tab_numeros_elements;
 988     delete [] tab_nombres_elements;
 989   }
 990
 991   END_OF_MED(LOC);
 992 }
 993
 994 //=======================================================================
 995 //function : getFamilies
 996 //purpose  : create families like MESH::createFamilies() but preserves
 997 //           the order of elements in GROUPs defined by constituent families
 998 //           order. Call it after getGroups()
 999 //=======================================================================
1000
1001 // void _intermediateMED::getFamilies(std::vector<FAMILY *> & _famCell,
1002 //                                    std::vector<FAMILY *> & _famFace,
1003 //                                    std::vector<FAMILY *> & _famEdge,
1004 //                                    std::vector<FAMILY *> & _famNode, MESH * _ptrMesh)
1005 // {
1006 //   const char * LOC = "_intermediateMED::getFamilies() : ";
1007 //   BEGIN_OF_MED(LOC);
1008
1009 //   int nbElemFam = 0, nbNodeFam = 0;
1010 //   std::map< GROUP*, vector< FAMILY * > > grpFamsMap;
1011 //   int dimension_maillage=maillage.rbegin()->dimension();
1012
1013 //   std::set<_maille>::const_iterator i=maillage.begin(); // iterateurs sur les mailles
1014 //   std::set<_maille>::const_iterator j=maillage.begin();
1015
1016 //   do
1017 //   {
1018 //     // make a family containing mailles shared by the same set of groups
1019 //     std::list<unsigned>&  grList = i->groupes;  // to define the family end
1020 //     int           dimension = i->dimension();        // to define the entity end
1021 //     medGeometryElement type = i->geometricType;
1022 //     medEntityMesh    entity = i->getEntity( dimension_maillage );
1023
1024 //     std::vector<medGeometryElement> tab_types_geometriques;
1025 //     std::vector<int> tab_index_types_geometriques;
1026 //     std::vector<int> tab_nombres_elements;
1027 //     std::vector<int> tab_numeros_elements;
1028
1029 //     int iMa = 1, nbtype = 0;
1030 //     tab_types_geometriques.push_back( type );
1031 //     tab_index_types_geometriques.push_back( iMa );
1032
1033 //     // scan family cells and fill the tab that are needed by the create a MED FAMILY
1034 //     while (i != maillage.end() &&
1035 //            i->groupes == grList &&
1036 //            i->dimension() == dimension)
1037 //     {
1038 //       if (type != i->geometricType) // si changement de type geometrique
1039 //       {
1040 //         tab_index_types_geometriques.push_back(iMa);
1041 //         tab_nombres_elements.push_back(nbtype);
1042 //         nbtype=0;
1043 //         type=i->geometricType;
1044 //         tab_types_geometriques.push_back(type); // stocke le nouveau type geometrique rencontre
1045 //       }
1046 //       ++nbtype;
1047 //       ++iMa;
1048 //       ++i;
1049 //     }
1050 //     tab_index_types_geometriques.push_back(iMa);
1051 //     tab_nombres_elements.push_back(nbtype); // n'a pas été stocké dans la boucle
1052
1053 //     tab_numeros_elements.resize( iMa - 1 );
1054 //     for ( iMa = 0; j != i; j++, iMa++ )
1055 //       tab_numeros_elements[ iMa ] = j->ordre;
1056
1057 //     int id = ( entity == MED_NODE ? ++nbNodeFam : -(++nbElemFam) );
1058
1059 //     ostringstream name;
1060 //     name << "FAM_" << id;
1061
1062 //     // create a empty MED FAMILY and fill it with the tabs we constructed
1063 //     FAMILY* newFam = new FAMILY();
1064 //     newFam->setTotalNumberOfElements( iMa );
1065 //     newFam->setName( name.str() );
1066 //     newFam->setMesh( _ptrMesh );
1067 //     newFam->setNumberOfGeometricType( tab_types_geometriques.size() );
1068 //     newFam->setGeometricType( &tab_types_geometriques[0] ); // we know the tab is not empy
1069 //     newFam->setNumberOfElements( &tab_nombres_elements[0] );
1070 //     newFam->setNumber( &tab_index_types_geometriques[0], &tab_numeros_elements[0] );
1071 //     newFam->setEntity( entity );
1072 //     newFam->setAll( false );
1073 //     newFam->setIdentifier( id );
1074 //     newFam->setNumberOfGroups( grList.size() );
1075
1076 //     // Update links between families and groups
1077 //     if ( ! grList.empty() )
1078 //     {
1079 //       std::string * groupNames = new string[ grList.size() ];
1080 //       std::list<unsigned>::iterator g = grList.begin();
1081 //       for ( int i = 0; g != grList.end(); ++g, ++i ) {
1082 //         GROUP * medGROUP = getGroup( *g );
1083 //         groupNames[ i ] = medGROUP->getName();
1084 //         grpFamsMap[ medGROUP ].push_back( newFam );
1085 //       }
1086 //       newFam->setGroupsNames(groupNames);
1087 //     }
1088 //     // store newFam
1089 //     std::vector<FAMILY*>* families = 0;
1090 //     switch ( entity )
1091 //     {
1092 //     case MED_CELL :
1093 //       families = & _famCell; break;
1094 //     case MED_FACE :
1095 //       families = & _famFace; break;
1096 //     case MED_EDGE :
1097 //       families = & _famEdge; break;
1098 //     case MED_NODE :
1099 //       families = & _famNode; break;
1100 //     }
1101 //     if ( families )
1102 //       families->push_back( newFam );
1103
1104 //   } while ( i != maillage.end() );
1105
1106 //   // update references in groups
1107 //   std::map< GROUP*, vector< FAMILY * > >::iterator gf = grpFamsMap.begin();
1108 //   for ( ; gf != grpFamsMap.end(); ++gf )
1109 //   {
1110 //     gf->first->setNumberOfFamilies( gf->second.size() );
1111 //     gf->first->setFamilies( gf->second );
1112 //   }
1113 //}
1114
1115 //=======================================================================
1116 //function : getGroup
1117 //purpose  :
1118 //=======================================================================
1119
1120 // GROUP * _intermediateMED::getGroup( int i )
1121 // {
1122 //   if ( i <(int) medGroupes.size() )
1123 //     return medGroupes[ i ];
1124 //   throw MEDEXCEPTION
1125 //     (LOCALIZED(STRING("_intermediateMED::getGroup(): WRONG GROUP INDEX: ")
1126 //                << medGroupes.size() << " <= " << i ));
1127 // }
1128
1129 //=======================================================================
1130 //function : getFields
1131 //purpose  :
1132 //=======================================================================
1133
1134 void _intermediateMED::getFields(std::list< FIELD_* >& theFields)
1135 {
1136   const char * LOC = "_intermediateMED::getFields() : ";
1137   BEGIN_OF_MED(LOC);
1138
1139   std::list< _fieldBase* >::const_iterator fIt = fields.begin();
1140   for ( ; fIt != fields.end(); fIt++ )
1141   {
1142     const _fieldBase* fb = *fIt;
1143     list<pair< FIELD_*, int> >  ff = fb->getField(groupes);
1144     list<pair< FIELD_*, int> >::iterator f_sup = ff.begin();
1145     for (int j = 1 ; f_sup != ff.end(); f_sup++, ++j )
1146     {
1147       FIELD_    * f = f_sup->first;
1148       SUPPORT * sup = groupes[ f_sup->second ].medGroup;
1149       if ( !sup )
1150         throw MEDEXCEPTION
1151           (LOCALIZED(STRING(LOC) <<"_intermediateMED::getFields(), NULL field support: "
1152                      << " group index: " << fb->_group_id));
1153       int nb_elems = sup->getNumberOfElements( MED_ALL_ELEMENTS );
1154       if ( nb_elems != f->getNumberOfValues() )
1155         throw MEDEXCEPTION
1156           (LOCALIZED(STRING("_intermediateMED::getFields(), field support size (")
1157                      << nb_elems  << ") != NumberOfValues (" << f->getNumberOfValues()));
1158       theFields.push_back( f );
1159       if ( sup->getName().empty() ) {
1160         ostringstream name;
1161         name << "GRP_" << f->getName() << "_" << j;
1162         sup->setName( name.str() );
1163       }
1164       f->setSupport( sup );
1165       //f->setIterationNumber( j );
1166       f->setOrderNumber( j );
1167     }
1168   }
1169   END_OF_MED(LOC);
1170 }
1171
1172 //=======================================================================
1173 //function : ~_intermediateMED
1174 //purpose  : destructor
1175 //=======================================================================
1176
1177 _intermediateMED::~_intermediateMED()
1178 {
1179   MESSAGE_MED( "~_intermediateMED()");
1180   std::list< _fieldBase* >::const_iterator fIt = fields.begin();
1181   for ( ; fIt != fields.end(); fIt++ )
1182     delete *fIt;
1183 }
1184
1185 //=======================================================================
1186 //function : getGroupIds
1187 //purpose  : return ids of main and/or sub- groups
1188 //=======================================================================
1189
1190 void _fieldBase::getGroupIds( std::set<int> & ids, bool all ) const
1191 {
1192   if ( hasCommonSupport() )
1193     ids.insert( _group_id );
1194   if ( all || !hasCommonSupport() ) {
1195     vector< _sub_data >::const_iterator sub = _sub.begin();
1196     for ( ; sub != _sub.end(); ++sub )
1197       ids.insert( sub->_supp_id );
1198   }
1199 }
1200
1201 //=======================================================================
1202 //function : hasSameComponentsBySupport
1203 //purpose  :
1204 //=======================================================================
1205
1206 bool _fieldBase::hasSameComponentsBySupport() const
1207 {
1208   vector< _sub_data >::const_iterator sub_data = _sub.begin();
1209   const _sub_data& first_sub_data = *sub_data;
1210   for ( ++sub_data ; sub_data != _sub.end(); ++sub_data )
1211   {
1212     if ( first_sub_data._comp_names != sub_data->_comp_names )
1213       return false; // diff names of components
1214
1215     if ( first_sub_data._nb_gauss != sub_data->_nb_gauss )
1216       return false; // diff nb of gauss points
1217   }
1218   return true;
1219 }
1220
1221 //=======================================================================
1222 //
1223 //function : mergeNodes
1224 //
1225 //=======================================================================
1226 namespace {
1227
1228 struct __NOEUD
1229 {
1230   typedef map<int,_noeud>::iterator Inoeud;
1231   Inoeud i_noeud;
1232
1233   __NOEUD() {}
1234   __NOEUD(Inoeud n): i_noeud(n) {}
1235   const double & operator[] (int i) const { return i_noeud->second.coord[i]; }
1236   double         operator[] (int i)       { return i_noeud->second.coord[i]; }
1237   int dim() const { return i_noeud->second.coord.size(); }
1238   int& num() { return i_noeud->second.number; }
1239   int id() const { return i_noeud->first; }
1240   bool isMerged() const { return i_noeud->second.number < 1; }
1241 };
1242 //-----------------------------------------------------------------------
1243 double DistanceL2(const __NOEUD &A,const __NOEUD &B)
1244 {
1245   if ( B.isMerged() ) return DBL_MAX;
1246   const double* cooA = &A[0], *cooB = &B[0], *cooEnd = cooA + A.dim();
1247   double dist, somme=0;
1248   while ( cooA < cooEnd ) {
1249     dist=((*cooA++) - (*cooB++));
1250     somme+=dist*dist;
1251   }
1252   return sqrt(somme);
1253 }
1254 //-----------------------------------------------------------------------
1255 struct __NUAGENOEUD
1256 {
1257   vector< __NOEUD > nodes;
1258   __NUAGENOEUD(_intermediateMED& imed);
1259   __NOEUD & operator [] (int i) { return nodes[i]; }
1260   int size() const { return nodes.size(); }
1261 };
1262 //-----------------------------------------------------------------------
1263 __NUAGENOEUD::__NUAGENOEUD(_intermediateMED& imed)
1264 {
1265   nodes.resize( imed.points.size() );
1266   map<int,_noeud>::iterator i_noeud = imed.points.begin();
1267   for( int i = 0; i_noeud!=imed.points.end(); ++i_noeud, ++i ) {
1268     i_noeud->second.number = i+1;
1269     nodes[ i ] = i_noeud;
1270   }
1271 }
1272 //-----------------------------------------------------------------------
1273 // mergeNodes
1274 template<int DIM> int mergeNodes(double            tolerance,
1275                                  _intermediateMED& imed,
1276                                  vector< int > &   /*newNodeIDs*/)
1277 {
1278   /*typedef dTree<__NOEUD,__NUAGENOEUD,DIM > DTree;
1279   __NUAGENOEUD aNUAGENOEUD( imed );
1280   DTree tree( &aNUAGENOEUD );
1281
1282 //   int maxID = imed.points.rbegin()->first;
1283 //   newNodeIDs.resize( maxID + 1, 0 );
1284
1285   int num = 0, nbRemoved = 0;
1286   int nbNodes = aNUAGENOEUD.size();
1287   for ( int i = 0; i < nbNodes; ++i )
1288   {
1289     __NOEUD & node = aNUAGENOEUD[i];
1290     int & number = node.num();
1291     if ( number < 0 )
1292       continue; // already merged
1293     number = ++num;
1294
1295     list<int> closeNumbers;
1296     int nbCoinsident = tree.get_all_close( node, tolerance, closeNumbers );
1297     if ( nbCoinsident > 1 ) // not only node it-self found
1298     {
1299       nbRemoved += nbCoinsident-1; // one of coincident nodes remains
1300       int id = node.id();
1301       list<int>::iterator n = closeNumbers.begin(), nEnd = closeNumbers.end();
1302       while ( n != nEnd ) {
1303         __NOEUD & coincNode = aNUAGENOEUD[ *n++ ];
1304         int coincID = coincNode.id();
1305         if ( coincID != id ) {
1306           coincNode.num() = -number;
1307           //newNodeIDs[ coincID ] = id;
1308         }
1309       }
1310     }
1311   }
1312   return nbRemoved;*/
1313   return 0;
1314 }
1315 //-----------------------------------------------------------------------
1316 // wrapper of _maille used after merging nodes to find equal mailles
1317 struct _mailleWrap
1318 {
1319   const _maille* _ma;
1320
1321   _mailleWrap(const _maille* ma=0): _ma(ma) { if (_ma) _ma->init(); }
1322   ~_mailleWrap()                            { if (_ma) _ma->init(); }
1323
1324   bool operator < (const _mailleWrap& mw) const
1325   {
1326     size_t l=_ma->sommets.size();
1327     const int* v1 = getSortedNodeNums(_ma);
1328     const int* v2 = getSortedNodeNums(mw._ma);
1329     for ( const int* vEnd = v1 + l; v1 < vEnd; ++v1, ++v2 )
1330       if(*v1 != *v2)
1331         return *v1 < *v2;
1332     return false; // cas d'égalité
1333   }
1334   static const int* getSortedNodeNums(const _maille* ma)
1335   {
1336     if ( !ma->sortedNodeIDs ) {
1337       size_t l=ma->sommets.size();
1338       ma->sortedNodeIDs = new int[ l ];
1339       for (size_t i=0; i!=l; ++i)
1340         ma->sortedNodeIDs[i]=ma->nodeNum(i);
1341       std::sort( ma->sortedNodeIDs, ma->sortedNodeIDs + l );
1342     }
1343     return ma->sortedNodeIDs;
1344   }
1345 };
1346
1347 } // namespace
1348
1349 //=======================================================================
1350 //function : mergeNodesAndElements
1351 //purpose  : optionally merge nodes and elements
1352 //=======================================================================
1353
1354 void _intermediateMED::mergeNodesAndElements(double tolerance)
1355 {
1356   //const char * LOC = "_intermediateMED::mergeNodesAndElements(): ";
1357   vector< int > newNodeIDs; // for i-th node id, id to replace with
1358
1359   int nbRemovedNodes = 0;
1360   const int spaceDimension=points.begin()->second.coord.size();
1361   if ( spaceDimension == 3 )
1362     nbRemovedNodes = mergeNodes<3>( tolerance, *this, newNodeIDs );
1363   else if ( spaceDimension == 2 )
1364     nbRemovedNodes = mergeNodes<2>( tolerance, *this, newNodeIDs );
1365
1366   myNodesNumerated = true;
1367
1368   if ( nbRemovedNodes == 0 )
1369     return;
1370
1371   // all commented code here was used to keep initial numeration but it was too slow
1372   //numerotationMaillage();
1373   treatGroupes();
1374
1375   nbRemovedByType[ MED_NONE   ] = nbRemovedNodes;
1376   nbRemovedByType[ MED_POINT1 ] = nbRemovedNodes;
1377
1378   bool hasPointMailles = false;
1379   _maillageByDimIterator entityMailles( *this, 0 );
1380   if ( const set<_maille> * points = entityMailles.nextType() ) {
1381     hasPointMailles = true;
1382     set<_maille>::const_iterator i, iEnd = points->end();
1383     for ( i = points->begin(); i != iEnd; ++i )
1384       i->setOrdre( i->sommets[0]->second.number ); // is merged if point is merged
1385   }
1386   const bool skipFirstType = ( hasPointMailles && hasMixedCells );
1387   // loop on entities
1388   for ( int dim = 1; dim <= 3; ++dim )
1389   {
1390     int iterDim = hasMixedCells ? -1 : dim;
1391     //int nbRemovedInEntity = 0;
1392
1393     // count total nb of mailles in entity
1394 //     int nbMailles = 0;
1395 //     entityMailles.init( iterDim, true );
1396 //     if ( skipFirstType ) entityMailles.nextType(); // merged POINT1's are same as nodes
1397 //     while ( const set<_maille> * typeMailles = entityMailles.nextType() )
1398 //       nbMailles += typeMailles->size();
1399
1400     // for each maille number, count shift due to removed mailles with lower numbers
1401     //vector<int> shift( nbMailles+1, 0 );
1402
1403     // merge and numerate mailles
1404     int num = 0;
1405     entityMailles.init( iterDim, true );
1406     if ( skipFirstType ) entityMailles.nextType(); // merged POINT1's are same as nodes
1407     while ( const set<_maille> * typeMailles = entityMailles.nextType() )
1408     {
1409       int nbRemoved = 0;
1410       set<_mailleWrap> maNodeNumSorted;
1411       pair< set<_mailleWrap>::const_iterator, bool > mw_isunique;
1412
1413       set<_maille>::const_iterator ma = typeMailles->begin(), maEnd = typeMailles->end();
1414       while ( ma != maEnd )
1415       {
1416         const _maille & m = *ma++;
1417         mw_isunique = maNodeNumSorted.insert( &m );
1418         if ( mw_isunique.second ) {
1419           m.setOrdre( ++num );
1420         }
1421         else {
1422           const _maille* equalMa = mw_isunique.first->_ma;
1423           //unsigned ordreToRemove = m.ordre();
1424           m.setMergedOrdre( equalMa->ordre() );
1425           nbRemoved++;
1426 //           while ( ordreToRemove <= nbMailles )
1427 //             shift[ ordreToRemove++ ]++;
1428         }
1429       }
1430
1431       if ( nbRemoved ) {
1432         nbRemovedByType[ entityMailles.type() ] = nbRemoved;
1433         //nbRemovedInEntity += nbRemoved;
1434       }
1435
1436       // update maille ordre
1437 //       if ( nbRemovedInEntity ) {
1438 //         for ( ma = typeMailles->begin(); ma != maEnd; ++ma ) {
1439 //           unsigned newOrdre = ma->ordre() - shift[ ma->ordre() ];
1440 //           if ( ma->isMerged() )
1441 //             ma->setMergedOrdre( newOrdre );
1442 //           else
1443 //             ma->setOrdre( newOrdre );
1444 //         }
1445 //       }
1446     }
1447   }
1448
1449   myMaillesNumerated = true;
1450
1451 }
1452
1453 //=======================================================================
1454 //function : getMeshDimension
1455 //purpose  : count mesh dimension
1456 //=======================================================================
1457
1458 int _intermediateMED::getMeshDimension() const
1459 {
1460   int dim = 0;
1461   _maillageByDimIterator allMailles( *this, -1, true );
1462   while ( allMailles.nextType() )
1463     dim = allMailles.dim();
1464   return dim;
1465 }
1466
1467 }
1468
1469 /////