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 100
  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: " << nb_subcomponents << endl <<
 130       "  nb comp: " << nb_components << endl <<
 131         "  group index: " << group_id << endl <<
 132           "  type: " << type << endl;
 133   os << "  comp names: ";
 134   for ( int i = 0; i < comp_names.size(); ++i )
 135     os << " |" << comp_names[ i ] << "|";
 136 }
 137
 138 std::ostream& MEDMEM::operator << (std::ostream& os, const _fieldBase * f)
 139 {
 140   f->dump( os );
 141   return os;
 142 }
 143
 144 std::ostream& MEDMEM::operator << (std::ostream& os, const _intermediateMED& mi)
 145 {
 146     os << "Set des " << mi.maillage.size() << " mailles : " << std::endl;
 147     std::set<_maille>::const_iterator i=mi.maillage.begin();
 148     int l;
 149     for( l = 0; l < DUMP_LINES_LIMIT && i!=mi.maillage.end(); ++i, ++l)
 150         os << setw(3) << l+1 << " " <<*i << std::endl;
 151     if ( l == DUMP_LINES_LIMIT )
 152       os << "   ... skip " << mi.maillage.size() - l << " mailles" << endl;
 153
 154     os << std::endl << "Vector des " << mi.groupes.size() << " groupes : " << std::endl;
 155     for (unsigned int k=0; k!=mi.groupes.size(); k++)
 156       os << std::setw(3) << k << " " << mi.groupes[k] << std::endl;
 157
 158     os << std::endl << "map des " << mi.points.size() << " noeuds : " << std::endl;
 159     std::map<int,_noeud>::const_iterator j=mi.points.begin();
 160     for( l = 0; l < DUMP_LINES_LIMIT && j!=mi.points.end(); ++j, ++l)
 161         os << j->second << std::endl;
 162     if ( l == DUMP_LINES_LIMIT )
 163       os << "   ... skip " << mi.points.size() - l << " noeuds" << endl;
 164
 165     os << endl << mi.fields.size() << " fields:" << endl;
 166     std::list<_fieldBase* >::const_iterator fIt = mi.fields.begin();
 167     for ( l = 0; fIt != mi.fields.end(); ++fIt, ++l )
 168       os << " - " << l+1 << " " << *fIt << endl;
 169
 170     return os;
 171 }
 172
 173
 174 //=======================================================================
 175 //function : treatGroupes
 176 //purpose  : detect groupes of mixed dimension and erase groupes that
 177 //           won't be converted
 178 //=======================================================================
 179
 180 void _intermediateMED::treatGroupes()
 181 {
 182   const char * LOC = "_intermediateMED::treatGroupes() : ";
 183   BEGIN_OF(LOC);
 184
 185   // --------------------
 186   // erase useless group
 187   // --------------------
 188
 189   // decrease hierarchical depth of subgroups
 190   vector<int>::iterator j;
 191   for (unsigned int i=0; i!=this->groupes.size(); ++i)
 192   {
 193     _groupe& grp = groupes[i];
 194     //INFOS( i << " " << grp.nom );
 195     j = grp.groupes.begin();
 196     while( j!=grp.groupes.end() ) {
 197       int grpInd = *j-1;
 198       if ( grpInd < 0 || grpInd >= groupes.size() ) {
 199         throw MEDEXCEPTION(LOCALIZED(STRING(LOC) << "Bad subgroup index: " << grpInd <<
 200                                      ", in " << i << " groupe.nom=" << grp.nom));
 201       }
 202       _groupe & sub_grp = groupes[ grpInd ];
 203       if ( !sub_grp.groupes.empty() ) {
 204         MESSAGE("High hierarchical depth of subgroups in group " << i );
 205         *j = sub_grp.groupes[0]; // replace j with its 1st subgroup
 206         // push back the rest subs
 207         for ( int k = 1; k < sub_grp.groupes.size(); ++k )
 208           grp.groupes.push_back( sub_grp.groupes[ k ]);
 209         // vector maybe is reallocated: restart iterator
 210         j = grp.groupes.begin();
 211       }
 212       else
 213         j++;
 214     }
 215     // remove empty sub-groupes
 216     j = grp.groupes.begin();
 217     while ( j!=grp.groupes.end() ) {
 218       if ( groupes[*j-1].empty() ) {
 219         grp.groupes.erase( j );
 220         j = grp.groupes.begin();
 221       }
 222       else
 223         j++;
 224     }
 225   }
 226   // get indices of groups that are field support -
 227   // do not erase them and their subgroups
 228   std::set<int> groups2convert;
 229   std::list< _fieldBase* >::const_iterator fIt = fields.begin();
 230   for ( ; fIt != fields.end(); ++fIt )
 231   {
 232     groups2convert.insert( (*fIt)->group_id );
 233     _groupe& grp = groupes[ (*fIt)->group_id ];
 234     for( j = grp.groupes.begin(); j!=grp.groupes.end(); ++j)
 235       groups2convert.insert( *j-1 );
 236   }
 237   // keep named groups and their subgroups
 238   for (unsigned int i=0; i!=this->groupes.size(); ++i)
 239   {
 240     _groupe& grp = groupes[i];
 241     if ( grp.empty() || grp.nom.empty() )
 242       continue;
 243     groups2convert.insert( i );
 244     for( j = grp.groupes.begin(); j!=grp.groupes.end(); ++j)
 245       groups2convert.insert( *j-1 );
 246   }
 247   // erase groups that are not in groups2convert
 248   for (unsigned int i=0; i!=this->groupes.size(); ++i)
 249   {
 250     if ( groups2convert.find( i ) == groups2convert.end() ) {
 251       _groupe& grp = groupes[i];
 252       grp.mailles.clear();
 253       grp.groupes.clear();
 254       //INFOS( "Erase " << i << "-th group " << grp.nom );
 255     }
 256   }
 257
 258   // ---------------------------------------------------
 259   // define if there are groups with mixed entity types
 260   // ---------------------------------------------------
 261
 262   hasMixedCells = false;
 263   for (unsigned int i=0; i!=this->groupes.size(); ++i)
 264   {
 265     _groupe& grp = groupes[i];
 266     if ( grp.groupes.empty() )
 267       continue;
 268
 269     // check if sub-groups have different dimension
 270     j = grp.groupes.begin();
 271     int dim = groupes[*j-1].mailles[0]->dimension();
 272     for( j++; !hasMixedCells && j!=grp.groupes.end(); ++j)
 273       hasMixedCells = ( dim != groupes[*j-1].mailles[0]->dimension() );
 274   }
 275
 276   if ( hasMixedCells )
 277     INFOS( "There will be groups of mixed dimention" );
 278   END_OF(LOC);
 279 }
 280
 281 void _intermediateMED::numerotationMaillage()
 282 {
 283   const char * LOC = "_intermediateMED::numerotationMaillage() : ";
 284   BEGIN_OF(LOC);
 285
 286   treatGroupes();
 287
 288   // numerotation mailles of type MED_POINT1 by node number
 289   std::set<_maille>::iterator i=maillage.begin();
 290   if ( i->geometricType == MED_POINT1 ) {
 291     numerotationPoints();
 292     while ( i!=maillage.end() && i->geometricType == MED_POINT1 ) {
 293       i->ordre = i->sommets[0]->second.number;
 294       i++;
 295     }
 296   }
 297   // numerotation des mailles par entité
 298     int i_maille=0;
 299     int dimension=i->dimension();
 300     for( ; i!=maillage.end(); ++i)
 301     {
 302         if ( !hasMixedCells && dimension!=i->dimension() ) // on change d'entite
 303         {
 304           MESSAGE( "NB dim " << dimension << " entities: " << i_maille);
 305             dimension=i->dimension();
 306             i_maille=0;
 307         }
 308         (*i).ordre=++i_maille;
 309     }
 310   END_OF(LOC);
 311 }
 312
 313 void _intermediateMED::numerotationPoints()
 314 {
 315     // Fonction de renumerotation des noeuds (necessaire quand il y a des trous dans la numerotation.
 316     int i_noeud=0;
 317     for( std::map<int,_noeud>::const_iterator i=points.begin(); i!=points.end(); ++i)
 318         i->second.number = ++i_noeud ;
 319 }
 320
 321
 322 /*!
 323  * \if developper
 324  * create a MED COORDINATE from the intermediate structure.
 325  * \endif
 326  */
 327 COORDINATE *
 328 _intermediateMED::getCoordinate(const string & coordinateSystem)
 329 {
 330     const medModeSwitch mode=MED_FULL_INTERLACE;
 331     int spaceDimension=points.begin()->second.coord.size();
 332     int numberOfNodes=points.size();
 333
 334
 335     // creation du tableau des coordonnees en mode MED_FULL_INTERLACE
 336     double * coord = new double[spaceDimension*numberOfNodes];
 337     int k=0;
 338     for( std::map<int,_noeud>::const_iterator i=points.begin(); i!=points.end(); ++i, ++k)
 339         std::copy(i->second.coord.begin(), i->second.coord.end(), coord+k*spaceDimension);
 340
 341     // creation de l'objet COORDINATE
 342     COORDINATE * coordinate = new COORDINATE(spaceDimension, numberOfNodes, mode);
 343     coordinate->setCoordinates(mode,coord);
 344     delete [] coord;
 345     coordinate->setCoordinatesSystem(coordinateSystem);
 346     return coordinate;
 347 }
 348
 349
 350   /*!
 351    * \if developper
 352    * create a MED CONNECTIVITY from the intermediate structure.
 353    * for memory optimisation, clear the intermediate structure (the "maillage" set )
 354    * \endif
 355    */
 356 CONNECTIVITY *
 357 _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
 369     medGeometryElement type=0; // variables de travail
 370     int nbtype=0;
 371     int dimension=0;
 372     bool first = true;
 373
 374     std::vector<medGeometryElement> vtype; // tableau de travail : stockage des types geometriques pour UNE entite
 375     std::vector<int> vcount; // tableau de travail : nombre d'elements pour chaque type geometrique de vtype
 376
 377     std::set<_maille>::const_iterator i=maillage.begin(); // iterateurs sur les mailles
 378     std::set<_maille>::const_iterator j=maillage.begin();
 379
 380     // renumerote les points de 1 a n (pour le cas ou certains points ne sont pas presents dans le maillage d'origine)
 381     numerotationPoints();
 382
 383     do
 384     {
 385         // boucle sur les entites de la structure maillage :
 386         //   - on parcourt une premiere fois avec i les mailles d'une entite pour récupérer
 387         //     des infos sur les types geometriques, leur nombre et le nombre d'elements.
 388         //   - on alloue la connectivite
 389         //   - on parcourt une deuxieme fois avec j pour lire les noeuds.
 390
 391
 392         type=i->geometricType; // init boucle for
 393         dimension=i->dimension();
 394         nbtype=0;
 395         vtype.push_back(type);
 396         // Boucle sur i de parcours des mailles d'une entite
 397         // Une entite se termine lorsqu'on atteint la fin de maillage ou lorsque la dimension des mailles change
 398         for( ; i!=maillage.end() && ( hasMixedCells || dimension==i->dimension()) ; ++i)
 399         {
 400             if (type != i->geometricType) // si changement de type geometrique
 401             {
 402                 vcount.push_back(nbtype); // stocke le nombre d'elements de type nbtype
 403                 nbtype=0;
 404                 type=i->geometricType;
 405                 vtype.push_back(type); // stocke le nouveau type geometrique rencontre
 406             }
 407
 408             ++nbtype;
 409         }
 410         vcount.push_back(dimension ? nbtype : numberOfNodes); // n'a pas été stocké dans la boucle
 411
 412         // Pour l'entite qu'on vient de parcourir, allocation des tableau et creation de la connectivite
 413 //      cout << "Dimension = " << dimension << endl;
 414 //      cout << "Nombre de type geometriques : " << vtype.size() << endl;
 415 //      for (unsigned k=0; k!=vtype.size(); ++k )
 416 //          cout << "  -> " << vtype[k] << " : " << vcount[k] << endl;
 417
 418         numberOfTypes=vtype.size(); // nombre de types de l'entite
 419
 420         if ( i==maillage.end() ) // cas de l'entite de plus haut niveau
 421             entity=MED_CELL;
 422         else if (dimension==2 )
 423             entity=MED_FACE;
 424         else if (dimension==1 )
 425             entity=MED_EDGE;
 426         else if (dimension==0 )
 427             entity=MED_NODE;
 428
 429         Connectivity = new CONNECTIVITY(numberOfTypes,entity);
 430         Connectivity->setEntityDimension(dimension);
 431         Connectivity->setNumberOfNodes(numberOfNodes);
 432
 433         types = new medGeometryElement[numberOfTypes];
 434         std::copy(vtype.begin(),vtype.end(),types);
 435         Connectivity->setGeometricTypes(types,entity);
 436         delete [] types;
 437
 438         count = new int[numberOfTypes+1];
 439         count[0]=1;
 440         for (unsigned int k=0; k!=vcount.size(); ++k )
 441           count[k+1]=count[k]+vcount[k];
 442         Connectivity->setCount (count, entity);
 443         SCRUTE( entity );
 444         SCRUTE(numberOfTypes);
 445         SCRUTE(count[numberOfTypes]-1);
 446         delete [] count;
 447
 448         for (int k=0; k!=numberOfTypes; ++k )
 449           {
 450             // pour chaque type géometrique k, copie des sommets dans connectivity et set dans Connectivity
 451             int nbSommetsParMaille = j->sommets.size();
 452             int n, nbSommets = vcount[k] * j->sommets.size();
 453             connectivity = new int[ nbSommets ];
 454             for (int l=0; l!=vcount[k]; ++l)
 455             {
 456                 if ( entity==MED_NODE )
 457                   connectivity[l] = l+1;
 458                 else
 459                 {
 460                   for ( n=0; n != nbSommetsParMaille; ++n) {
 461                     connectivity[nbSommetsParMaille*l+n] =
 462                       j->sommets[ j->reverse ? nbSommetsParMaille-n-1 : n ]->second.number;
 463                   }
 464                 // DO NOT ERASE, maillage will be used while fields construction
 465                 //maillage.erase(j);    ; // dangereux, mais optimise la mémoire consommée!
 466                   j++;
 467                 }
 468             }
 469
 470             Connectivity->setNodal  (connectivity, entity, vtype[k]);
 471             delete [] connectivity;
 472           }
 473
 474         if ( entity==MED_NODE )
 475           j = i;
 476         else if (i!=j)
 477             throw MEDEXCEPTION(LOCALIZED(STRING(LOC) << "Erreur de lecture des mailles ."));
 478
 479         if ( ! first)
 480             Connectivity->setConstituent (Constituent);
 481         else
 482             first = false;
 483         Constituent = Connectivity; // stocke Connectivity pour utilisation dans setConstituent lors de la boucle suivante
 484
 485         vtype.clear();
 486         vcount.clear();
 487
 488     }
 489     while ( i != maillage.end() );
 490
 491     END_OF(LOC);
 492     return Connectivity;
 493 }
 494
 495
 496   /*!
 497    * \if developper
 498    * fill the arguments vector of groups from the intermediate structure.
 499    * This function must be called before getConnectivity()
 500    * \endif
 501    */
 502 void
 503 _intermediateMED::getGroups(std::vector<GROUP *> & _groupCell, std::vector<GROUP *> & _groupFace, std::vector<GROUP *> & _groupEdge, std::vector<GROUP *> & _groupNode, MESH * _ptrMesh)
 504 {
 505   const char * LOC = "_intermediateMED::getGroups() : ";
 506   BEGIN_OF(LOC);
 507   if (maillage.size() == 0) {
 508     INFOS( "Erreur : il n'y a plus de mailles");
 509     return;
 510   }
 511
 512   // get indices of groups that are field support - do not skip them
 513   std::set<int> support_groups;
 514   std::list< _fieldBase* >::const_iterator fIt = fields.begin();
 515   for ( ; fIt != fields.end(); ++fIt )
 516     support_groups.insert( (*fIt)->group_id );
 517
 518   numerotationMaillage(); // Renumerotation des mailles par entite
 519
 520   int dimension_maillage=maillage.rbegin()->dimension();
 521
 522   for (unsigned int i=0; i!=this->groupes.size(); ++i)
 523   {
 524     _groupe& grp = groupes[i];
 525     // si le groupe est vide, ou si le groupe n'est pas nommé : on passe au suivant
 526     if ( grp.empty() ||
 527         ( grp.nom.empty() && support_groups.find( i ) == support_groups.end() )) {
 528       if ( !grp.nom.empty() )
 529         INFOS("Skip group " << grp.nom );
 530       medGroupes.push_back( NULL );
 531       continue;
 532     }
 533
 534     // Build a set of mailles: sort mailles by type and exclude maille doubling
 535     typedef set< set<_maille>::iterator, _mailleIteratorCompare > TMailleSet;
 536     TMailleSet mailleSet;
 537     if( grp.groupes.size() ) {// le groupe i contient des sous-maillages
 538       int nb_elem = 0;
 539       for( vector<int>::iterator j=grp.groupes.begin(); j!=grp.groupes.end(); ++j)
 540       {
 541         nb_elem += groupes[*j-1].mailles.size();
 542         _groupe::mailleIter maIt=groupes[*j-1].mailles.begin();
 543         for( ; maIt!=groupes[*j-1].mailles.end(); ++maIt) {
 544 //           TMailleSet::const_iterator ma_it = mailleSet.find( *maIt );
 545 //           if ( ma_it != mailleSet.end() ) {
 546 //             MESSAGE("EQUAL ELEMS: " << *ma_it << " AND " << *maIt);
 547 //           }
 548 //           else
 549             mailleSet.insert( *maIt );
 550         }
 551       }
 552       if ( nb_elem != mailleSet.size() ) {
 553         INFOS("Self intersecting group: " << i+1 << " <" << grp.nom << ">"
 554               << ", mailleSet.size = " << mailleSet.size() << ", sum nb elems = " << nb_elem);
 555         for( vector<int>::iterator j=grp.groupes.begin(); j!=grp.groupes.end(); ++j)
 556           INFOS(" in sub-group "<<  *j << " <" << groupes[*j-1].nom << "> "
 557                 << groupes[*j-1].mailles.size() << " mailles of type "
 558                 << groupes[*j-1].mailles[0]->geometricType);
 559       }
 560     }
 561     else {
 562       _groupe::mailleIter maIt=grp.mailles.begin();
 563       for(; maIt!=grp.mailles.end(); ++maIt)
 564         mailleSet.insert( *maIt );
 565       if ( grp.mailles.size() != mailleSet.size() )
 566         INFOS( "Self intersecting group: " << i+1 << " <" << grp.nom << ">"
 567               << ", mailleSet.size = " << mailleSet.size() << ", nb elems = " << grp.mailles.size());
 568     }
 569
 570     // 1. Build a map _maille* -> index in MEDMEM::GROUP.getNumber(MED_ALL_ELEMENTS).
 571     // It is used while fields building.
 572     // 2. make mailles know the groups they belong to, that is used in getFamilies()
 573     TMailleSet::iterator maIt = mailleSet.begin();
 574     int iMa;
 575     for ( iMa = 0; maIt != mailleSet.end(); maIt++ ) {
 576       grp.relocMap.insert( make_pair( &(**maIt), ++iMa ));
 577       (*maIt)->groupes.push_back( i );
 578     }
 579     ASSERT( iMa == grp.relocMap.size() );
 580
 581     int nb_geometric_types=1;
 582     TMailleSet::iterator j=mailleSet.begin();
 583     // initialise groupe_entity a l'entite de la premiere maille du groupe
 584     medEntityMesh groupe_entity = (**mailleSet.rbegin()).getEntity(dimension_maillage);
 585     if ( hasMixedCells )
 586       groupe_entity = MED_CELL;
 587     medGeometryElement geometrictype=(**j).geometricType;
 588
 589     //Parcours des mailles (a partir de la deuxieme) pour compter les types geometriques
 590     for ( ++j ; j!=mailleSet.end(); ++j )
 591     {
 592       //Compte nombre de types geometriques
 593       if ( (**j).geometricType != geometrictype ) // si on change de type geometrique
 594       {
 595         nb_geometric_types++;
 596         geometrictype=(**j).geometricType;
 597       }
 598     }
 599
 600     MED_EN::medGeometryElement * tab_types_geometriques = new MED_EN::medGeometryElement[nb_geometric_types];
 601     int * tab_index_types_geometriques = new int[nb_geometric_types+1];
 602     int * tab_numeros_elements = new int[mailleSet.size()];
 603     int * tab_nombres_elements = new int[nb_geometric_types];
 604
 605     //Remplit tableaux entree des methodes set
 606     int indice_mailles=0/*, maxOrdre = -1*/;
 607     j=mailleSet.begin();
 608     geometrictype=(**j).geometricType;
 609     tab_index_types_geometriques[0]=1;
 610     int indice_types_geometriques=1;
 611     tab_types_geometriques[0]=geometrictype;
 612     //parcours des mailles du groupe
 613     for (  ; j!=mailleSet.end(); ++j , ++indice_mailles)
 614     {
 615       const _maille& ma = **j;
 616       tab_numeros_elements[indice_mailles]= ma.ordre;
 617 //       if ( maxOrdre < tab_numeros_elements[indice_mailles] )
 618 //         maxOrdre = tab_numeros_elements[indice_mailles];
 619       if (ma.geometricType != geometrictype)
 620       {
 621         tab_index_types_geometriques[indice_types_geometriques]=indice_mailles+1;
 622         geometrictype=ma.geometricType;
 623         tab_types_geometriques[indice_types_geometriques]=geometrictype;
 624         ++indice_types_geometriques;
 625       }
 626     }
 627     tab_index_types_geometriques[indice_types_geometriques]=indice_mailles+1;
 628     for (int k=0; k != nb_geometric_types; ++k)
 629     {
 630       tab_nombres_elements[k] = tab_index_types_geometriques[k+1]-tab_index_types_geometriques[k];
 631     }
 632     //INFOS( "MAX ORDRE in grp " << grp.nom << " entity " << groupe_entity << " : " << maxOrdre);
 633
 634     //Determination type entite du groupe
 635     vector <GROUP *> * vect_group;
 636     switch ( groupe_entity )
 637     {
 638     case MED_CELL :
 639       vect_group= & _groupCell;
 640       break;
 641     case MED_FACE :
 642       vect_group= & _groupFace;
 643       break;
 644     case MED_EDGE :
 645       vect_group= & _groupEdge;
 646       break;
 647     case MED_NODE :
 648       vect_group= & _groupNode;
 649       break;
 650     }
 651     //Creation nouveau groupe MED
 652     GROUP * new_group = new GROUP();
 653     //Appel methodes set
 654     new_group->setTotalNumberOfElements(mailleSet.size());
 655     new_group->setName(grp.nom);
 656     new_group->setMesh(_ptrMesh);
 657     new_group->setNumberOfGeometricType(nb_geometric_types);
 658     new_group->setGeometricType(tab_types_geometriques);
 659     new_group->setNumberOfElements(tab_nombres_elements);
 660     new_group->setNumber(tab_index_types_geometriques,tab_numeros_elements);
 661     new_group->setEntity(groupe_entity);
 662     new_group->setAll(mailleSet.size() == maillage.size());
 663
 664     vector<int> nbGaussPnt( nb_geometric_types, 1 );
 665     new_group->setNumberOfGaussPoint( &nbGaussPnt[0] );
 666
 667     vect_group->push_back(new_group);
 668     medGroupes.push_back( new_group );
 669
 670     delete [] tab_types_geometriques;
 671     delete [] tab_index_types_geometriques;
 672     delete [] tab_numeros_elements;
 673     delete [] tab_nombres_elements;
 674   }
 675   SCRUTE( medGroupes.size() );
 676
 677   END_OF(LOC);
 678 }
 679
 680 //=======================================================================
 681 //function : getFamilies
 682 //purpose  : create families like MESH::createFamilies() but preserves
 683 //           the order of elements in GROUPs defined by constituent families
 684 //           order. Call it after getGroups()
 685 //=======================================================================
 686
 687 void _intermediateMED::getFamilies(std::vector<FAMILY *> & _famCell,
 688                                    std::vector<FAMILY *> & _famFace,
 689                                    std::vector<FAMILY *> & _famEdge,
 690                                    std::vector<FAMILY *> & _famNode, MESH * _ptrMesh)
 691 {
 692   const char * LOC = "_intermediateMED::getFamilies() : ";
 693   BEGIN_OF(LOC);
 694
 695   int nbElemFam = 0, nbNodeFam = 0;
 696   std::map< GROUP*, vector< FAMILY * > > grpFamsMap;
 697   int dimension_maillage=maillage.rbegin()->dimension();
 698
 699   std::set<_maille>::const_iterator i=maillage.begin(); // iterateurs sur les mailles
 700   std::set<_maille>::const_iterator j=maillage.begin();
 701
 702   do
 703   {
 704     // make a family containing mailles shared by the same set of groups
 705     std::list<unsigned>&  grList = i->groupes;  // to define the family end
 706     int           dimension = i->dimension();        // to define the entity end
 707     medGeometryElement type = i->geometricType;
 708     medEntityMesh    entity = i->getEntity( dimension_maillage );
 709
 710     std::vector<medGeometryElement> tab_types_geometriques;
 711     std::vector<int> tab_index_types_geometriques;
 712     std::vector<int> tab_nombres_elements;
 713     std::vector<int> tab_numeros_elements;
 714
 715     int iMa = 1, nbtype = 0;
 716     tab_types_geometriques.push_back( type );
 717     tab_index_types_geometriques.push_back( iMa );
 718
 719     // scan family cells and fill the tab that are needed by the create a MED FAMILY
 720     while (i != maillage.end() &&
 721            i->groupes == grList &&
 722            i->dimension() == dimension)
 723     {
 724       if (type != i->geometricType) // si changement de type geometrique
 725       {
 726         tab_index_types_geometriques.push_back(iMa);
 727         tab_nombres_elements.push_back(nbtype);
 728         nbtype=0;
 729         type=i->geometricType;
 730         tab_types_geometriques.push_back(type); // stocke le nouveau type geometrique rencontre
 731       }
 732       ++nbtype;
 733       ++iMa;
 734       ++i;
 735     }
 736     tab_index_types_geometriques.push_back(iMa);
 737     tab_nombres_elements.push_back(nbtype); // n'a pas été stocké dans la boucle
 738
 739     tab_numeros_elements.resize( iMa - 1 );
 740     for ( iMa = 0; j != i; j++, iMa++ )
 741       tab_numeros_elements[ iMa ] = j->ordre;
 742
 743     int id = ( entity == MED_NODE ? ++nbNodeFam : -(++nbElemFam) );
 744
 745     ostringstream name;
 746     name << "FAM_" << id;
 747
 748     // create a empty MED FAMILY and fill it with the tabs we constructed
 749     FAMILY* newFam = new FAMILY();
 750     newFam->setTotalNumberOfElements( iMa );
 751     newFam->setName( name.str() );
 752     newFam->setMesh( _ptrMesh );
 753     newFam->setNumberOfGeometricType( tab_types_geometriques.size() );
 754     newFam->setGeometricType( &tab_types_geometriques[0] ); // we know the tab is not empy
 755     newFam->setNumberOfElements( &tab_nombres_elements[0] );
 756     newFam->setNumber( &tab_index_types_geometriques[0], &tab_numeros_elements[0] );
 757     newFam->setEntity( entity );
 758     newFam->setAll( false );
 759     newFam->setIdentifier( id );
 760     newFam->setNumberOfGroups( grList.size() );
 761
 762     // Update links between families and groups
 763     if ( ! grList.empty() )
 764     {
 765       std::string * groupNames = new string[ grList.size() ];
 766       std::list<unsigned>::iterator g = grList.begin();
 767       for ( int i = 0; g != grList.end(); ++g, ++i ) {
 768         GROUP * medGROUP = getGroup( *g );
 769         groupNames[ i ] = medGROUP->getName();
 770         grpFamsMap[ medGROUP ].push_back( newFam );
 771       }
 772       newFam->setGroupsNames(groupNames);
 773     }
 774     // store newFam
 775     std::vector<FAMILY*>* families = 0;
 776     switch ( entity )
 777     {
 778     case MED_CELL :
 779       families = & _famCell; break;
 780     case MED_FACE :
 781       families = & _famFace; break;
 782     case MED_EDGE :
 783       families = & _famEdge; break;
 784     case MED_NODE :
 785       families = & _famNode; break;
 786     }
 787     if ( families )
 788       families->push_back( newFam );
 789
 790   } while ( i != maillage.end() );
 791
 792   // update references in groups
 793   std::map< GROUP*, vector< FAMILY * > >::iterator gf = grpFamsMap.begin();
 794   for ( ; gf != grpFamsMap.end(); ++gf )
 795   {
 796     gf->first->setNumberOfFamilies( gf->second.size() );
 797     gf->first->setFamilies( gf->second );
 798   }
 799 }
 800
 801 //=======================================================================
 802 //function : getGroup
 803 //purpose  :
 804 //=======================================================================
 805
 806 GROUP * _intermediateMED::getGroup( int i )
 807 {
 808   if ( i < medGroupes.size() )
 809     return medGroupes[ i ];
 810   throw MEDEXCEPTION
 811     (LOCALIZED(STRING("_intermediateMED::getGroup(): WRONG GROUP INDEX: ")
 812                << medGroupes.size() << " <= " << i ));
 813 }
 814
 815 //=======================================================================
 816 //function : getFields
 817 //purpose  :
 818 //=======================================================================
 819
 820 void _intermediateMED::getFields(std::list< FIELD_* >& theFields)
 821 {
 822   const char * LOC = "_intermediateMED::getFields() : ";
 823   BEGIN_OF(LOC);
 824   std::list< _fieldBase* >::const_iterator fIt = fields.begin();
 825   for ( ; fIt != fields.end(); fIt++ )
 826   {
 827     const _fieldBase* fb = *fIt;
 828     SUPPORT* sup = getGroup( fb->group_id );
 829     if ( !sup )
 830       throw MEDEXCEPTION
 831         (LOCALIZED(STRING(LOC) <<"_intermediateMED::getFields(), NULL field support: "
 832                    << " group index: " << fb->group_id));
 833     int nb_elems = sup->getNumberOfElements( MED_ALL_ELEMENTS );
 834
 835     std::list< FIELD_* > ff = fb->getField(groupes);
 836     std::list< FIELD_* >::iterator it = ff.begin();
 837     for (int j = 1 ; it != ff.end(); it++, ++j )
 838     {
 839       FIELD_* f = *it;
 840       if ( nb_elems != f->getNumberOfValues() )
 841         throw MEDEXCEPTION
 842           (LOCALIZED(STRING("_intermediateMED::getFields(), field support size (")
 843                      << nb_elems  << ") != NumberOfValues (" << f->getNumberOfValues()));
 844       theFields.push_back( f );
 845       f->setSupport( sup );
 846       //f->setIterationNumber( j );
 847       f->setOrderNumber( j );
 848     }
 849   }
 850   END_OF(LOC);
 851 }
 852
 853 _intermediateMED::~_intermediateMED()
 854 {
 855   MESSAGE( "~_intermediateMED()");
 856   std::list< _fieldBase* >::const_iterator fIt = fields.begin();
 857   for ( ; fIt != fields.end(); fIt++ )
 858     delete *fIt;
 859 }
 860
 861
 862 /////