src/MEDMEM/MEDMEM_DriverTools.cxx

   1 //  Copyright (C) 2007-2008  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 #include "MEDMEM_DriverTools.hxx"
  23 #include "MEDMEM_STRING.hxx"
  24 #include "MEDMEM_Exception.hxx"
  25 #include "MEDMEM_Mesh.hxx"
  26 #include "MEDMEM_Group.hxx"
  27 #include "MEDMEM_Field.hxx"
  28 #include "MEDMEM_InterpolationHighLevelObjects.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), _ordre(ma._ordre),
  43     reverse(ma.reverse), sortedNodeIDs(0)
  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 < sommets.size() );
 112   int i2 = ( i + 1 == 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( int i=1; i!=ma.sommets.size(); ++i)
 163         os << ", " << ma.nodeNum( i );
 164     os << " > sortedNodeIDs: ";
 165     if ( ma.sortedNodeIDs ) {
 166       os << "< ";
 167       for( int 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 != 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   // STANDARD types connectivity
 563
 564   // loop on entities
 565   for ( int dim = 0; dim <= 3; ++dim )
 566   {
 567     // skip nodes and elements of <dimension_maillage - 2> or less dimension
 568     // Unfortunately, it is impossible because of MESH::createFamilies() that requires
 569     // presence of connectivity even for nodes!
 570     //int dimension_maillage_moin_2=maillage.rbegin()->dimension() - 2;
 571
 572     // tableau de travail : nombre d'elements pour chaque type geometrique
 573     vector<int> count;
 574     count.reserve( maillageByType.size() );
 575     count.push_back( 1 );
 576
 577     // tableau de travail : stockage des types geometriques pour UNE entite
 578     vector<medGeometryElement> types;
 579     types.reserve( maillageByType.size() );
 580
 581     // iterator returning mailles of each type of an entity,
 582     // but if hasMixedCells, we iterate on all types at every dim, since
 583     // in this case we store POINT1 elems as MED_NODE and
 584     // elems of all the rest types as MED_CELL
 585     int iterDim = hasMixedCells ? -1 : dim;
 586     _maillageByDimIterator entityMailles( *this, iterDim );
 587
 588     // count nb of types and nb mailles of each type
 589     int dimension=0;
 590     if ( dim == 0 ) {
 591       if ( entityMailles.nextType() && entityMailles.dim() == 0 )
 592       {
 593         count.push_back( count.back() + numberOfNodes );
 594         types.push_back( MED_POINT1 );
 595       }
 596     }
 597     else {
 598       while ( entityMailles.nextType() )
 599       {
 600         if ( entityMailles.dim() > 3 ) break; // ignore poly
 601
 602         dimension = entityMailles.dim();
 603         if ( dimension == 0 ) continue; // if hasMixedCells, iterator returns all types
 604
 605         count.push_back( count.back() + entityMailles.sizeWithoutMerged() );
 606         types.push_back( entityMailles.type() );
 607       }
 608     }
 609     int numberOfTypes = types.size(); // nombre de types de l'entite
 610     if ( numberOfTypes == 0 )
 611       continue;
 612
 613     if ( dimension == meshDim ) entity=MED_CELL;
 614     else if (dimension==2 )     entity=MED_FACE;
 615     else if (dimension==1 )     entity=MED_EDGE;
 616     else if (dimension==0 )     entity=MED_NODE;
 617
 618     Connectivity = new CONNECTIVITY ( numberOfTypes, entity );
 619     Connectivity->setEntityDimension( dimension );
 620     Connectivity->setNumberOfNodes  ( numberOfNodes );
 621     Connectivity->setGeometricTypes ( &types[0], entity);
 622     Connectivity->setCount          ( &count[0], entity );
 623
 624     int prevNbElems = 1; // in previous type
 625     for (int k=0; k!=numberOfTypes; ++k )
 626     {
 627       set<_maille> & typeMailles = maillageByType[ types[k] ];
 628       i = typeMailles.begin(), iEnd = typeMailles.end();
 629       // copie des sommets dans connectivity et set dans Connectivity
 630       int nbSommetsParMaille = i->sommets.size();
 631       int nbMailles = count[k+1]-count[k];
 632       int nbSommets = nbMailles * nbSommetsParMaille;
 633       int* connectivity = new int[ nbSommets ];
 634       if ( entity==MED_NODE ) {
 635         for (int l=0; l!=nbSommets; ++l) {
 636           connectivity[l] = l+1;
 637         }
 638       }
 639       else {
 640         for ( ; i != iEnd; ++i ) { // loop on elements of geom type
 641           if ( i->isMerged() )
 642             continue;
 643           int* mailleConn = connectivity + nbSommetsParMaille * ( i->ordre() - prevNbElems );
 644           if ( i->reverse )
 645             for ( int n=nbSommetsParMaille-1; n!=-1; --n)
 646               *mailleConn++ = i->nodeNum( n );
 647           else
 648             for ( int n=0; n != nbSommetsParMaille; ++n)
 649               *mailleConn++ = i->nodeNum( n );
 650         }
 651         // DO NOT ERASE, maillage will be used while fields construction
 652         //maillage.erase(j);    ; // dangereux, mais optimise la mémoire consommée!
 653       }
 654       Connectivity->setNodal (connectivity, entity, types[k]);
 655       delete [] connectivity;
 656       prevNbElems += nbMailles;
 657     }
 658
 659     if ( Constituent )
 660       Connectivity->setConstituent (Constituent);
 661     // stocke Connectivity pour utilisation dans setConstituent lors de la boucle suivante
 662     Constituent = Connectivity;
 663
 664     if ( entity == MED_CELL )
 665       break; // necessary if hasMixedCells
 666   }
 667
 668   // POLYGONAL connectivity
 669
 670   set<_maille > & polygones = maillageByType[ MED_POLYGON ];
 671   if ( !polygones.empty() )
 672   {
 673     // create connectivity if necessary
 674     entity = ( meshDim == 2 ) ? MED_CELL : MED_FACE;
 675     if ( !Connectivity || Connectivity->getEntity() > entity ) {
 676       Connectivity = new CONNECTIVITY ( 0, entity );
 677       Connectivity->setEntityDimension( 2 );
 678       Connectivity->setNumberOfNodes  ( numberOfNodes );
 679       if ( Constituent )
 680         Connectivity->setConstituent (Constituent);
 681       Constituent = Connectivity;
 682     }
 683
 684     // put polygones in order of increasing number
 685     int numShift = 1 + Connectivity->getNumberOf( entity, MED_ALL_ELEMENTS );
 686     int nbPoly = polygones.size() - nbMerged( MED_POLYGON );
 687     vector<const _maille*> orderedPoly( nbPoly );
 688     for ( i = polygones.begin(), iEnd = polygones.end(); i != iEnd; ++i )
 689       if ( !i->isMerged() )
 690         orderedPoly[ i->ordre() - numShift ] = &(*i);
 691
 692     // make index
 693     vector<int> polyIndex;
 694     polyIndex.reserve( nbPoly + 1 );
 695     vector<const _maille*>::iterator poly = orderedPoly.begin(), polyEnd = orderedPoly.end();
 696     for ( polyIndex.push_back( 1 ); poly != polyEnd; ++poly)
 697       polyIndex.push_back( polyIndex.back() + (*poly)->sommets.size() );
 698
 699     // make connectivity
 700     int nbNodes = polyIndex.back() - 1;
 701     vector<int> polyConn( nbNodes );
 702     vector<int>::iterator conn = polyConn.begin();
 703     for ( poly = orderedPoly.begin(); poly != polyEnd; ++poly) {
 704       for ( int j = 0, nbNodes = (*poly)->sommets.size(); j < nbNodes; ++j )
 705         *conn++ = (*poly)->nodeNum( j );
 706     }
 707     Connectivity->setPolygonsConnectivity(MED_NODAL, entity,
 708                                           &polyConn[0], &polyIndex[0],
 709                                           polyConn.size(), nbPoly );
 710   }
 711
 712   // POLYHEDRAL connectivity
 713
 714   set<_maille > & pHedra = maillageByType[ MED_POLYHEDRA ];
 715   if ( !pHedra.empty() )
 716   {
 717     if ( pHedra.size() != polyherdalNbFaceNodes.size() )
 718       throw MEDEXCEPTION (LOCALIZED(STRING(LOC) << "Missing info on polyhedron faces"));
 719
 720     // create connectivity if necessary
 721     entity = MED_CELL;
 722     if ( !Connectivity || Connectivity->getEntity() != entity ) {
 723       Connectivity = new CONNECTIVITY ( 0, entity );
 724       Connectivity->setEntityDimension( 3 );
 725       Connectivity->setNumberOfNodes  ( numberOfNodes );
 726       if ( Constituent )
 727         Connectivity->setConstituent (Constituent);
 728     }
 729     typedef TPolyherdalNbFaceNodes::iterator TPolyFaNoIter;
 730     TPolyFaNoIter polyFaNo, polyFaNoEnd = polyherdalNbFaceNodes.end();
 731
 732     // put poly's in order of increasing number
 733     int numShift = 1 + Connectivity->getNumberOf( entity, MED_ALL_ELEMENTS );
 734     int nbPoly = pHedra.size() - nbMerged( MED_POLYHEDRA );
 735     vector<TPolyFaNoIter> orderedPolyFaNo( nbPoly );
 736     for ( polyFaNo = polyherdalNbFaceNodes.begin(); polyFaNo != polyFaNoEnd; ++polyFaNo )
 737       if ( !polyFaNo->first->isMerged() )
 738         orderedPolyFaNo[ polyFaNo->first->ordre() - numShift ] = polyFaNo;
 739
 740     vector<TPolyFaNoIter>::iterator pfnIter, pfnEnd = orderedPolyFaNo.end();
 741
 742     // make index pointing to faces of a polyhedron
 743     vector<int> polyIndex;
 744     polyIndex.reserve( nbPoly + 1 );
 745     polyIndex.push_back( 1 );
 746     for ( pfnIter = orderedPolyFaNo.begin(); pfnIter != pfnEnd; ++pfnIter) {
 747       int nbFaces = (*pfnIter)->second.size();
 748       polyIndex.push_back( polyIndex.back() + nbFaces );
 749     }
 750
 751     // make face index pointing to nodes of a face
 752     int nbFaces = polyIndex.back() - 1;
 753     vector<int> faceIndex;
 754     faceIndex.reserve( polyIndex.back() );
 755     faceIndex.push_back( 1 );
 756     for ( pfnIter = orderedPolyFaNo.begin(); pfnIter != pfnEnd; ++pfnIter) {
 757       vector<int> & faceNodes = (*pfnIter)->second;
 758       vector<int>::iterator nbNodes = faceNodes.begin(), nbEnd = faceNodes.end();
 759       for ( ; nbNodes != nbEnd; ++nbNodes )
 760         faceIndex.push_back( faceIndex.back() + *nbNodes );
 761     }
 762
 763     // make connectivity
 764     int nbNodes = faceIndex.back() - 1;
 765     vector<int> polyConn( nbNodes );
 766     vector<int>::iterator conn = polyConn.begin();
 767     for ( pfnIter = orderedPolyFaNo.begin(); pfnIter != pfnEnd; ++pfnIter) {
 768       const _maille * poly = (*pfnIter)->first;
 769       for ( int j = 0, nbNodes = poly->sommets.size(); j < nbNodes; ++j )
 770         *conn++ = poly->nodeNum( j );
 771     }
 772     Connectivity->setPolyhedronConnectivity(MED_NODAL,
 773                                             &polyConn[0], &polyIndex[0],
 774                                             polyConn.size(), nbPoly,
 775                                             &faceIndex[0], nbFaces );
 776   }
 777
 778   END_OF_MED(LOC);
 779   return Connectivity;
 780 }
 781
 782
 783 /*!
 784  * \if developper
 785  * fill the arguments vector of groups from the intermediate structure.
 786  * This function must be called before getConnectivity()
 787  * \endif
 788  */
 789 void _intermediateMED::getGroups(vector<GROUP *> & _groupCell,
 790                                  vector<GROUP *> & _groupFace,
 791                                  vector<GROUP *> & _groupEdge,
 792                                  vector<GROUP *> & _groupNode, MESH * _ptrMesh)
 793 {
 794   const char* LOC = "_intermediateMED::getGroups() : ";
 795   BEGIN_OF_MED(LOC);
 796
 797   //medGroupes.resize( groupes.size(), 0 );
 798   if (maillageByType.size() == 0) {
 799     INFOS_MED( "Erreur : il n'y a plus de mailles");
 800     return;
 801   }
 802
 803   // get indices of groups that are field support - do not skip them
 804   std::set<int> support_groups;
 805   std::list< _fieldBase* >::const_iterator fIt = fields.begin();
 806   for ( ; fIt != fields.end(); ++fIt )
 807     (*fIt)->getGroupIds( support_groups, false );
 808
 809   numerotationMaillage(); // Renumerotation des mailles par entite
 810
 811   int dimension_maillage=getMeshDimension();
 812
 813   for (size_t i=0; i!=this->groupes.size(); ++i)
 814   {
 815     _groupe& grp = groupes[i];
 816     if ( grp.medGroup /*medGroupes[ i ]*/ )
 817       continue; // grp already converted into a MED group
 818
 819     bool isFieldSupport = ( support_groups.find( i ) != support_groups.end() );
 820     // si le groupe est vide, ou si le groupe n'est pas nommé : on passe au suivant
 821     if ( grp.empty() || ( grp.nom.empty() && !isFieldSupport )) {
 822       if ( !grp.nom.empty() ) {
 823         INFOS_MED("Skip group " << i << " " << grp.nom);
 824 //       } else {
 825 //         MESSAGE_MED("Skip group " << i << " <" << grp.nom << "> " << ( grp.empty() ? ": empty" : ""));
 826       }
 827       continue;
 828     }
 829
 830     // Build a set of mailles: sort mailles by type and exclude maille doubling
 831     bool isSelfIntersect = false;
 832     typedef set< set<_maille>::iterator, _mailleIteratorCompare > TMailleSet;
 833     TMailleSet mailleSet;
 834     if( grp.groupes.size() ) {// le groupe i contient des sous-maillages
 835       int nb_elem = 0;
 836       for( vector<int>::iterator j=grp.groupes.begin(); j!=grp.groupes.end(); ++j)
 837       {
 838         nb_elem += groupes[*j-1].mailles.size();
 839         _groupe::TMailleIter maIt=groupes[*j-1].mailles.begin();
 840         for( ; maIt!=groupes[*j-1].mailles.end(); ++maIt) {
 841 //           TMailleSet::const_iterator ma_it = mailleSet.find( *maIt );
 842 //           if ( ma_it != mailleSet.end() ) {
 843 //             MESSAGE_MED("EQUAL ELEMS: " << *ma_it << " AND " << *maIt);
 844 //           }
 845 //           else
 846             mailleSet.insert( *maIt );
 847         }
 848       }
 849       if ( nb_elem != mailleSet.size() ) { // Self intersecting compound group
 850         isSelfIntersect = true;
 851         INFOS_MED("Self intersecting group: " << i << " <" << grp.nom << ">"
 852               << ", mailleSet.size = " << mailleSet.size() << ", sum nb elems = " << nb_elem);
 853         for( vector<int>::iterator j=grp.groupes.begin(); j!=grp.groupes.end(); ++j)
 854           INFOS_MED(" in sub-group "<<  *j << " <" << groupes[*j-1].nom << "> "
 855                 << groupes[*j-1].mailles.size() << " mailles of type "
 856                 << groupes[*j-1].mailles[0]->geometricType);
 857         // if a group serve as a support to a field, then the _field is to be converted
 858         // into several MEDMEM::FIELDs, one per sub-group; if we already skipped some sub-groups,
 859         // we roll back the loop on groups to create MED groups for skipped ones
 860         if ( isFieldSupport ) {
 861           if ( grp.nom.empty() ) // clear group existing for the sake of field only
 862             grp.groupes.clear();
 863           std::set<int> sub_grps;
 864           for ( fIt = fields.begin(); fIt != fields.end(); ++fIt ) {
 865             _fieldBase * field = (*fIt);
 866             if ( field->_group_id == i ) {
 867               field->_group_id = -1; // -> a field by support
 868               field->getGroupIds( sub_grps, false );
 869             }
 870           }
 871           if ( (int)i > *sub_grps.begin() ) { // roll back
 872             support_groups.erase( i );
 873             support_groups.insert( sub_grps.begin(), sub_grps.end() );
 874             i = *sub_grps.begin() - 1;
 875             continue;
 876           }
 877         }
 878       }
 879     }
 880     else {
 881       _groupe::TMailleIter maIt=grp.mailles.begin();
 882       for(; maIt!=grp.mailles.end(); ++maIt)
 883         mailleSet.insert( *maIt );
 884       if ( grp.mailles.size() != mailleSet.size() )
 885         INFOS_MED( "Self intersecting group: " << i << " <" << grp.nom << ">"
 886                << ", mailleSet.size = " << mailleSet.size()
 887                << ", nb elems = " << grp.mailles.size());
 888     }
 889     // MEDMEM does not allow constituents of <mesh_dim>-2 and less dimension
 890     // but do not skip node group
 891     int group_min_dim = (**mailleSet.begin()).dimensionWithPoly();
 892     int group_max_dim = (**(--mailleSet.end())).dimensionWithPoly();
 893     if ( group_max_dim != 0 && group_min_dim <= dimension_maillage - 2  ) {
 894       MESSAGE_MED("Skip group: " << i << " <" << grp.nom << "> - too small dimension: "
 895               << group_min_dim);
 896       continue;
 897     }
 898
 899     // initialise groupe_entity a l'entite de la premiere maille du groupe
 900     medEntityMesh groupe_entity = (**mailleSet.rbegin()).getEntity(dimension_maillage);
 901     if ( hasMixedCells && group_min_dim > 0 )
 902       groupe_entity = MED_CELL;
 903
 904     // total nb of elements in mesh of groupe_entity
 905     int totalNbElements = 0;
 906     if ( groupe_entity == MED_NODE )
 907       totalNbElements = points.size() - nbMerged( MED_POINT1 );
 908     else {
 909       int entityDim = hasMixedCells ? -1 : group_min_dim;
 910       _maillageByDimIterator allMailles( *this, entityDim, true );
 911       while ( allMailles.nextType() )
 912         if ( allMailles.dim() > 0 )
 913           totalNbElements += allMailles.sizeWithoutMerged();
 914     }
 915     const bool isOnAll = ( mailleSet.size() == totalNbElements );
 916
 917     // if !isOnAll, build a map _maille::ordre() -> index in GROUP.getNumber(MED_ALL_ELEMENTS).
 918     // It is used while fields building.
 919     if ( !isOnAll || isSelfIntersect ) {
 920       TMailleSet::iterator maIt = mailleSet.begin();
 921       for ( int iMa = 0; maIt != mailleSet.end(); maIt++ )
 922         grp.relocMap.insert( make_pair( (*maIt)->ordre(), ++iMa ));
 923     }
 924     //Parcours des mailles (a partir de la deuxieme) pour compter les types geometriques
 925     int nb_geometric_types=1;
 926     TMailleSet::iterator j=mailleSet.begin();
 927     medGeometryElement geometrictype=(**j).geometricType;
 928     for ( ++j ; j!=mailleSet.end(); ++j )
 929     {
 930       //Compte nombre de types geometriques
 931       if ( (**j).geometricType != geometrictype ) // si on change de type geometrique
 932       {
 933         nb_geometric_types++;
 934         geometrictype=(**j).geometricType;
 935       }
 936     }
 937
 938     MED_EN::medGeometryElement * tab_types_geometriques = new MED_EN::medGeometryElement[nb_geometric_types];
 939     int * tab_index_types_geometriques = new int[nb_geometric_types+1];
 940     int * tab_numeros_elements = new int[mailleSet.size()];
 941     int * tab_nombres_elements = new int[nb_geometric_types];
 942
 943     //Remplit tableaux entree des methodes set
 944     int indice_mailles=0/*, maxOrdre = -1*/;
 945     j=mailleSet.begin();
 946     geometrictype=(**j).geometricType;
 947     tab_index_types_geometriques[0]=1;
 948     int indice_types_geometriques=1;
 949     tab_types_geometriques[0]=geometrictype;
 950     //parcours des mailles du groupe
 951     for (  ; j!=mailleSet.end(); ++j , ++indice_mailles)
 952     {
 953       const _maille& ma = **j;
 954       tab_numeros_elements[indice_mailles]= ma.ordre();
 955 //       if ( maxOrdre < tab_numeros_elements[indice_mailles] )
 956 //         maxOrdre = tab_numeros_elements[indice_mailles];
 957       if (ma.geometricType != geometrictype)
 958       {
 959         tab_index_types_geometriques[indice_types_geometriques]=indice_mailles+1;
 960         geometrictype=ma.geometricType;
 961         tab_types_geometriques[indice_types_geometriques]=geometrictype;
 962         ++indice_types_geometriques;
 963       }
 964     }
 965     tab_index_types_geometriques[indice_types_geometriques]=indice_mailles+1;
 966     for (int k=0; k != nb_geometric_types; ++k)
 967     {
 968       tab_nombres_elements[k] = tab_index_types_geometriques[k+1]-tab_index_types_geometriques[k];
 969     }
 970     //INFOS_MED( "MAX ORDRE in grp " << grp.nom << " entity " << groupe_entity << " : " << maxOrdre);
 971
 972     //Determination type entite du groupe
 973     vector <GROUP *> * vect_group;
 974     switch ( groupe_entity )
 975     {
 976     case MED_CELL :
 977       vect_group= & _groupCell;
 978       break;
 979     case MED_FACE :
 980       vect_group= & _groupFace;
 981       break;
 982     case MED_EDGE :
 983       vect_group= & _groupEdge;
 984       break;
 985     case MED_NODE :
 986       vect_group= & _groupNode;
 987       break;
 988     }
 989
 990     //Creation nouveau groupe MED
 991     GROUP * new_group = new GROUP();
 992     grp.medGroup = new_group;
 993     //Appel methodes set
 994     //new_group->setTotalNumberOfElements(mailleSet.size());
 995     new_group->setName(grp.nom);
 996     new_group->setMesh(_ptrMesh);
 997     new_group->setNumberOfGeometricType(nb_geometric_types);
 998     new_group->setGeometricType(tab_types_geometriques);
 999     new_group->setNumberOfElements(tab_nombres_elements);
1000     new_group->setNumber(tab_index_types_geometriques,tab_numeros_elements);
1001     new_group->setEntity(groupe_entity);
1002     new_group->setAll( isOnAll );
1003
1004     vect_group->push_back(new_group);
1005
1006     delete [] tab_types_geometriques;
1007     delete [] tab_index_types_geometriques;
1008     delete [] tab_numeros_elements;
1009     delete [] tab_nombres_elements;
1010   }
1011
1012   END_OF_MED(LOC);
1013 }
1014
1015 //=======================================================================
1016 //function : getFamilies
1017 //purpose  : create families like MESH::createFamilies() but preserves
1018 //           the order of elements in GROUPs defined by constituent families
1019 //           order. Call it after getGroups()
1020 //=======================================================================
1021
1022 // void _intermediateMED::getFamilies(std::vector<FAMILY *> & _famCell,
1023 //                                    std::vector<FAMILY *> & _famFace,
1024 //                                    std::vector<FAMILY *> & _famEdge,
1025 //                                    std::vector<FAMILY *> & _famNode, MESH * _ptrMesh)
1026 // {
1027 //   const char * LOC = "_intermediateMED::getFamilies() : ";
1028 //   BEGIN_OF_MED(LOC);
1029
1030 //   int nbElemFam = 0, nbNodeFam = 0;
1031 //   std::map< GROUP*, vector< FAMILY * > > grpFamsMap;
1032 //   int dimension_maillage=maillage.rbegin()->dimension();
1033
1034 //   std::set<_maille>::const_iterator i=maillage.begin(); // iterateurs sur les mailles
1035 //   std::set<_maille>::const_iterator j=maillage.begin();
1036
1037 //   do
1038 //   {
1039 //     // make a family containing mailles shared by the same set of groups
1040 //     std::list<unsigned>&  grList = i->groupes;  // to define the family end
1041 //     int           dimension = i->dimension();        // to define the entity end
1042 //     medGeometryElement type = i->geometricType;
1043 //     medEntityMesh    entity = i->getEntity( dimension_maillage );
1044
1045 //     std::vector<medGeometryElement> tab_types_geometriques;
1046 //     std::vector<int> tab_index_types_geometriques;
1047 //     std::vector<int> tab_nombres_elements;
1048 //     std::vector<int> tab_numeros_elements;
1049
1050 //     int iMa = 1, nbtype = 0;
1051 //     tab_types_geometriques.push_back( type );
1052 //     tab_index_types_geometriques.push_back( iMa );
1053
1054 //     // scan family cells and fill the tab that are needed by the create a MED FAMILY
1055 //     while (i != maillage.end() &&
1056 //            i->groupes == grList &&
1057 //            i->dimension() == dimension)
1058 //     {
1059 //       if (type != i->geometricType) // si changement de type geometrique
1060 //       {
1061 //         tab_index_types_geometriques.push_back(iMa);
1062 //         tab_nombres_elements.push_back(nbtype);
1063 //         nbtype=0;
1064 //         type=i->geometricType;
1065 //         tab_types_geometriques.push_back(type); // stocke le nouveau type geometrique rencontre
1066 //       }
1067 //       ++nbtype;
1068 //       ++iMa;
1069 //       ++i;
1070 //     }
1071 //     tab_index_types_geometriques.push_back(iMa);
1072 //     tab_nombres_elements.push_back(nbtype); // n'a pas été stocké dans la boucle
1073
1074 //     tab_numeros_elements.resize( iMa - 1 );
1075 //     for ( iMa = 0; j != i; j++, iMa++ )
1076 //       tab_numeros_elements[ iMa ] = j->ordre;
1077
1078 //     int id = ( entity == MED_NODE ? ++nbNodeFam : -(++nbElemFam) );
1079
1080 //     ostringstream name;
1081 //     name << "FAM_" << id;
1082
1083 //     // create a empty MED FAMILY and fill it with the tabs we constructed
1084 //     FAMILY* newFam = new FAMILY();
1085 //     newFam->setTotalNumberOfElements( iMa );
1086 //     newFam->setName( name.str() );
1087 //     newFam->setMesh( _ptrMesh );
1088 //     newFam->setNumberOfGeometricType( tab_types_geometriques.size() );
1089 //     newFam->setGeometricType( &tab_types_geometriques[0] ); // we know the tab is not empy
1090 //     newFam->setNumberOfElements( &tab_nombres_elements[0] );
1091 //     newFam->setNumber( &tab_index_types_geometriques[0], &tab_numeros_elements[0] );
1092 //     newFam->setEntity( entity );
1093 //     newFam->setAll( false );
1094 //     newFam->setIdentifier( id );
1095 //     newFam->setNumberOfGroups( grList.size() );
1096
1097 //     // Update links between families and groups
1098 //     if ( ! grList.empty() )
1099 //     {
1100 //       std::string * groupNames = new string[ grList.size() ];
1101 //       std::list<unsigned>::iterator g = grList.begin();
1102 //       for ( int i = 0; g != grList.end(); ++g, ++i ) {
1103 //         GROUP * medGROUP = getGroup( *g );
1104 //         groupNames[ i ] = medGROUP->getName();
1105 //         grpFamsMap[ medGROUP ].push_back( newFam );
1106 //       }
1107 //       newFam->setGroupsNames(groupNames);
1108 //     }
1109 //     // store newFam
1110 //     std::vector<FAMILY*>* families = 0;
1111 //     switch ( entity )
1112 //     {
1113 //     case MED_CELL :
1114 //       families = & _famCell; break;
1115 //     case MED_FACE :
1116 //       families = & _famFace; break;
1117 //     case MED_EDGE :
1118 //       families = & _famEdge; break;
1119 //     case MED_NODE :
1120 //       families = & _famNode; break;
1121 //     }
1122 //     if ( families )
1123 //       families->push_back( newFam );
1124
1125 //   } while ( i != maillage.end() );
1126
1127 //   // update references in groups
1128 //   std::map< GROUP*, vector< FAMILY * > >::iterator gf = grpFamsMap.begin();
1129 //   for ( ; gf != grpFamsMap.end(); ++gf )
1130 //   {
1131 //     gf->first->setNumberOfFamilies( gf->second.size() );
1132 //     gf->first->setFamilies( gf->second );
1133 //   }
1134 //}
1135
1136 //=======================================================================
1137 //function : getGroup
1138 //purpose  :
1139 //=======================================================================
1140
1141 // GROUP * _intermediateMED::getGroup( int i )
1142 // {
1143 //   if ( i <(int) medGroupes.size() )
1144 //     return medGroupes[ i ];
1145 //   throw MEDEXCEPTION
1146 //     (LOCALIZED(STRING("_intermediateMED::getGroup(): WRONG GROUP INDEX: ")
1147 //                << medGroupes.size() << " <= " << i ));
1148 // }
1149
1150 //=======================================================================
1151 //function : getFields
1152 //purpose  :
1153 //=======================================================================
1154
1155 void _intermediateMED::getFields(std::list< FIELD_* >& theFields)
1156 {
1157   const char * LOC = "_intermediateMED::getFields() : ";
1158   BEGIN_OF_MED(LOC);
1159
1160   std::list< _fieldBase* >::const_iterator fIt = fields.begin();
1161   for ( ; fIt != fields.end(); fIt++ )
1162   {
1163     const _fieldBase* fb = *fIt;
1164     list<pair< FIELD_*, int> >  ff = fb->getField(groupes);
1165     list<pair< FIELD_*, int> >::iterator f_sup = ff.begin();
1166     for (int j = 1 ; f_sup != ff.end(); f_sup++, ++j )
1167     {
1168       FIELD_    * f = f_sup->first;
1169       SUPPORT * sup = groupes[ f_sup->second ].medGroup;
1170       if ( !sup )
1171         throw MEDEXCEPTION
1172           (LOCALIZED(STRING(LOC) <<"_intermediateMED::getFields(), NULL field support: "
1173                      << " group index: " << fb->_group_id));
1174       int nb_elems = sup->getNumberOfElements( MED_ALL_ELEMENTS );
1175       if ( nb_elems != f->getNumberOfValues() )
1176         throw MEDEXCEPTION
1177           (LOCALIZED(STRING("_intermediateMED::getFields(), field support size (")
1178                      << nb_elems  << ") != NumberOfValues (" << f->getNumberOfValues()));
1179       theFields.push_back( f );
1180       if ( sup->getName().empty() ) {
1181         ostringstream name;
1182         name << "GRP_" << f->getName() << "_" << j;
1183         sup->setName( name.str() );
1184       }
1185       f->setSupport( sup );
1186       //f->setIterationNumber( j );
1187       f->setOrderNumber( j );
1188     }
1189   }
1190   END_OF_MED(LOC);
1191 }
1192
1193 //=======================================================================
1194 //function : ~_intermediateMED
1195 //purpose  : destructor
1196 //=======================================================================
1197
1198 _intermediateMED::~_intermediateMED()
1199 {
1200   MESSAGE_MED( "~_intermediateMED()");
1201   std::list< _fieldBase* >::const_iterator fIt = fields.begin();
1202   for ( ; fIt != fields.end(); fIt++ )
1203     delete *fIt;
1204 }
1205
1206 //=======================================================================
1207 //function : getGroupIds
1208 //purpose  : return ids of main and/or sub- groups
1209 //=======================================================================
1210
1211 void _fieldBase::getGroupIds( std::set<int> & ids, bool all ) const
1212 {
1213   if ( hasCommonSupport() )
1214     ids.insert( _group_id );
1215   if ( all || !hasCommonSupport() ) {
1216     vector< _sub_data >::const_iterator sub = _sub.begin();
1217     for ( ; sub != _sub.end(); ++sub )
1218       ids.insert( sub->_supp_id );
1219   }
1220 }
1221
1222 //=======================================================================
1223 //function : hasSameComponentsBySupport
1224 //purpose  :
1225 //=======================================================================
1226
1227 bool _fieldBase::hasSameComponentsBySupport() const
1228 {
1229   vector< _sub_data >::const_iterator sub_data = _sub.begin();
1230   const _sub_data& first_sub_data = *sub_data;
1231   for ( ++sub_data ; sub_data != _sub.end(); ++sub_data )
1232   {
1233     if ( first_sub_data._comp_names != sub_data->_comp_names )
1234       return false; // diff names of components
1235
1236     if ( first_sub_data._nb_gauss != sub_data->_nb_gauss )
1237       return false; // diff nb of gauss points
1238   }
1239   return true;
1240 }
1241
1242 //=======================================================================
1243 //
1244 //function : mergeNodes
1245 //
1246 //=======================================================================
1247 namespace {
1248
1249 struct __NOEUD
1250 {
1251   typedef map<int,_noeud>::iterator Inoeud;
1252   Inoeud i_noeud;
1253
1254   __NOEUD() {}
1255   __NOEUD(Inoeud n): i_noeud(n) {}
1256   const double & operator[] (int i) const { return i_noeud->second.coord[i]; }
1257   double         operator[] (int i)       { return i_noeud->second.coord[i]; }
1258   int dim() const { return i_noeud->second.coord.size(); }
1259   int& num() { return i_noeud->second.number; }
1260   int id() const { return i_noeud->first; }
1261   bool isMerged() const { return i_noeud->second.number < 1; }
1262 };
1263 //-----------------------------------------------------------------------
1264 double DistanceL2(const __NOEUD &A,const __NOEUD &B)
1265 {
1266   if ( B.isMerged() ) return DBL_MAX;
1267   const double* cooA = &A[0], *cooB = &B[0], *cooEnd = cooA + A.dim();
1268   double dist, somme=0;
1269   while ( cooA < cooEnd ) {
1270     dist=((*cooA++) - (*cooB++));
1271     somme+=dist*dist;
1272   }
1273   return sqrt(somme);
1274 }
1275 //-----------------------------------------------------------------------
1276 struct __NUAGENOEUD
1277 {
1278   vector< __NOEUD > nodes;
1279   __NUAGENOEUD(_intermediateMED& imed);
1280   __NOEUD & operator [] (int i) { return nodes[i]; }
1281   int size() const { return nodes.size(); }
1282 };
1283 //-----------------------------------------------------------------------
1284 __NUAGENOEUD::__NUAGENOEUD(_intermediateMED& imed)
1285 {
1286   nodes.resize( imed.points.size() );
1287   map<int,_noeud>::iterator i_noeud = imed.points.begin();
1288   for( int i = 0; i_noeud!=imed.points.end(); ++i_noeud, ++i ) {
1289     i_noeud->second.number = i+1;
1290     nodes[ i ] = i_noeud;
1291   }
1292 }
1293 //-----------------------------------------------------------------------
1294 // mergeNodes
1295 template<int DIM> int mergeNodes(double            tolerance,
1296                                  _intermediateMED& imed,
1297                                  vector< int > &   /*newNodeIDs*/)
1298 {
1299   typedef dTree<__NOEUD,__NUAGENOEUD,DIM > DTree;
1300   __NUAGENOEUD aNUAGENOEUD( imed );
1301   DTree tree( &aNUAGENOEUD );
1302
1303 //   int maxID = imed.points.rbegin()->first;
1304 //   newNodeIDs.resize( maxID + 1, 0 );
1305
1306   int num = 0, nbRemoved = 0;
1307   int nbNodes = aNUAGENOEUD.size();
1308   for ( int i = 0; i < nbNodes; ++i )
1309   {
1310     __NOEUD & node = aNUAGENOEUD[i];
1311     int & number = node.num();
1312     if ( number < 0 )
1313       continue; // already merged
1314     number = ++num;
1315
1316     list<int> closeNumbers;
1317     int nbCoinsident = tree.get_all_close( node, tolerance, closeNumbers );
1318     if ( nbCoinsident > 1 ) // not only node it-self found
1319     {
1320       nbRemoved += nbCoinsident-1; // one of coincident nodes remains
1321       int id = node.id();
1322       list<int>::iterator n = closeNumbers.begin(), nEnd = closeNumbers.end();
1323       while ( n != nEnd ) {
1324         __NOEUD & coincNode = aNUAGENOEUD[ *n++ ];
1325         int coincID = coincNode.id();
1326         if ( coincID != id ) {
1327           coincNode.num() = -number;
1328           //newNodeIDs[ coincID ] = id;
1329         }
1330       }
1331     }
1332   }
1333   return nbRemoved;
1334 }
1335 //-----------------------------------------------------------------------
1336 // wrapper of _maille used after merging nodes to find equal mailles
1337 struct _mailleWrap
1338 {
1339   const _maille* _ma;
1340
1341   _mailleWrap(const _maille* ma=0): _ma(ma) { if (_ma) _ma->init(); }
1342   ~_mailleWrap()                            { if (_ma) _ma->init(); }
1343
1344   bool operator < (const _mailleWrap& mw) const
1345   {
1346     size_t l=_ma->sommets.size();
1347     const int* v1 = getSortedNodeNums(_ma);
1348     const int* v2 = getSortedNodeNums(mw._ma);
1349     for ( const int* vEnd = v1 + l; v1 < vEnd; ++v1, ++v2 )
1350       if(*v1 != *v2)
1351         return *v1 < *v2;
1352     return false; // cas d'égalité
1353   }
1354   static const int* getSortedNodeNums(const _maille* ma)
1355   {
1356     if ( !ma->sortedNodeIDs ) {
1357       size_t l=ma->sommets.size();
1358       ma->sortedNodeIDs = new int[ l ];
1359       for (size_t i=0; i!=l; ++i)
1360         ma->sortedNodeIDs[i]=ma->nodeNum(i);
1361       std::sort( ma->sortedNodeIDs, ma->sortedNodeIDs + l );
1362     }
1363     return ma->sortedNodeIDs;
1364   }
1365 };
1366
1367 } // namespace
1368
1369 //=======================================================================
1370 //function : mergeNodesAndElements
1371 //purpose  : optionally merge nodes and elements
1372 //=======================================================================
1373
1374 void _intermediateMED::mergeNodesAndElements(double tolerance)
1375 {
1376   //const char * LOC = "_intermediateMED::mergeNodesAndElements(): ";
1377   vector< int > newNodeIDs; // for i-th node id, id to replace with
1378
1379   int nbRemovedNodes = 0;
1380   const int spaceDimension=points.begin()->second.coord.size();
1381   if ( spaceDimension == 3 )
1382     nbRemovedNodes = mergeNodes<3>( tolerance, *this, newNodeIDs );
1383   else if ( spaceDimension == 2 )
1384     nbRemovedNodes = mergeNodes<2>( tolerance, *this, newNodeIDs );
1385
1386   myNodesNumerated = true;
1387
1388   if ( nbRemovedNodes == 0 )
1389     return;
1390
1391   // all commented code here was used to keep initial numeration but it was too slow
1392   //numerotationMaillage();
1393   treatGroupes();
1394
1395   nbRemovedByType[ MED_NONE   ] = nbRemovedNodes;
1396   nbRemovedByType[ MED_POINT1 ] = nbRemovedNodes;
1397
1398   bool hasPointMailles = false;
1399   _maillageByDimIterator entityMailles( *this, 0 );
1400   if ( const set<_maille> * points = entityMailles.nextType() ) {
1401     hasPointMailles = true;
1402     set<_maille>::const_iterator i, iEnd = points->end();
1403     for ( i = points->begin(); i != iEnd; ++i )
1404       i->setOrdre( i->sommets[0]->second.number ); // is merged if point is merged
1405   }
1406   const bool skipFirstType = ( hasPointMailles && hasMixedCells );
1407   // loop on entities
1408   for ( int dim = 1; dim <= 3; ++dim )
1409   {
1410     int iterDim = hasMixedCells ? -1 : dim;
1411     //int nbRemovedInEntity = 0;
1412
1413     // count total nb of mailles in entity
1414 //     int nbMailles = 0;
1415 //     entityMailles.init( iterDim, true );
1416 //     if ( skipFirstType ) entityMailles.nextType(); // merged POINT1's are same as nodes
1417 //     while ( const set<_maille> * typeMailles = entityMailles.nextType() )
1418 //       nbMailles += typeMailles->size();
1419
1420     // for each maille number, count shift due to removed mailles with lower numbers
1421     //vector<int> shift( nbMailles+1, 0 );
1422
1423     // merge and numerate mailles
1424     int num = 0;
1425     entityMailles.init( iterDim, true );
1426     if ( skipFirstType ) entityMailles.nextType(); // merged POINT1's are same as nodes
1427     while ( const set<_maille> * typeMailles = entityMailles.nextType() )
1428     {
1429       int nbRemoved = 0;
1430       set<_mailleWrap> maNodeNumSorted;
1431       pair< set<_mailleWrap>::const_iterator, bool > mw_isunique;
1432
1433       set<_maille>::const_iterator ma = typeMailles->begin(), maEnd = typeMailles->end();
1434       while ( ma != maEnd )
1435       {
1436         const _maille & m = *ma++;
1437         mw_isunique = maNodeNumSorted.insert( &m );
1438         if ( mw_isunique.second ) {
1439           m.setOrdre( ++num );
1440         }
1441         else {
1442           const _maille* equalMa = mw_isunique.first->_ma;
1443           //unsigned ordreToRemove = m.ordre();
1444           m.setMergedOrdre( equalMa->ordre() );
1445           nbRemoved++;
1446 //           while ( ordreToRemove <= nbMailles )
1447 //             shift[ ordreToRemove++ ]++;
1448         }
1449       }
1450
1451       if ( nbRemoved ) {
1452         nbRemovedByType[ entityMailles.type() ] = nbRemoved;
1453         //nbRemovedInEntity += nbRemoved;
1454       }
1455
1456       // update maille ordre
1457 //       if ( nbRemovedInEntity ) {
1458 //         for ( ma = typeMailles->begin(); ma != maEnd; ++ma ) {
1459 //           unsigned newOrdre = ma->ordre() - shift[ ma->ordre() ];
1460 //           if ( ma->isMerged() )
1461 //             ma->setMergedOrdre( newOrdre );
1462 //           else
1463 //             ma->setOrdre( newOrdre );
1464 //         }
1465 //       }
1466     }
1467   }
1468
1469   myMaillesNumerated = true;
1470
1471 }
1472
1473 //=======================================================================
1474 //function : getMeshDimension
1475 //purpose  : count mesh dimension
1476 //=======================================================================
1477
1478 int _intermediateMED::getMeshDimension() const
1479 {
1480   int dim = 0;
1481   _maillageByDimIterator allMailles( *this, -1, true );
1482   while ( allMailles.nextType() )
1483     dim = allMailles.dim();
1484   return dim;
1485 }
1486
1487 }
1488
1489 /////