src/MEDMEM/MEDMEM_DriverTools.cxx

   1 #include "MEDMEM_DriverTools.hxx"
   2 #include "MEDMEM_STRING.hxx"
   3 #include "MEDMEM_Exception.hxx"
   4 #include "MEDMEM_Mesh.hxx"
   5 #include "MEDMEM_Group.hxx"
   6 #include <iomanip>
   7 #include <algorithm>
   8
   9 using namespace std;
  10 using namespace MEDMEM;
  11 using namespace MED_EN;
  12
  13 #define DUMP_LINES_LIMIT 20
  14
  15 // Cet opérateur permet d'ordonner les mailles dans un set suivant l'ordre requis par MED
  16 bool _maille::operator < (const _maille& ma) const
  17 {
  18     // si le type géométrique differe, la comparaison est basée dessus
  19     // sinon on se base sur une comparaison des numéros de sommets
  20     if(geometricType==ma.geometricType)
  21     {
  22         // construction de deux vecteur temporaire contenant les numeros de sommets
  23         // pour faire le tri et les comparaisons
  24         size_t l=sommets.size();
  25         std::vector<int> v1(l);
  26         std::vector<int> v2(l);
  27         for (unsigned int i=0; i!=l; ++i)
  28           {
  29             v1[i]=sommets[i]->second.number;
  30             v2[i]=ma.sommets[i]->second.number;
  31           }
  32         std::sort(v1.begin(), v1.end());
  33         std::sort(v2.begin(), v2.end());
  34         for(std::vector<int>::const_iterator i1=v1.begin(), i2=v2.begin(); i1!=v1.end(); ++i1, ++i2)
  35             if(*i1 != *i2)
  36                 return *i1 < *i2;
  37         return false; // cas d'égalité
  38     }
  39     else
  40         return geometricType<ma.geometricType;
  41 };
  42 _link _maille::link(int i) const
  43 {
  44   ASSERT ( i >= 0 && i < sommets.size() );
  45   int i2 = ( i + 1 == sommets.size() ) ? 0 : i + 1;
  46   if ( reverse )
  47     return make_pair( sommets[i2]->first, sommets[i]->first );
  48   else
  49     return make_pair( sommets[i]->first, sommets[i2]->first );
  50 }
  51
  52 // retourne l'entité d'une maille en fonction de la dimension du maillage.
  53 MED_EN::medEntityMesh _maille::getEntity(const int meshDimension) const throw (MEDEXCEPTION)
  54 {
  55   const char * LOC = "_maille::getEntity(const int meshDimension)" ;
  56 //   BEGIN_OF(LOC);
  57
  58   int mailleDimension = this->dimension();
  59   medEntityMesh entity;
  60   if (mailleDimension == meshDimension)
  61     entity = MED_CELL;
  62   else
  63     switch (mailleDimension)
  64       {
  65       case 0 :
  66         entity = MED_NODE;
  67         break;
  68       case 1 :
  69         entity = MED_EDGE;
  70         break;
  71       case 2 :
  72         entity = MED_FACE;
  73         break;
  74       default :
  75         throw MEDEXCEPTION(LOCALIZED(STRING(LOC) << "Impossible de determiner l'entite de la maille."));
  76       }
  77 return entity;
  78
  79 //END_OF(LOC);
  80 };
  81
  82 std::ostream& MEDMEM::operator << (std::ostream& os, const _maille& ma)
  83 {
  84     os << "maille " << ma.ordre << " (" << ma.geometricType << ") : < ";
  85     std::vector< std::map<int,_noeud>::iterator >::const_iterator i=ma.sommets.begin();
  86     os << (*i++)->second.number ;
  87     for( ; i!=ma.sommets.end(); ++i)
  88         os << ", " << (*i)->second.number;
  89     os << " >";
  90     return os;
  91 }
  92
  93 std::ostream& MEDMEM::operator << (std::ostream& os, const _groupe& gr)
  94 {
  95     os << "--- Groupe " << gr.nom << " --- " << std::endl ;
  96     os << " -> liste des sous-groupes : ";
  97     for( std::vector<int>::const_iterator i=gr.groupes.begin(); i!=gr.groupes.end(); ++i)
  98         os << *i << " ";
  99     os << std::endl << " -> liste des "<< gr.mailles.size() << " mailles : " << std::endl;
 100     _groupe::mailleIter i=gr.mailles.begin();
 101     int l;
 102     for(l = 0; l < DUMP_LINES_LIMIT && i!=gr.mailles.end(); i++, l++)
 103         os << setw(3) << l+1 << " " << *(*i) << std::endl;
 104     if ( l == DUMP_LINES_LIMIT )
 105       os << "   ... skip " << gr.mailles.size() - l << " mailles" << endl;
 106     os << " relocMap, size=" << gr.relocMap.size() << endl;
 107     map<const _maille*,int>::const_iterator it = gr.relocMap.begin();
 108     for ( l = 0; l < DUMP_LINES_LIMIT && it != gr.relocMap.end(); ++it, ++l )
 109       os << " (" << it->first << "," << it->second << ")";
 110     if ( gr.relocMap.size() > 0 )
 111       os << endl;
 112     return os;
 113 }
 114
 115 std::ostream& MEDMEM::operator << (std::ostream& os, const _noeud& no)
 116 {
 117     os << "noeud " << no.number << " : < ";
 118     std::vector<double>::const_iterator i=no.coord.begin();
 119     os << *i++ ;
 120     for( ; i!=no.coord.end(); ++i)
 121         os << ", " << *i;
 122     os << " >";
 123     return os;
 124 }
 125
 126 void MEDMEM::_fieldBase::dump(std::ostream& os) const
 127 {
 128   os << "field " << "<" << _name << ">" << endl <<
 129     "  nb sub: " << _sub.size() << endl <<
 130     "  group index: " << _group_id << endl <<
 131     "  type: " << _type << endl;
 132   os << "  subcomponents:" << endl;
 133   vector< _sub_data >::const_iterator sub_data = _sub.begin();
 134   for ( int i_sub = 1; sub_data != _sub.end(); ++sub_data, i_sub ) {
 135     os << "    group index: " << sub_data->_supp_id <<
 136       ", " << sub_data->nbComponents() << " comp names: ";
 137     for ( int i_comp = 0; i_comp < sub_data->nbComponents(); ++i_comp )
 138       os << " |" << sub_data->_comp_names[ i_comp ] << "|";
 139     os << endl;
 140   }
 141 }
 142
 143 std::ostream& MEDMEM::operator << (std::ostream& os, const _fieldBase * f)
 144 {
 145   f->dump( os );
 146   return os;
 147 }
 148
 149 std::ostream& MEDMEM::operator << (std::ostream& os, const _intermediateMED& mi)
 150 {
 151     os << "Set des " << mi.maillage.size() << " mailles : " << std::endl;
 152     std::set<_maille>::const_iterator i=mi.maillage.begin();
 153     int l;
 154     for( l = 0; l < DUMP_LINES_LIMIT && i!=mi.maillage.end(); ++i, ++l)
 155         os << setw(3) << l+1 << " " <<*i << std::endl;
 156     if ( l == DUMP_LINES_LIMIT )
 157       os << "   ... skip " << mi.maillage.size() - l << " mailles" << endl;
 158
 159     os << std::endl << "Vector des " << mi.groupes.size() << " groupes : " << std::endl;
 160     for (unsigned int k=0; k!=mi.groupes.size(); k++)
 161       os << std::setw(3) << k << " " << mi.groupes[k] << std::endl;
 162
 163     os << std::endl << "map des " << mi.points.size() << " noeuds : " << std::endl;
 164     std::map<int,_noeud>::const_iterator j=mi.points.begin();
 165     for( l = 0; l < DUMP_LINES_LIMIT && j!=mi.points.end(); ++j, ++l)
 166         os << j->second << std::endl;
 167     if ( l == DUMP_LINES_LIMIT )
 168       os << "   ... skip " << mi.points.size() - l << " noeuds" << endl;
 169
 170     os << endl << mi.fields.size() << " fields:" << endl;
 171     std::list<_fieldBase* >::const_iterator fIt = mi.fields.begin();
 172     for ( l = 0; fIt != mi.fields.end(); ++fIt, ++l )
 173       os << " - " << l+1 << " " << *fIt << endl;
 174
 175     return os;
 176 }
 177
 178
 179 //=======================================================================
 180 //function : treatGroupes
 181 //purpose  : detect groupes of mixed dimension and erase groupes that
 182 //           won't be converted
 183 //=======================================================================
 184
 185 void _intermediateMED::treatGroupes()
 186 {
 187   const char * LOC = "_intermediateMED::treatGroupes() : ";
 188   BEGIN_OF(LOC);
 189
 190   // --------------------
 191   // erase useless group
 192   // --------------------
 193
 194   // decrease hierarchical depth of subgroups
 195   vector<int>::iterator j;
 196   for (unsigned int i=0; i!=this->groupes.size(); ++i)
 197   {
 198     _groupe& grp = groupes[i];
 199     //INFOS( i << " " << grp.nom );
 200     j = grp.groupes.begin();
 201     while( j!=grp.groupes.end() ) {
 202       int grpInd = *j-1;
 203       if ( grpInd < 0 || grpInd >= groupes.size() ) {
 204         throw MEDEXCEPTION(LOCALIZED(STRING(LOC) << "Bad subgroup index: " << grpInd <<
 205                                      ", in " << i << " groupe.nom=" << grp.nom));
 206       }
 207       _groupe & sub_grp = groupes[ grpInd ];
 208       if ( !sub_grp.groupes.empty() ) {
 209         MESSAGE("High hierarchical depth of subgroups in group " << i );
 210         *j = sub_grp.groupes[0]; // replace j with its 1st subgroup
 211         // push back the rest subs
 212         for ( int k = 1; k < sub_grp.groupes.size(); ++k )
 213           grp.groupes.push_back( sub_grp.groupes[ k ]);
 214         // vector maybe is reallocated: restart iterator
 215         j = grp.groupes.begin();
 216       }
 217       else
 218         j++;
 219     }
 220     // remove empty sub-groupes
 221     j = grp.groupes.begin();
 222     while ( j!=grp.groupes.end() ) {
 223       if ( groupes[*j-1].empty() ) {
 224         grp.groupes.erase( j );
 225         j = grp.groupes.begin();
 226       }
 227       else
 228         j++;
 229     }
 230   }
 231   // get indices of groups that are field support -
 232   // do not erase them and their subgroups
 233   std::set<int> groups2convert;
 234   std::list< _fieldBase* >::const_iterator fIt = fields.begin();
 235   for ( ; fIt != fields.end(); ++fIt )
 236     (*fIt)->getGroupIds( groups2convert, true );
 237
 238   // keep named groups and their subgroups
 239   for (unsigned int i=0; i!=this->groupes.size(); ++i)
 240   {
 241     _groupe& grp = groupes[i];
 242     if ( grp.empty() || grp.nom.empty() )
 243       continue;
 244     groups2convert.insert( i );
 245     for( j = grp.groupes.begin(); j!=grp.groupes.end(); ++j)
 246       groups2convert.insert( *j-1 );
 247   }
 248   // erase groups that are not in groups2convert
 249   for (unsigned int i=0; i!=this->groupes.size(); ++i)
 250   {
 251     if ( groups2convert.find( i ) == groups2convert.end() ) {
 252       _groupe& grp = groupes[i];
 253       grp.mailles.clear();
 254       grp.groupes.clear();
 255       MESSAGE( "Erase " << i << "-th group " << grp.nom );
 256     }
 257   }
 258
 259   // ---------------------------------------------------
 260   // define if there are groups with mixed entity types
 261   // ---------------------------------------------------
 262
 263   hasMixedCells = false;
 264   for (unsigned int i=0; i!=this->groupes.size(); ++i)
 265   {
 266     _groupe& grp = groupes[i];
 267     if ( grp.groupes.empty() )
 268       continue;
 269
 270     // check if sub-groups have different dimension
 271     j = grp.groupes.begin();
 272     int dim = groupes[*j-1].mailles[0]->dimension();
 273     for( j++; !hasMixedCells && j!=grp.groupes.end(); ++j)
 274       hasMixedCells = ( dim != groupes[*j-1].mailles[0]->dimension() );
 275   }
 276
 277   if ( hasMixedCells )
 278     INFOS( "There will be groups of mixed dimention" );
 279   END_OF(LOC);
 280 }
 281
 282 void _intermediateMED::numerotationMaillage()
 283 {
 284   const char * LOC = "_intermediateMED::numerotationMaillage() : ";
 285   BEGIN_OF(LOC);
 286
 287   treatGroupes();
 288
 289   // numerotation mailles of type MED_POINT1 by node number
 290   std::set<_maille>::iterator i=maillage.begin();
 291   if ( i->geometricType == MED_POINT1 ) {
 292     numerotationPoints();
 293     while ( i!=maillage.end() && i->geometricType == MED_POINT1 ) {
 294       i->ordre = i->sommets[0]->second.number;
 295       i++;
 296     }
 297   }
 298   // numerotation des mailles par entité
 299     int i_maille=0;
 300     int dimension=i->dimension();
 301     for( ; i!=maillage.end(); ++i)
 302     {
 303         if ( !hasMixedCells && dimension!=i->dimension() ) // on change d'entite
 304         {
 305           MESSAGE( "NB dim " << dimension << " entities: " << i_maille);
 306             dimension=i->dimension();
 307             i_maille=0;
 308         }
 309         (*i).ordre=++i_maille;
 310     }
 311   END_OF(LOC);
 312 }
 313
 314 void _intermediateMED::numerotationPoints()
 315 {
 316     // Fonction de renumerotation des noeuds (necessaire quand il y a des trous dans la numerotation.
 317     int i_noeud=0;
 318     for( std::map<int,_noeud>::const_iterator i=points.begin(); i!=points.end(); ++i)
 319         i->second.number = ++i_noeud ;
 320 }
 321
 322
 323 /*!
 324  * \if developper
 325  * create a MED COORDINATE from the intermediate structure.
 326  * \endif
 327  */
 328 COORDINATE *
 329 _intermediateMED::getCoordinate(const string & coordinateSystem)
 330 {
 331     const medModeSwitch mode=MED_FULL_INTERLACE;
 332     int spaceDimension=points.begin()->second.coord.size();
 333     int numberOfNodes=points.size();
 334
 335
 336     // creation du tableau des coordonnees en mode MED_FULL_INTERLACE
 337     double * coord = new double[spaceDimension*numberOfNodes];
 338     int k=0;
 339     for( std::map<int,_noeud>::const_iterator i=points.begin(); i!=points.end(); ++i, ++k)
 340         std::copy(i->second.coord.begin(), i->second.coord.end(), coord+k*spaceDimension);
 341
 342     // creation de l'objet COORDINATE
 343     COORDINATE * coordinate = new COORDINATE(spaceDimension, numberOfNodes, mode);
 344     coordinate->setCoordinates(mode,coord);
 345     delete [] coord;
 346     coordinate->setCoordinatesSystem(coordinateSystem);
 347     return coordinate;
 348 }
 349
 350
 351   /*!
 352    * \if developper
 353    * create a MED CONNECTIVITY from the intermediate structure.
 354    * for memory optimisation, clear the intermediate structure (the "maillage" set )
 355    * \endif
 356    */
 357 CONNECTIVITY * _intermediateMED::getConnectivity()
 358 {
 359     const char * LOC = "_intermediateMED::getConnectivity() : ";
 360     BEGIN_OF(LOC);
 361     int numberOfNodes=points.size(); // number of nodes in the mesh
 362     int numberOfTypes=0;
 363     medEntityMesh entity;
 364     medGeometryElement * types=NULL; // pointeurs pour allouer les structures MED
 365     int * count=NULL;
 366     int * connectivity=NULL;
 367     CONNECTIVITY *Connectivity, *Constituent;
 368     int dimension_maillage_moin_2=maillage.rbegin()->dimension() - 2;
 369
 370     medGeometryElement type=0; // variables de travail
 371     int nbtype=0;
 372     int dimension=0;
 373     bool first = true;
 374
 375     std::vector<medGeometryElement> vtype; // tableau de travail : stockage des types geometriques pour UNE entite
 376     std::vector<int> vcount; // tableau de travail : nombre d'elements pour chaque type geometrique de vtype
 377
 378     std::set<_maille>::const_iterator i, j; // iterateurs sur les mailles
 379
 380     // skip nodes and elements of <dimension_maillage - 2> or less dimension
 381     for ( i = maillage.begin(); i != maillage.end(); ++i )
 382       //if ( i->geometricType != MED_POINT1 && i->dimension() > dimension_maillage_moin_2 )
 383         break;
 384     j = i;
 385
 386     // renumerote les points de 1 a n (pour le cas ou certains points ne sont pas presents dans le maillage d'origine)
 387     numerotationPoints();
 388
 389     do
 390     {
 391         // boucle sur les entites de la structure maillage :
 392         //   - on parcourt une premiere fois avec i les mailles d'une entite pour récupérer
 393         //     des infos sur les types geometriques, leur nombre et le nombre d'elements.
 394         //   - on alloue la connectivite
 395         //   - on parcourt une deuxieme fois avec j pour lire les noeuds.
 396
 397
 398         type=i->geometricType; // init boucle for
 399         dimension=i->dimension();
 400         nbtype=0;
 401         vtype.push_back(type);
 402         // Boucle sur i de parcours des mailles d'une entite
 403         // Une entite se termine lorsqu'on atteint la fin de maillage ou lorsque la dimension des mailles change
 404         for( ; i!=maillage.end() && ( hasMixedCells || dimension==i->dimension()) ; ++i)
 405         {
 406             if (type != i->geometricType) // si changement de type geometrique
 407             {
 408                 vcount.push_back(nbtype); // stocke le nombre d'elements de type nbtype
 409                 nbtype=0;
 410                 type=i->geometricType;
 411                 vtype.push_back(type); // stocke le nouveau type geometrique rencontre
 412             }
 413
 414             ++nbtype;
 415         }
 416         vcount.push_back(dimension ? nbtype : numberOfNodes); // n'a pas été stocké dans la boucle
 417
 418         // Pour l'entite qu'on vient de parcourir, allocation des tableau et creation de la connectivite
 419 //      cout << "Dimension = " << dimension << endl;
 420 //      cout << "Nombre de type geometriques : " << vtype.size() << endl;
 421 //      for (unsigned k=0; k!=vtype.size(); ++k )
 422 //          cout << "  -> " << vtype[k] << " : " << vcount[k] << endl;
 423
 424         numberOfTypes=vtype.size(); // nombre de types de l'entite
 425
 426         if ( i==maillage.end() ) // cas de l'entite de plus haut niveau
 427             entity=MED_CELL;
 428         else if (dimension==2 )
 429             entity=MED_FACE;
 430         else if (dimension==1 )
 431             entity=MED_EDGE;
 432         else if (dimension==0 )
 433             entity=MED_NODE;
 434
 435         Connectivity = new CONNECTIVITY(numberOfTypes,entity);
 436         Connectivity->setEntityDimension(dimension);
 437         Connectivity->setNumberOfNodes(numberOfNodes);
 438
 439         types = new medGeometryElement[numberOfTypes];
 440         std::copy(vtype.begin(),vtype.end(),types);
 441         Connectivity->setGeometricTypes(types,entity);
 442         delete [] types;
 443
 444         count = new int[numberOfTypes+1];
 445         count[0]=1;
 446         for (unsigned int k=0; k!=vcount.size(); ++k )
 447           count[k+1]=count[k]+vcount[k];
 448         Connectivity->setCount (count, entity);
 449         SCRUTE( entity );
 450         SCRUTE(numberOfTypes);
 451         SCRUTE(count[numberOfTypes]-1);
 452         delete [] count;
 453
 454         for (int k=0; k!=numberOfTypes; ++k )
 455           {
 456             // pour chaque type géometrique k, copie des sommets dans connectivity et set dans Connectivity
 457             int nbSommetsParMaille = j->sommets.size();
 458             int n, nbSommets = vcount[k] * j->sommets.size();
 459             connectivity = new int[ nbSommets ];
 460             for (int l=0; l!=vcount[k]; ++l)
 461             {
 462                 if ( entity==MED_NODE )
 463                   connectivity[l] = l+1;
 464                 else
 465                 {
 466                   for ( n=0; n != nbSommetsParMaille; ++n) {
 467                     connectivity[nbSommetsParMaille*l+n] =
 468                       j->sommets[ j->reverse ? nbSommetsParMaille-n-1 : n ]->second.number;
 469                   }
 470                 // DO NOT ERASE, maillage will be used while fields construction
 471                 //maillage.erase(j);    ; // dangereux, mais optimise la mémoire consommée!
 472                   j++;
 473                 }
 474             }
 475
 476             Connectivity->setNodal  (connectivity, entity, vtype[k]);
 477             delete [] connectivity;
 478           }
 479
 480         if ( entity==MED_NODE )
 481           j = i;
 482         else if (i!=j)
 483             throw MEDEXCEPTION(LOCALIZED(STRING(LOC) << "Erreur de lecture des mailles ."));
 484
 485         if ( ! first)
 486             Connectivity->setConstituent (Constituent);
 487         else
 488             first = false;
 489         Constituent = Connectivity; // stocke Connectivity pour utilisation dans setConstituent lors de la boucle suivante
 490
 491         vtype.clear();
 492         vcount.clear();
 493
 494     }
 495     while ( i != maillage.end() );
 496
 497     END_OF(LOC);
 498     return Connectivity;
 499 }
 500
 501
 502   /*!
 503    * \if developper
 504    * fill the arguments vector of groups from the intermediate structure.
 505    * This function must be called before getConnectivity()
 506    * \endif
 507    */
 508 void
 509 _intermediateMED::getGroups(vector<GROUP *> & _groupCell,
 510                             vector<GROUP *> & _groupFace,
 511                             vector<GROUP *> & _groupEdge,
 512                             vector<GROUP *> & _groupNode, MESH * _ptrMesh)
 513 {
 514   const char * LOC = "_intermediateMED::getGroups() : ";
 515   BEGIN_OF(LOC);
 516
 517   medGroupes.resize( groupes.size(), 0 );
 518   if (maillage.size() == 0) {
 519     INFOS( "Erreur : il n'y a plus de mailles");
 520     return;
 521   }
 522
 523   // get indices of groups that are field support - do not skip them
 524   std::set<int> support_groups;
 525   std::list< _fieldBase* >::const_iterator fIt = fields.begin();
 526   for ( ; fIt != fields.end(); ++fIt )
 527     (*fIt)->getGroupIds( support_groups, false );
 528
 529   numerotationMaillage(); // Renumerotation des mailles par entite
 530
 531   int dimension_maillage=maillage.rbegin()->dimension();
 532
 533   for (unsigned int i=0; i!=this->groupes.size(); ++i)
 534   {
 535     _groupe& grp = groupes[i];
 536     if ( medGroupes[ i ] )
 537       continue; // grp already converted into a MED group
 538
 539     bool isFieldSupport = ( support_groups.find( i ) != support_groups.end() );
 540     // si le groupe est vide, ou si le groupe n'est pas nommé : on passe au suivant
 541     if ( grp.empty() || ( grp.nom.empty() && !isFieldSupport )) {
 542       if ( !grp.nom.empty() ) {
 543         INFOS("Skip group " << i << grp.nom);
 544 //       } else {
 545 //         MESSAGE("Skip group " << i << " <" << grp.nom << "> " << ( grp.empty() ? ": empty" : ""));
 546       }
 547       continue;
 548     }
 549
 550     // Build a set of mailles: sort mailles by type and exclude maille doubling
 551     typedef set< set<_maille>::iterator, _mailleIteratorCompare > TMailleSet;
 552     TMailleSet mailleSet;
 553     if( grp.groupes.size() ) {// le groupe i contient des sous-maillages
 554       int nb_elem = 0;
 555       for( vector<int>::iterator j=grp.groupes.begin(); j!=grp.groupes.end(); ++j)
 556       {
 557         nb_elem += groupes[*j-1].mailles.size();
 558         _groupe::mailleIter maIt=groupes[*j-1].mailles.begin();
 559         for( ; maIt!=groupes[*j-1].mailles.end(); ++maIt) {
 560 //           TMailleSet::const_iterator ma_it = mailleSet.find( *maIt );
 561 //           if ( ma_it != mailleSet.end() ) {
 562 //             MESSAGE("EQUAL ELEMS: " << *ma_it << " AND " << *maIt);
 563 //           }
 564 //           else
 565             mailleSet.insert( *maIt );
 566         }
 567       }
 568       if ( nb_elem != mailleSet.size() ) { // Self intersecting compound group
 569         INFOS("Self intersecting group: " << i << " <" << grp.nom << ">"
 570               << ", mailleSet.size = " << mailleSet.size() << ", sum nb elems = " << nb_elem);
 571         for( vector<int>::iterator j=grp.groupes.begin(); j!=grp.groupes.end(); ++j)
 572           INFOS(" in sub-group "<<  *j << " <" << groupes[*j-1].nom << "> "
 573                 << groupes[*j-1].mailles.size() << " mailles of type "
 574                 << groupes[*j-1].mailles[0]->geometricType);
 575         // if a group serve as a support to a field, then the _field is to be converted
 576         // into several MEDMEM::FIELDs, one per sub-group; if we already skipped some sub-groups,
 577         // we roll back the loop on groups to create MED groups for skipped ones
 578         if ( isFieldSupport ) {
 579           if ( grp.nom.empty() ) // clear group existing for the sake of field only
 580             grp.groupes.clear();
 581           std::set<int> sub_grps;
 582           for ( fIt = fields.begin(); fIt != fields.end(); ++fIt ) {
 583             _fieldBase * field = (*fIt);
 584             if ( field->_group_id == i ) {
 585               field->_group_id = -1; // -> a field by support
 586               field->getGroupIds( sub_grps, false );
 587             }
 588           }
 589           if ( i > *sub_grps.begin() ) { // roll back
 590             support_groups.erase( i );
 591             support_groups.insert( sub_grps.begin(), sub_grps.end() );
 592             i = *sub_grps.begin() - 1;
 593             continue;
 594           }
 595         }
 596       }
 597     }
 598     else {
 599       _groupe::mailleIter maIt=grp.mailles.begin();
 600       for(; maIt!=grp.mailles.end(); ++maIt)
 601         mailleSet.insert( *maIt );
 602       if ( grp.mailles.size() != mailleSet.size() )
 603         INFOS( "Self intersecting group: " << i << " <" << grp.nom << ">"
 604                << ", mailleSet.size = " << mailleSet.size()
 605                << ", nb elems = " << grp.mailles.size());
 606     }
 607     // MEDMEM does not allow constituents of <mesh_dim>-2 and less dimension
 608     // but do not skip node group
 609     int group_min_dim = (**mailleSet.begin()).dimension();
 610     int group_max_dim = (**(--mailleSet.end())).dimension();
 611     if ( group_max_dim != 0 && group_min_dim <= dimension_maillage - 2  ) {
 612       INFOS("Skip group: " << i << " <" << grp.nom << "> - too small dimension: "
 613             << group_min_dim);
 614       continue;
 615     }
 616
 617     // 1. Build a map _maille* -> index in MEDMEM::GROUP.getNumber(MED_ALL_ELEMENTS).
 618     // It is used while fields building.
 619     // 2. make mailles know the groups they belong to, that is used in getFamilies()
 620     TMailleSet::iterator maIt = mailleSet.begin();
 621     int iMa;
 622     for ( iMa = 0; maIt != mailleSet.end(); maIt++ ) {
 623       grp.relocMap.insert( make_pair( &(**maIt), ++iMa ));
 624       (*maIt)->groupes.push_back( i );
 625     }
 626     ASSERT( iMa == grp.relocMap.size() );
 627
 628     int nb_geometric_types=1;
 629     TMailleSet::iterator j=mailleSet.begin();
 630     // initialise groupe_entity a l'entite de la premiere maille du groupe
 631     medEntityMesh groupe_entity = (**mailleSet.rbegin()).getEntity(dimension_maillage);
 632     if ( hasMixedCells )
 633       groupe_entity = MED_CELL;
 634     medGeometryElement geometrictype=(**j).geometricType;
 635
 636     //Parcours des mailles (a partir de la deuxieme) pour compter les types geometriques
 637     for ( ++j ; j!=mailleSet.end(); ++j )
 638     {
 639       //Compte nombre de types geometriques
 640       if ( (**j).geometricType != geometrictype ) // si on change de type geometrique
 641       {
 642         nb_geometric_types++;
 643         geometrictype=(**j).geometricType;
 644       }
 645     }
 646
 647     MED_EN::medGeometryElement * tab_types_geometriques = new MED_EN::medGeometryElement[nb_geometric_types];
 648     int * tab_index_types_geometriques = new int[nb_geometric_types+1];
 649     int * tab_numeros_elements = new int[mailleSet.size()];
 650     int * tab_nombres_elements = new int[nb_geometric_types];
 651
 652     //Remplit tableaux entree des methodes set
 653     int indice_mailles=0/*, maxOrdre = -1*/;
 654     j=mailleSet.begin();
 655     geometrictype=(**j).geometricType;
 656     tab_index_types_geometriques[0]=1;
 657     int indice_types_geometriques=1;
 658     tab_types_geometriques[0]=geometrictype;
 659     //parcours des mailles du groupe
 660     for (  ; j!=mailleSet.end(); ++j , ++indice_mailles)
 661     {
 662       const _maille& ma = **j;
 663       tab_numeros_elements[indice_mailles]= ma.ordre;
 664 //       if ( maxOrdre < tab_numeros_elements[indice_mailles] )
 665 //         maxOrdre = tab_numeros_elements[indice_mailles];
 666       if (ma.geometricType != geometrictype)
 667       {
 668         tab_index_types_geometriques[indice_types_geometriques]=indice_mailles+1;
 669         geometrictype=ma.geometricType;
 670         tab_types_geometriques[indice_types_geometriques]=geometrictype;
 671         ++indice_types_geometriques;
 672       }
 673     }
 674     tab_index_types_geometriques[indice_types_geometriques]=indice_mailles+1;
 675     for (int k=0; k != nb_geometric_types; ++k)
 676     {
 677       tab_nombres_elements[k] = tab_index_types_geometriques[k+1]-tab_index_types_geometriques[k];
 678     }
 679     //INFOS( "MAX ORDRE in grp " << grp.nom << " entity " << groupe_entity << " : " << maxOrdre);
 680
 681     //Determination type entite du groupe
 682     vector <GROUP *> * vect_group;
 683     switch ( groupe_entity )
 684     {
 685     case MED_CELL :
 686       vect_group= & _groupCell;
 687       break;
 688     case MED_FACE :
 689       vect_group= & _groupFace;
 690       break;
 691     case MED_EDGE :
 692       vect_group= & _groupEdge;
 693       break;
 694     case MED_NODE :
 695       vect_group= & _groupNode;
 696       break;
 697     }
 698     //Creation nouveau groupe MED
 699     GROUP * new_group = new GROUP();
 700     medGroupes[ i ] = new_group;
 701     //Appel methodes set
 702     new_group->setTotalNumberOfElements(mailleSet.size());
 703     new_group->setName(grp.nom);
 704     new_group->setMesh(_ptrMesh);
 705     new_group->setNumberOfGeometricType(nb_geometric_types);
 706     new_group->setGeometricType(tab_types_geometriques);
 707     new_group->setNumberOfElements(tab_nombres_elements);
 708     new_group->setNumber(tab_index_types_geometriques,tab_numeros_elements);
 709     new_group->setEntity(groupe_entity);
 710     new_group->setAll(mailleSet.size() == maillage.size());
 711
 712     vector<int> nbGaussPnt( nb_geometric_types, 1 );
 713     new_group->setNumberOfGaussPoint( &nbGaussPnt[0] );
 714
 715     vect_group->push_back(new_group);
 716
 717     delete [] tab_types_geometriques;
 718     delete [] tab_index_types_geometriques;
 719     delete [] tab_numeros_elements;
 720     delete [] tab_nombres_elements;
 721   }
 722
 723   END_OF(LOC);
 724 }
 725
 726 //=======================================================================
 727 //function : getFamilies
 728 //purpose  : create families like MESH::createFamilies() but preserves
 729 //           the order of elements in GROUPs defined by constituent families
 730 //           order. Call it after getGroups()
 731 //=======================================================================
 732
 733 void _intermediateMED::getFamilies(std::vector<FAMILY *> & _famCell,
 734                                    std::vector<FAMILY *> & _famFace,
 735                                    std::vector<FAMILY *> & _famEdge,
 736                                    std::vector<FAMILY *> & _famNode, MESH * _ptrMesh)
 737 {
 738   const char * LOC = "_intermediateMED::getFamilies() : ";
 739   BEGIN_OF(LOC);
 740
 741   int nbElemFam = 0, nbNodeFam = 0;
 742   std::map< GROUP*, vector< FAMILY * > > grpFamsMap;
 743   int dimension_maillage=maillage.rbegin()->dimension();
 744
 745   std::set<_maille>::const_iterator i=maillage.begin(); // iterateurs sur les mailles
 746   std::set<_maille>::const_iterator j=maillage.begin();
 747
 748   do
 749   {
 750     // make a family containing mailles shared by the same set of groups
 751     std::list<unsigned>&  grList = i->groupes;  // to define the family end
 752     int           dimension = i->dimension();        // to define the entity end
 753     medGeometryElement type = i->geometricType;
 754     medEntityMesh    entity = i->getEntity( dimension_maillage );
 755
 756     std::vector<medGeometryElement> tab_types_geometriques;
 757     std::vector<int> tab_index_types_geometriques;
 758     std::vector<int> tab_nombres_elements;
 759     std::vector<int> tab_numeros_elements;
 760
 761     int iMa = 1, nbtype = 0;
 762     tab_types_geometriques.push_back( type );
 763     tab_index_types_geometriques.push_back( iMa );
 764
 765     // scan family cells and fill the tab that are needed by the create a MED FAMILY
 766     while (i != maillage.end() &&
 767            i->groupes == grList &&
 768            i->dimension() == dimension)
 769     {
 770       if (type != i->geometricType) // si changement de type geometrique
 771       {
 772         tab_index_types_geometriques.push_back(iMa);
 773         tab_nombres_elements.push_back(nbtype);
 774         nbtype=0;
 775         type=i->geometricType;
 776         tab_types_geometriques.push_back(type); // stocke le nouveau type geometrique rencontre
 777       }
 778       ++nbtype;
 779       ++iMa;
 780       ++i;
 781     }
 782     tab_index_types_geometriques.push_back(iMa);
 783     tab_nombres_elements.push_back(nbtype); // n'a pas été stocké dans la boucle
 784
 785     tab_numeros_elements.resize( iMa - 1 );
 786     for ( iMa = 0; j != i; j++, iMa++ )
 787       tab_numeros_elements[ iMa ] = j->ordre;
 788
 789     int id = ( entity == MED_NODE ? ++nbNodeFam : -(++nbElemFam) );
 790
 791     ostringstream name;
 792     name << "FAM_" << id;
 793
 794     // create a empty MED FAMILY and fill it with the tabs we constructed
 795     FAMILY* newFam = new FAMILY();
 796     newFam->setTotalNumberOfElements( iMa );
 797     newFam->setName( name.str() );
 798     newFam->setMesh( _ptrMesh );
 799     newFam->setNumberOfGeometricType( tab_types_geometriques.size() );
 800     newFam->setGeometricType( &tab_types_geometriques[0] ); // we know the tab is not empy
 801     newFam->setNumberOfElements( &tab_nombres_elements[0] );
 802     newFam->setNumber( &tab_index_types_geometriques[0], &tab_numeros_elements[0] );
 803     newFam->setEntity( entity );
 804     newFam->setAll( false );
 805     newFam->setIdentifier( id );
 806     newFam->setNumberOfGroups( grList.size() );
 807
 808     // Update links between families and groups
 809     if ( ! grList.empty() )
 810     {
 811       std::string * groupNames = new string[ grList.size() ];
 812       std::list<unsigned>::iterator g = grList.begin();
 813       for ( int i = 0; g != grList.end(); ++g, ++i ) {
 814         GROUP * medGROUP = getGroup( *g );
 815         groupNames[ i ] = medGROUP->getName();
 816         grpFamsMap[ medGROUP ].push_back( newFam );
 817       }
 818       newFam->setGroupsNames(groupNames);
 819     }
 820     // store newFam
 821     std::vector<FAMILY*>* families = 0;
 822     switch ( entity )
 823     {
 824     case MED_CELL :
 825       families = & _famCell; break;
 826     case MED_FACE :
 827       families = & _famFace; break;
 828     case MED_EDGE :
 829       families = & _famEdge; break;
 830     case MED_NODE :
 831       families = & _famNode; break;
 832     }
 833     if ( families )
 834       families->push_back( newFam );
 835
 836   } while ( i != maillage.end() );
 837
 838   // update references in groups
 839   std::map< GROUP*, vector< FAMILY * > >::iterator gf = grpFamsMap.begin();
 840   for ( ; gf != grpFamsMap.end(); ++gf )
 841   {
 842     gf->first->setNumberOfFamilies( gf->second.size() );
 843     gf->first->setFamilies( gf->second );
 844   }
 845 }
 846
 847 //=======================================================================
 848 //function : getGroup
 849 //purpose  :
 850 //=======================================================================
 851
 852 GROUP * _intermediateMED::getGroup( int i )
 853 {
 854   if ( i < medGroupes.size() )
 855     return medGroupes[ i ];
 856   throw MEDEXCEPTION
 857     (LOCALIZED(STRING("_intermediateMED::getGroup(): WRONG GROUP INDEX: ")
 858                << medGroupes.size() << " <= " << i ));
 859 }
 860
 861 //=======================================================================
 862 //function : getFields
 863 //purpose  :
 864 //=======================================================================
 865
 866 void _intermediateMED::getFields(std::list< FIELD_* >& theFields)
 867 {
 868   const char * LOC = "_intermediateMED::getFields() : ";
 869   BEGIN_OF(LOC);
 870   std::list< _fieldBase* >::const_iterator fIt = fields.begin();
 871   for ( ; fIt != fields.end(); fIt++ )
 872   {
 873     const _fieldBase* fb = *fIt;
 874     list<pair< FIELD_*, int> >  ff = fb->getField(groupes);
 875     list<pair< FIELD_*, int> >::iterator f_sup = ff.begin();
 876     for (int j = 1 ; f_sup != ff.end(); f_sup++, ++j )
 877     {
 878       FIELD_* f = f_sup->first;
 879       SUPPORT* sup = getGroup( f_sup->second );
 880       if ( !sup )
 881         throw MEDEXCEPTION
 882           (LOCALIZED(STRING(LOC) <<"_intermediateMED::getFields(), NULL field support: "
 883                      << " group index: " << fb->_group_id));
 884       int nb_elems = sup->getNumberOfElements( MED_ALL_ELEMENTS );
 885       if ( nb_elems != f->getNumberOfValues() )
 886         throw MEDEXCEPTION
 887           (LOCALIZED(STRING("_intermediateMED::getFields(), field support size (")
 888                      << nb_elems  << ") != NumberOfValues (" << f->getNumberOfValues()));
 889       theFields.push_back( f );
 890       if ( sup->getName().empty() )
 891         sup->setName( "GRP_" + f->getName() );
 892       f->setSupport( sup );
 893       //f->setIterationNumber( j );
 894       f->setOrderNumber( j );
 895     }
 896   }
 897   END_OF(LOC);
 898 }
 899
 900 _intermediateMED::~_intermediateMED()
 901 {
 902   MESSAGE( "~_intermediateMED()");
 903   std::list< _fieldBase* >::const_iterator fIt = fields.begin();
 904   for ( ; fIt != fields.end(); fIt++ )
 905     delete *fIt;
 906 }
 907
 908 //=======================================================================
 909 //function : getGroupIds
 910 //purpose  : return ids of main and/or sub- groups
 911 //=======================================================================
 912
 913 void _fieldBase::getGroupIds( std::set<int> & ids, bool all ) const
 914 {
 915   if ( hasCommonSupport() )
 916     ids.insert( _group_id );
 917   if ( all || !hasCommonSupport() ) {
 918     vector< _sub_data >::const_iterator sub = _sub.begin();
 919     for ( ; sub != _sub.end(); ++sub )
 920       ids.insert( sub->_supp_id );
 921   }
 922 }
 923
 924 //=======================================================================
 925 //function : hasSameComponentsBySupport
 926 //purpose  :
 927 //=======================================================================
 928
 929 bool _fieldBase::hasSameComponentsBySupport() const
 930 {
 931   vector< _sub_data >::const_iterator sub_data = _sub.begin();
 932   const _sub_data& first_sub_data = *sub_data;
 933   for ( ++sub_data ; sub_data != _sub.end(); ++sub_data )
 934   {
 935     if ( first_sub_data.nbComponents() != sub_data->nbComponents() )
 936       return false;
 937     for ( int iComp = 0; iComp < first_sub_data.nbComponents(); ++iComp )
 938       if ( first_sub_data._comp_names[ iComp ] != sub_data->_comp_names[ iComp ])
 939         return false;
 940   }
 941   return true;
 942 }
 943
 944 /////