src/MEDMEM/MEDMEM_GibiMeshDriver.cxx

   1 #include <algorithm>
   2 #include <queue>
   3
   4 #include "MEDMEM_GibiMeshDriver.hxx"
   5
   6 #include "MEDMEM_DriversDef.hxx"
   7
   8 #include "MEDMEM_Med.hxx"
   9 #include "MEDMEM_Family.hxx"
  10 #include "MEDMEM_Field.hxx"
  11 #include "MEDMEM_Group.hxx"
  12 #include "MEDMEM_Coordinate.hxx"
  13 #include "MEDMEM_Connectivity.hxx"
  14 #include "MEDMEM_Mesh.hxx"
  15 #include "MEDMEM_CellModel.hxx"
  16 #include "MEDMEM_define.hxx"
  17 #include "MEDMEM_DriverTools.hxx"
  18
  19 #include <stdio.h>
  20 #include <fcntl.h>
  21 #ifdef WNT
  22 #include <io.h>
  23 #else
  24 #include <unistd.h>
  25 #endif
  26
  27 #include <float.h>
  28
  29 /////
  30 using namespace std;
  31 using namespace MED_EN;
  32 using namespace MEDMEM;
  33 /////
  34
  35
  36 // Every memory allocation made in the MedDriver members function are desallocated in the Mesh destructor
  37
  38 /////
  39 const size_t GIBI_MESH_DRIVER::nb_geometrie_gibi;
  40
  41 const medGeometryElement GIBI_MESH_DRIVER::geomGIBItoMED[nb_geometrie_gibi] =
  42      {   /*1 */ MED_POINT1 ,/*2 */ MED_SEG2   ,/*3 */ MED_SEG3   ,/*4 */ MED_TRIA3  ,/*5 */ MED_NONE   ,
  43        /*6 */ MED_TRIA6  ,/*7 */ MED_NONE   ,/*8 */ MED_QUAD4  ,/*9 */ MED_NONE   ,/*10*/ MED_QUAD8  ,
  44        /*11*/ MED_NONE   ,/*12*/ MED_NONE   ,/*13*/ MED_NONE   ,/*14*/ MED_HEXA8  ,/*15*/ MED_HEXA20 ,
  45        /*16*/ MED_PENTA6 ,/*17*/ MED_PENTA15,/*18*/ MED_NONE   ,/*19*/ MED_NONE   ,/*20*/ MED_NONE   ,
  46        /*21*/ MED_NONE   ,/*22*/ MED_NONE   ,/*23*/ MED_TETRA4 ,/*24*/ MED_TETRA10,/*25*/ MED_PYRA5  ,
  47        /*26*/ MED_PYRA13 ,/*27*/ MED_NONE   ,/*28*/ MED_NONE   ,/*29*/ MED_NONE   ,/*30*/ MED_NONE   ,
  48        /*31*/ MED_NONE   ,/*32*/ MED_NONE   ,/*33*/ MED_NONE   ,/*34*/ MED_NONE   ,/*35*/ MED_NONE   ,
  49        /*36*/ MED_NONE   ,/*37*/ MED_NONE   ,/*38*/ MED_NONE   ,/*39*/ MED_NONE   ,/*40*/ MED_NONE   ,
  50        /*41*/ MED_NONE   ,/*42*/ MED_NONE   ,/*43*/ MED_NONE   ,/*44*/ MED_NONE   ,/*45*/ MED_NONE   ,
  51        /*46*/ MED_NONE   ,/*47*/ MED_NONE   };
  52
  53 //=======================================================================
  54 //function : gibi2medGeom
  55 //purpose  :
  56 //=======================================================================
  57
  58 medGeometryElement GIBI_MESH_DRIVER::gibi2medGeom( size_t gibiTypeNb )
  59 {
  60   if ( gibiTypeNb < 1 || gibiTypeNb > 47 )
  61     return MED_NONE;
  62
  63   return geomGIBItoMED[ gibiTypeNb - 1 ];
  64 }
  65
  66 //=======================================================================
  67 //function : med2gibiGeom
  68 //purpose  :
  69 //=======================================================================
  70
  71 int GIBI_MESH_DRIVER::med2gibiGeom( medGeometryElement medGeomType )
  72 {
  73   for ( int i = 0; i < nb_geometrie_gibi; i++ )
  74     if ( geomGIBItoMED[ i ] == medGeomType )
  75       return i + 1;
  76
  77   return -1;
  78 }
  79
  80 //=======================================================================
  81 //function : getGroupId
  82 //purpose  :
  83 //=======================================================================
  84
  85 static int getGroupId(const vector<int>& support_ids, _intermediateMED*  medi)
  86 {
  87   int group_id = 0;
  88   vector<int>::const_iterator sb = support_ids.begin(), se = support_ids.end();
  89   if (support_ids.size() == 1 || // one or equal support ids
  90       *std::max_element( sb, se ) == *std::min_element( sb, se ))
  91     group_id = support_ids[0] - 1;
  92   else {
  93     // try to find an existing group with the same sub-groups
  94     set<int> sup_set;
  95     for ( ; sb != se; sb++ )
  96       sup_set.insert( *sb );
  97     for ( group_id = 0; group_id < medi->groupes.size(); ++group_id ) {
  98       if (sup_set.size() == medi->groupes[ group_id ].groupes.size() &&
  99           std::equal (sup_set.begin(), sup_set.end(),
 100                       medi->groupes[ group_id ].groupes.begin()))
 101         break;
 102     }
 103     if ( group_id == medi->groupes.size() )
 104     {
 105       // no such a group, add a new one
 106       medi->groupes.push_back( _groupe() );
 107       _groupe& new_grp = medi->groupes.back();
 108       //new_grp.nom = string( group_id % 10 + 1, 'G' );
 109       new_grp.groupes.reserve( sup_set.size() );
 110       for ( set<int>::iterator it = sup_set.begin(); it != sup_set.end(); it++ ) {
 111         new_grp.groupes.push_back( *it );
 112         //new_grp.nom += "_" + medi->groupes[ *it - 1 ].nom;
 113       }
 114     }
 115   }
 116   return group_id;
 117 }
 118
 119 //=======================================================================
 120 //function : read
 121 //purpose  :
 122 //=======================================================================
 123
 124 #define GIBI_EQUAL(var_str, stat_str) \
 125   (strncmp (var_str, stat_str, sizeof(stat_str)-1) == 0)
 126 #define DUMP_LINE_NB " on line " << _lineNb
 127
 128 bool GIBI_MESH_RDONLY_DRIVER::readFile (_intermediateMED* medi, bool readFields )
 129 {
 130   const char * LOC = "GIBI_MESH_RDONLY_DRIVER::readFile() : " ;
 131   BEGIN_OF(LOC);
 132
 133   // LECTURE DES DONNEES DS FICHIER GIBI
 134
 135   enum Readable_Piles {
 136     PILE_SOUS_MAILLAGE=1,
 137     PILE_NODES_FIELD  =2,
 138     PILE_NOEUDS       =32,
 139     PILE_COORDONNEES  =33,
 140     PILE_FIELD        =39,
 141     PILE_LAST_READABLE=39
 142     };
 143   Readable_Piles readable_Piles [] = {
 144     PILE_SOUS_MAILLAGE,
 145     PILE_NODES_FIELD,
 146     PILE_NOEUDS,
 147     PILE_COORDONNEES,
 148     PILE_FIELD,
 149     PILE_LAST_READABLE
 150     };
 151   char* ligne; // pour lire une ligne
 152   const char* enregistrement_type=" ENREGISTREMENT DE TYPE";
 153   vector<int> numero_noeuds; // tableau de travail (indices)
 154   set<int> donePiles; // already read piles
 155   unsigned space_dimension = 0;
 156
 157   while ( getNextLine(ligne, false)) // boucle externe de recherche de "ENREGISTREMENT DE TYPE"
 158   {
 159     if ( !GIBI_EQUAL( ligne, enregistrement_type ))
 160       continue; // "ENREGISTREMENT DE TYPE" non trouvé -> on lit la ligne suivante
 161
 162     // lecture du numéro d'enregistrement
 163     int numero_enregistrement = atoi( ligne + strlen(enregistrement_type) + 1 );
 164
 165     enum { ENREG_TYPE_2=2, ENREG_TYPE_4=4}; // énumération des types d'enregistrement traités
 166     int numero_pile, nb_objets_nommes, nb_objets, nb_indices;
 167     vector<int> indices_objets_nommes;
 168     vector<string> objets_nommes;
 169
 170     if (numero_enregistrement == ENREG_TYPE_4)
 171     {
 172       getNextLine(ligne);
 173       const char* s = " NIVEAU  15 NIVEAU ERREUR   0 DIMENSION";
 174       space_dimension = atoi( ligne + strlen( s ) + 1 );
 175       if ( !GIBI_EQUAL( ligne, " NIVEAU" ) || space_dimension < 1 ) {
 176         INFOS( " Could not read file: syntax error in type 4 record");
 177         return false;
 178       }
 179     }
 180     else if (numero_enregistrement == ENREG_TYPE_2 )
 181     {
 182       if ( space_dimension == 0 ) {
 183         INFOS( "Missing ENREGISTREMENT DE TYPE   4");
 184         return false;
 185       }
 186       // FORMAT(' PILE NUMERO',I4,'NBRE OBJETS NOMMES',I8,'NBRE OBJETS',I8)
 187       getNextLine(ligne);
 188       const char *s1 = " PILE NUMERO", *s2 = "NBRE OBJETS NOMMES", *s3 = "NBRE OBJETS";
 189       if ( ! GIBI_EQUAL( ligne, s1 ) ) {
 190         INFOS( " Could not read file: error in type 2 record. " << ligne);
 191         return false;
 192       }
 193       ligne = ligne + strlen(s1);
 194       numero_pile = atoi( ligne );
 195       ligne = ligne + 4 + strlen(s2);
 196       nb_objets_nommes = atoi( ligne );
 197       ligne = ligne + 8 + strlen(s3);
 198       nb_objets = atoi( ligne );
 199       if ( nb_objets_nommes<0 || nb_objets<0  ) {
 200         INFOS(" Could not read file: " << nb_objets << " " <<nb_objets_nommes);
 201         return false;
 202       }
 203       if ( !donePiles.insert( numero_pile ).second ) // piles may repeat
 204         continue;
 205
 206       if ( numero_pile > PILE_LAST_READABLE )
 207         break; // stop file reading
 208
 209       // skip not readable piles
 210       int i = -1;
 211       while ( readable_Piles[ ++i ] != PILE_LAST_READABLE )
 212         if ( readable_Piles[ i ] == numero_pile )
 213           break;
 214       if ( readable_Piles[ i ] != numero_pile )
 215         continue;
 216
 217       // lecture des objets nommés et de leurs indices
 218       objets_nommes.resize(nb_objets_nommes);
 219       indices_objets_nommes.resize(nb_objets_nommes);
 220       for ( initNameReading( nb_objets_nommes ); more(); next() ) {
 221         objets_nommes[ index() ] = getName();
 222       }
 223       for ( initIntReading( nb_objets_nommes ); more(); next() )
 224         indices_objets_nommes[ index() ] = getInt();
 225
 226       // boucle interne : lecture de la pile
 227
 228       MESSAGE(LOC << "---- Traitement pile " << numero_pile);
 229
 230       // -----------------------------------
 231       //                        MESH GROUPS
 232       // -----------------------------------
 233
 234       if (numero_pile == PILE_SOUS_MAILLAGE )
 235       {
 236         map<int,int> strangeGroupType;
 237         medi->groupes.reserve(nb_objets*2); // fields may add some groups
 238         for (int objet=0; objet!=nb_objets; ++objet) // pour chaque groupe
 239         {
 240           initIntReading( 5 );
 241           unsigned type_geom_castem = getInt(); next();
 242           unsigned nb_sous_maillage = getInt(); next();
 243           unsigned nb_reference     = getInt(); next();
 244           unsigned nb_noeud         = getInt(); next();
 245           unsigned nb_elements      = getInt();
 246
 247           // le cas type_geom_castem=0 correspond aux maillages composites
 248           if (type_geom_castem<0) {
 249             INFOS(" Error while reading file, bad geometric type:" << type_geom_castem);
 250             return false;
 251           }
 252
 253           medi->groupes.push_back(_groupe());
 254           _groupe & groupe = medi->groupes.back();
 255
 256           // si le groupe se compose de sous-maillages (ie groupe composite)
 257           if (type_geom_castem==0 && nb_sous_maillage>0)
 258           {
 259             // lecture des indices des sous-maillages, stockage.
 260             // les mailles correspondant a ces sous_maillages seront inserees a la fin du case
 261             groupe.groupes.resize( nb_sous_maillage );
 262             for ( initIntReading( nb_sous_maillage ); more(); next() ) {
 263               groupe.groupes[ index() ] = getInt();
 264             }
 265             if ( readFields )
 266               std::sort( groupe.groupes.begin(), groupe.groupes.end() );
 267           }
 268           // lecture des references (non utilisé pour MED)
 269           for ( i = 0; i < nb_reference; i += 10 ) {// FORMAT(10I8)
 270             getNextLine(ligne);
 271           }
 272           // lecture des couleurs (non utilisé pour MED)
 273           for ( i = 0; i < nb_elements; i += 10 ) {
 274             getNextLine(ligne);
 275           }
 276           // not a composit group
 277           if (type_geom_castem>0 && nb_sous_maillage==0)
 278           {
 279             medGeometryElement medType = gibi2medGeom(type_geom_castem);
 280             bool goodType = ( medType!=MED_NONE );
 281             if ( !goodType )
 282               strangeGroupType.insert( make_pair( objet, type_geom_castem ));
 283
 284             pair<set<_maille>::iterator,bool> p;
 285             pair<map<int,_noeud>::iterator,bool> p_no;
 286             _noeud no;
 287             no.coord.resize(space_dimension);
 288             _maille ma( medType, nb_noeud );
 289             ma.sommets.resize(nb_noeud);
 290             if ( goodType )
 291               groupe.mailles.resize( nb_elements );
 292
 293             // lecture pour chaque maille des sommets et insertions
 294             initIntReading( nb_elements * nb_noeud );
 295             if ( !goodType ) {
 296               while ( more() )
 297                 next();
 298             }
 299             else {
 300               for ( i = 0; i < nb_elements; ++i )
 301               {
 302                 for (unsigned n = 0; n < nb_noeud; ++n, next() )
 303                 {
 304                   if ( !more() ) {
 305                     INFOS( " Error while reading elem nodes ");
 306                     return false;
 307                   }
 308                   no.number = getInt();
 309                   p_no=medi->points.insert(make_pair(no.number, no));
 310                   ma.sommets[n]=p_no.first;
 311                 }
 312                 p=medi->maillage.insert(ma);
 313                 groupe.mailles[i] = p.first; // on stocke dans le groupe un iterateur sur la maille
 314               }
 315             }
 316           }
 317         } // loop on groups
 318
 319         // set group names
 320         for (i=0; i!=nb_objets_nommes; ++i) {
 321           int grpID = indices_objets_nommes[i];
 322           _groupe & grp = medi->groupes[ grpID-1 ];
 323           if ( !grp.nom.empty() ) // a group has several names
 324           { // create a group with subgroup grp and named grp.nom
 325             medi->groupes.push_back(_groupe());
 326             medi->groupes.back().groupes.push_back( grpID );
 327             medi->groupes.back().nom = grp.nom;
 328           }
 329           grp.nom=objets_nommes[i];
 330           map<int,int>::iterator it = strangeGroupType.find( grpID - 1 );
 331           if ( it != strangeGroupType.end() ) {
 332             //INFOS( "Skip " << grp.nom << " of not supported CASTEM type: " << it->second );
 333           }
 334         }
 335
 336       }// Fin case PILE_SOUS_MAILLAGE
 337
 338       // ---------------------------------
 339       //                            NODES
 340       // ---------------------------------
 341
 342       else if ( numero_pile == PILE_NOEUDS )
 343       {
 344         getNextLine( ligne );
 345         std::vector<int> place_noeuds;
 346         nb_indices = atoi ( ligne );
 347         if (nb_indices != nb_objets)
 348         {
 349           INFOS("Erreur de lecture dans enregistrement de pile " << PILE_NOEUDS);
 350           return false;
 351         }
 352
 353         place_noeuds.resize(nb_objets);
 354         for ( initIntReading( nb_objets ); more(); next() )
 355           place_noeuds[ index() ] = getInt();
 356         int max=(* std::max_element(place_noeuds.begin(),place_noeuds.end()));
 357
 358         // numero_noeuds contient pour chacun des max noeuds qu'on va lire dans le case PILE_COORDONNEES
 359         // son indice dans la connectivite du maillage. Cet indice correspond egalement a la cle du map
 360         // medi->points ou l'on stocke les noeuds.
 361         numero_noeuds.resize(max,-1);
 362         for (unsigned i=0; i!=place_noeuds.size(); ++i)
 363           numero_noeuds[place_noeuds[i]-1]=i+1;
 364       }
 365
 366       // ---------------------------------------
 367       //                            COORDINATES
 368       // ---------------------------------------
 369
 370       else if ( numero_pile == PILE_COORDONNEES )
 371       {
 372         getNextLine( ligne );
 373         unsigned nb_reels = atoi( ligne );
 374         // PROVISOIRE : certains fichier gibi n`ont
 375         if (nb_reels < numero_noeuds.size()*(space_dimension)) {
 376           INFOS("Erreur de lecture dans enregistrement de pile " << PILE_COORDONNEES);
 377           return false;
 378         }
 379         initDoubleReading( nb_reels );
 380         map< int, _noeud >::iterator pIt;
 381         for (unsigned i=0; i!=numero_noeuds.size(); ++i)
 382         {
 383           // si le noeud est utilisé dans le maillage,
 384           //on lit ses coordonnées et on les stocke dans la structure
 385           if (( numero_noeuds[i] != -1 ) &&
 386               (( pIt = medi->points.find(numero_noeuds[i])) != medi->points.end()))
 387           {
 388             for (unsigned j=0; j!=space_dimension; ++j, next())
 389               pIt->second.coord[j] = getDouble();
 390             next(); // on ne conserve pas la densite
 391           }
 392           else // sinon, on passe au noeud suivant
 393           {
 394             for (unsigned j=0; j!=space_dimension+1; ++j)
 395               next();
 396           }
 397         }
 398       }
 399
 400       // ---------------------------------------
 401       //                            NODE FIELDS
 402       // ---------------------------------------
 403
 404       else if ( numero_pile == PILE_NODES_FIELD && readFields )
 405       {
 406         vector< _fieldBase* > fields( nb_objets );
 407         for (int objet=0; objet!=nb_objets; ++objet) // pour chaque field
 408         {
 409           initIntReading( 4 ); // nb subcomponents, nb components, IFOUR, nb attributes
 410           int i_sub,  nb_sub  = getInt(); next();
 411           int i_comp, nb_comp = getInt(); next(); next();
 412           int nb_attr = getInt();
 413           if ( nb_sub < 0 || nb_comp < 0 ) {
 414             INFOS("Error of field reading: wrong nb of components "
 415                   << nb_sub << " " << nb_comp << DUMP_LINE_NB);
 416             return false;
 417           }
 418           _field<double>* fdouble = 0;
 419           if ( nb_sub > 0 && nb_comp > 0 ) {
 420             fdouble = new _field<double>(MED_REEL64);
 421             medi->fields.push_back( fdouble );
 422             fdouble->nb_subcomponents = nb_sub;
 423             fdouble->nb_components    = nb_comp;
 424             fdouble->comp_names.resize( nb_comp );
 425           }
 426           fields[ objet ] = fdouble;
 427
 428           // loop on subcomponents of a field, each of which refers to
 429           // a certain support and has the same number of components
 430           vector<int>     support_ids( nb_sub );
 431           vector<int>     nb_values( nb_sub );
 432           vector<string>  comp_names( nb_comp );
 433           initIntReading( nb_sub * 3 );
 434           for ( i_sub = 0; i_sub < nb_sub; ++i_sub )
 435           {
 436             support_ids[ i_sub ] = -getInt(); next(); // reference to support
 437             if ( support_ids[ i_sub ] < 1 || support_ids[ i_sub ] > medi->groupes.size() ) {
 438               INFOS("Error of field reading: wrong mesh reference "<< support_ids[ i_sub ]);
 439               return false;
 440             }
 441             nb_values[ i_sub ] = getInt(); next();    // nb points
 442             if ( nb_values[ i_sub ] < 0 ) {
 443               INFOS(" Wrong nb of points: " << nb_values[ i_sub ] );
 444               return false;
 445             }
 446             next(); // skip nb of comp names equal to nb_comp
 447           }
 448           initNameReading( nb_sub * nb_comp, 4 );
 449           for ( i_sub = 0; i_sub < nb_sub; ++i_sub )
 450           {
 451             for ( i_comp = 0; i_comp < nb_comp; ++i_comp, next() ) {
 452               ASSERT( more() );
 453               fdouble->comp_names[ i_comp ] = getName();// component names
 454             }
 455           }
 456
 457           for ( initIntReading( nb_sub ); more(); next() )// nb harmonics ( ignored )
 458             ;
 459
 460           getNextLine( ligne );                             // TYPE ( ignored )
 461           getNextLine( ligne );                             // TITRE ( ignored )
 462           for ( initIntReading( nb_attr ); more(); next() ) // attributes ( ignored )
 463             ;
 464
 465           for ( i_sub = 0; i_sub < nb_sub; ++i_sub )
 466           {
 467             // loop on components: read values
 468             initDoubleReading( nb_values[ i_sub ] * nb_comp );
 469             for ( i_comp = 0; i_comp < nb_comp; ++i_comp ) {
 470               vector<double> & vals =
 471                 fdouble->addComponent( support_ids[ i_sub ], nb_values[ i_sub ] );
 472               for ( int i = 0; more() && i < nb_values[ i_sub ]; next(), ++i ) {
 473                 vals[ i ] = getDouble();
 474               }
 475             }
 476           } // loop on subcomponents of a field
 477
 478           // supporting group id
 479           if ( fdouble )
 480             fdouble->group_id = getGroupId( support_ids, medi );
 481
 482         } // end loop on field objects
 483
 484         // set field names
 485         for ( i = 0; i < nb_objets_nommes; ++i ) {
 486           int fieldIndex = indices_objets_nommes[ i ];
 487           if ( fields[ fieldIndex - 1 ] )
 488             fields[ fieldIndex - 1 ]->name = objets_nommes[ i ];
 489         }
 490
 491       }  // Fin numero_pile == PILE_NODES_FIELD
 492
 493       // -------------------------------------------------
 494       //                                           FIELDS
 495       // -------------------------------------------------
 496
 497       else if ( numero_pile == PILE_FIELD && readFields )
 498       {
 499         vector< _fieldBase* > fields( nb_objets, (_fieldBase*)NULL );
 500         for (int objet=0; objet!=nb_objets; ++objet) // pour chaque field
 501         {
 502           initIntReading( 1 );
 503           int i_sub, nb_sub = getInt(); // <nb_sub> -1 6 <title length>
 504           if ( nb_sub < 1 ) {
 505             INFOS("Error of field reading: wrong nb of subcomponents " << nb_sub);
 506             return false;
 507           }
 508           getNextLine( ligne ); // title
 509           // look for a line starting with '-' : <reference to support>
 510           do {
 511             initIntReading( nb_sub * 9 );
 512           } while ( getInt() >= 0 );
 513
 514           int i_comp, nb_comp = 0;
 515           vector<int> support_ids( nb_sub );
 516           for ( i_sub = 0; i_sub < nb_sub; ++i_sub )
 517           {
 518             support_ids[ i_sub ] = -getInt(); next(); // <reference to support>
 519             next();                                   // skip <address>
 520             nb_comp = getInt(); next();               // <nb of components in the field>
 521             for ( i = 0; i < 6; ++i )                 //  0 0 0 -1 0 2
 522               next();
 523             if ( support_ids[ i_sub ] < 1 || support_ids[ i_sub ] > medi->groupes.size() ) {
 524               INFOS("Error of field reading: wrong mesh reference "<< support_ids[ i_sub ]);
 525               return false;
 526             }
 527             if ( nb_comp < 1 ) {
 528               INFOS("Error of field reading: wrong nb of components " << nb_comp);
 529               return false;
 530             }
 531           }
 532           for ( initNameReading( nb_sub, 17 ); more(); next() )
 533             ; // dummy strings
 534           for ( initNameReading( nb_sub ); more(); next() )
 535             ; // dummy strings
 536
 537           // loop on subcomponents of a field, each of which refers to
 538           // a certain support and has the same number of components
 539           _field<double>* fdouble = 0;
 540           _field<int>*    fint    = 0;
 541           vector<string> comp_names( nb_comp ), comp_type( nb_comp );
 542           for ( i_sub = 0; i_sub < nb_sub; ++ i_sub )
 543           {
 544             getNextLine( ligne );  // skip nb_comp addresses of MELVAL structure
 545
 546             // component names
 547             for ( initNameReading( nb_comp ); more(); next() )
 548               comp_names[ index() ] = getName();
 549
 550             // component type
 551             for ( initNameReading( nb_comp, 17 ); more(); next() ) { // 17 is name width
 552               comp_type[ index() ] = getName();
 553               // component types must be the same
 554               if ( index() > 0 && comp_type[ index() ] != comp_type[ index() - 1] ) {
 555                 INFOS( "Error of field reading: diff component types <"
 556                       << comp_type[ index() ] << "> != <" << comp_type[ index() - 1 ] << ">");
 557                 return false;
 558               }
 559             }
 560             // now type is known, create a field
 561             bool isReal = ( comp_type[0] == "REAL*8" );
 562             if ( !fdouble && !fint ) {
 563               _fieldBase * fbase;
 564               if ( !isReal ) {
 565                 fbase = fint = new _field<int>(MED_INT32);
 566                 INFOS( "Warning: read NOT REAL field, type <" << comp_type[0] << ">"
 567                       << DUMP_LINE_NB);
 568               }
 569               else
 570                 fbase = fdouble = new _field<double>(MED_REEL64);
 571               medi->fields.push_back( fbase ); // medi->fields is a std::list
 572               fields[ objet ]         = fbase;
 573               fbase->nb_subcomponents = nb_sub;
 574               fbase->nb_components    = nb_comp;
 575               fbase->comp_names       = comp_names;
 576               fbase->group_id         = getGroupId( support_ids, medi );
 577             }
 578
 579             // loop on components: read values
 580             for ( i_comp = 0; i_comp < nb_comp; ++i_comp ) {
 581               // nb of values
 582               initIntReading( 2 );
 583               int nb_val_by_elem = getInt(); next();
 584               int nb_values      = getInt();
 585               if ( nb_val_by_elem != 1 ) {
 586                 INFOS("Error of field reading: nb of values by element "
 587                       << " != 1 : " << nb_val_by_elem << DUMP_LINE_NB );
 588                 return false;
 589               }
 590               if ( isReal ) {
 591                 vector<double> & vals =
 592                   fdouble->addComponent( support_ids[ i_sub ], nb_values );
 593                 for ( initDoubleReading( nb_values ); more(); next()) {
 594                   vals[ index() ] = getDouble();
 595                 }
 596               }
 597               else {
 598                 vector<int> & vals =
 599                   fint->addComponent( support_ids[ i_sub ], nb_values );
 600                 for ( initIntReading( nb_values ); more(); next() ) {
 601                   vals[ index() ] = getInt();
 602                 }
 603               }
 604             }
 605           } // loop on subcomponents of a field
 606
 607         } // end loop on field objects
 608
 609         // set field names
 610         for ( i = 0; i < nb_objets_nommes; ++i ) {
 611           int fieldIndex = indices_objets_nommes[ i ];
 612           fields[ fieldIndex - 1 ]->name = objets_nommes[ i ];
 613         }
 614
 615       } // numero_pile == PILE_FIELD && readFields
 616
 617       else if ( numero_pile >= PILE_LAST_READABLE )
 618         break; // stop file reading
 619
 620     } // Fin case ENREG_TYPE_2
 621   } //  fin de la boucle while de lecture externe
 622
 623   // check if all needed piles present
 624   if ( donePiles.find( PILE_SOUS_MAILLAGE ) != donePiles.end() )
 625   {
 626     if (donePiles.find( PILE_NOEUDS ) == donePiles.end() ) {
 627       INFOS( " Missing pile " << PILE_NOEUDS );
 628       return false;
 629     }
 630     if (donePiles.find( PILE_COORDONNEES ) == donePiles.end()) {
 631       INFOS( " Missing pile " << PILE_COORDONNEES );
 632       return false;
 633     }
 634   }
 635
 636   END_OF(LOC);
 637   return true;
 638 }
 639
 640 GIBI_MESH_DRIVER::GIBI_MESH_DRIVER():
 641        GENDRIVER(),
 642        _ptrMesh(( MESH *)MED_NULL),
 643        // A VOIR _medIdt(MED_INVALID),
 644        _meshName("")
 645 {
 646   MESSAGE("GIBI_MESH_DRIVER()");
 647 }
 648
 649 GIBI_MESH_DRIVER::GIBI_MESH_DRIVER(const string & fileName,
 650                                    MESH * ptrMesh,
 651                                    MED_EN::med_mode_acces accessMode):
 652   GENDRIVER(fileName,accessMode),
 653   _ptrMesh(ptrMesh)
 654   // A VOIR _medIdt(MED_INVALID),
 655 {
 656   MESSAGE( "GIBI_MESH_DRIVER(" << fileName <<","<<accessMode );
 657 //   _meshName=fileName.substr(0,fileName.rfind("."));
 658     // mesh name construction from fileName
 659     const string ext=".sauve"; // expected extension
 660     string::size_type pos=fileName.find(ext,0);
 661     string::size_type pos1=fileName.rfind('/');
 662     _meshName = string(fileName,pos1+1,pos-pos1-1); //get rid of directory & extension
 663     SCRUTE(_meshName);
 664 }
 665
 666 GIBI_MESH_DRIVER::GIBI_MESH_DRIVER(const GIBI_MESH_DRIVER & driver):
 667   GENDRIVER(driver),
 668   _ptrMesh(driver._ptrMesh),
 669   // A VOIR _medIdt(MED_INVALID),
 670   _meshName(driver._meshName)
 671 {
 672   MESSAGE("GIBI_MESH_DRIVER(const GIBI_MESH_DRIVER & driver)");
 673 }
 674
 675 GIBI_MESH_DRIVER::~GIBI_MESH_DRIVER()
 676 {
 677   MESSAGE("~GIBI_MESH_DRIVER()");
 678 }
 679 void    GIBI_MESH_DRIVER::setMeshName(const string & meshName) { _meshName = meshName; };
 680 string  GIBI_MESH_DRIVER::getMeshName() const { return _meshName; };
 681
 682
 683 //---------------------------------- RDONLY PART -------------------------------------------------------------
 684
 685 GIBI_MESH_RDONLY_DRIVER::GIBI_MESH_RDONLY_DRIVER():
 686        GIBI_MESH_DRIVER(),
 687        _File (-1),_start(0L),_ptr  (0L),_eptr (0L)
 688 {
 689 }
 690 GIBI_MESH_RDONLY_DRIVER::GIBI_MESH_RDONLY_DRIVER(const string & fileName,MESH * ptrMesh):
 691        GIBI_MESH_DRIVER(fileName,ptrMesh,MED_RDONLY),
 692        _File (-1),_start(0L),_ptr  (0L),_eptr (0L)
 693 {
 694     MESSAGE("GIBI_MESH_RDONLY_DRIVER::GIBI_MESH_RDONLY_DRIVER"
 695             "(const string & fileName, MESH * ptrMesh) has been created, "
 696             << fileName << ", " << MED_RDONLY);
 697 }
 698 GIBI_MESH_RDONLY_DRIVER::GIBI_MESH_RDONLY_DRIVER(const GIBI_MESH_RDONLY_DRIVER & driver):
 699 GIBI_MESH_DRIVER(driver)
 700 {
 701 }
 702 GIBI_MESH_RDONLY_DRIVER::~GIBI_MESH_RDONLY_DRIVER()
 703 {
 704   BEGIN_OF( "~GIBI_MESH_RDONLY_DRIVER()");
 705   if (_File >= 0)
 706   {
 707     ::close (_File);
 708     if (_start != 0L)
 709       delete [] _start;
 710   }
 711   MESSAGE("GIBI_MESH_RDONLY_DRIVER::~GIBI_MESH_RDONLY_DRIVER() has been destroyed");
 712 }
 713 GENDRIVER * GIBI_MESH_RDONLY_DRIVER::copy(void) const
 714 {
 715   return new GIBI_MESH_RDONLY_DRIVER(*this);
 716 }
 717
 718 //=======================================================================
 719 //function : open
 720 //purpose  :
 721 //=======================================================================
 722
 723 const int GIBI_MaxOutputLen   = 150;
 724 const int GIBI_BufferSize     = 16184; // for non-stream input
 725
 726 void GIBI_MESH_RDONLY_DRIVER::open()
 727   //     throw (MEDEXCEPTION)
 728 {
 729   const char * LOC = "GIBI_MESH_RDONLY_DRIVER::open()" ;
 730   BEGIN_OF(LOC);
 731
 732 //   MED_EN::med_mode_acces aMode = getAccessMode();
 733 //   if ( aMode != MED_EN::MED_LECT && aMode != MED_EN::MED_REMP )
 734 //     throw MEDEXCEPTION(LOCALIZED(STRING(LOC) << " Bad file mode access ! " << aMode));
 735
 736 #ifdef WNT
 737   _File = ::_open (_fileName.c_str(), _O_RDONLY|_O_BINARY);
 738 #else
 739   _File = ::open (_fileName.c_str(), O_RDONLY);
 740 #endif
 741   if (_File >= 0)
 742   {
 743     _start = new char [GIBI_BufferSize];
 744     _ptr   = _start;
 745     _eptr  = _start;
 746     _status = MED_OPENED;
 747     _lineNb = 0;
 748   }
 749   else
 750   {
 751     _status = MED_CLOSED;
 752     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" Could not open file "<<_fileName
 753                                  << " fd: " << _File));
 754   }
 755   END_OF(LOC);
 756 }
 757
 758 //=======================================================================
 759 //function : close
 760 //purpose  :
 761 //=======================================================================
 762
 763 void GIBI_MESH_RDONLY_DRIVER::close()
 764 {
 765     const char * LOC = "GIBI_MESH_DRIVER::close() " ;
 766     BEGIN_OF(LOC);
 767     if ( _status == MED_OPENED)
 768     {
 769       if (_File >= 0)
 770       {
 771         ::close (_File);
 772         if (_start != 0L)
 773           delete [] _start;
 774         _File = -1;
 775       }
 776       _status = MED_CLOSED;
 777     }
 778     END_OF(LOC);
 779 }
 780
 781 //=======================================================================
 782 //function : getLine
 783 //purpose  :
 784 //=======================================================================
 785
 786 bool GIBI_MESH_RDONLY_DRIVER::getLine(char* & aLine)
 787 {
 788   bool aResult = true;
 789   // Check the state of the buffer;
 790   // if there is too little left, read the next portion of data
 791   int nBytesRest = _eptr - _ptr;
 792   if (nBytesRest < GIBI_MaxOutputLen)
 793   {
 794     if (nBytesRest > 0) {
 795       memcpy (_start, _ptr, nBytesRest);
 796     } else
 797       nBytesRest = 0;
 798     _ptr = _start;
 799     const int nBytesRead = ::read (_File,
 800                                    &_start [nBytesRest],
 801                                    GIBI_BufferSize - nBytesRest);
 802     nBytesRest += nBytesRead;
 803     _eptr = &_start [nBytesRest];
 804   }
 805   // Check the buffer for the end-of-line
 806   char * ptr = _ptr;
 807   while (~0)
 808   {
 809     // Check for end-of-the-buffer, the ultimate criterion for termination
 810     if (ptr >= _eptr)
 811     {
 812       if (nBytesRest <= 0)
 813         aResult = false;
 814       else
 815         _eptr[-1] = '\0';
 816       break;
 817     }
 818     // seek the line-feed character
 819     if (ptr[0] == '\n')
 820     {
 821       if (ptr[-1] == '\r')
 822         ptr[-1] = '\0';
 823       ptr[0] = '\0';
 824       ++ptr;
 825       break;
 826     }
 827     ++ptr;
 828   }
 829   // Output the result
 830   aLine = _ptr;
 831   _ptr = ptr;
 832   _lineNb++;
 833
 834   return aResult;
 835 }
 836
 837 //=======================================================================
 838 //function : init
 839 //purpose  :
 840 //=======================================================================
 841
 842 void GIBI_MESH_RDONLY_DRIVER::init( int nbToRead, int nbPosInLine, int width, int shift )
 843 {
 844   _nbToRead = nbToRead;
 845   _nbPosInLine = nbPosInLine;
 846   _width = width;
 847   _shift = shift;
 848   _iPos = _iRead = 0;
 849   if ( _nbToRead ) {
 850     getNextLine( _curPos );
 851     _curPos = _curPos + _shift;
 852   }
 853   else
 854     _curPos = 0;
 855 }
 856
 857 //=======================================================================
 858 //function : next
 859 //purpose  : line getting
 860 //=======================================================================
 861
 862 void GIBI_MESH_RDONLY_DRIVER::next()
 863 {
 864   if ( !more() ) throw MEDEXCEPTION(LOCALIZED("!more()"));
 865   ++_iRead;
 866   ++_iPos;
 867   if ( _iRead < _nbToRead ) {
 868     if ( _iPos >= _nbPosInLine ) {
 869       getNextLine( _curPos );
 870       _curPos = _curPos + _shift;
 871       _iPos = 0;
 872     }
 873     else
 874       _curPos = _curPos + _width + _shift;
 875   }
 876   else
 877     _curPos = 0;
 878 }
 879
 880 //=======================================================================
 881 //function : getName
 882 //purpose  : names reading
 883 //=======================================================================
 884
 885 string GIBI_MESH_RDONLY_DRIVER::getName() const
 886 {
 887   int len = _width;
 888   while (( _curPos[len-1] == ' ' || _curPos[len-1] == 0) && len > 0 )
 889     len--;
 890   return string( _curPos, len );
 891 }
 892
 893 //=======================================================================
 894 //function : read
 895 //purpose  :
 896 //=======================================================================
 897
 898 void GIBI_MESH_RDONLY_DRIVER::read(void) throw (MEDEXCEPTION)
 899 {
 900   const char * LOC = "_GIBI_RDONLY_DRIVER::read() : " ;
 901
 902   if (_status!=MED_OPENED)
 903     throw MEDEXCEPTION(LOCALIZED(STRING(LOC) << "The _idt of file " << _fileName << " is : "
 904                                  <<  " (the file is not opened)." )) ;
 905   if ( ! _ptrMesh->isEmpty() )
 906     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Mesh object not empty : can't fill it!"));
 907
 908   _intermediateMED medi;
 909   try {
 910     if ( readFile( &medi, false )) {
 911       // impression résultats
 912       MESSAGE(LOC << "GIBI_MESH_RDONLY_DRIVER::read : RESULTATS STRUCTURE INTERMEDIAIRES : ");
 913       MESSAGE(LOC <<  medi );
 914
 915       fillMesh( &medi );
 916     }
 917   }
 918   catch (MEDEXCEPTION &ex)
 919   {
 920     INFOS( ex.what() );
 921   }
 922
 923 }
 924
 925 //=======================================================================
 926 //function : getReverseVector
 927 //purpose  :
 928 //=======================================================================
 929
 930 static void getReverseVector (const medGeometryElement type,
 931                               vector<pair<int,int> > & swapVec )
 932 {
 933   BEGIN_OF("void getReverseVector()");
 934   swapVec.clear();
 935
 936   switch ( type ) {
 937   case MED_TETRA4:
 938     swapVec.resize(1);
 939     swapVec[0] = make_pair( 1, 2 );
 940     break;
 941   case MED_PYRA5:
 942     swapVec.resize(1);
 943     swapVec[0] = make_pair( 1, 3 );
 944     break;
 945   case MED_PENTA6:
 946     swapVec.resize(2);
 947     swapVec[0] = make_pair( 1, 2 );
 948     swapVec[1] = make_pair( 4, 5 );
 949     break;
 950   case MED_HEXA8:
 951     swapVec.resize(2);
 952     swapVec[0] = make_pair( 1, 3 );
 953     swapVec[1] = make_pair( 5, 7 );
 954     break;
 955   case MED_TETRA10:
 956     swapVec.resize(3);
 957     swapVec[0] = make_pair( 1, 2 );
 958     swapVec[1] = make_pair( 4, 5 );
 959     swapVec[2] = make_pair( 8, 9 );
 960     break;
 961   case MED_PYRA13:
 962     swapVec.resize(4);
 963     swapVec[0] = make_pair( 1, 3 );
 964     swapVec[1] = make_pair( 5, 8 );
 965     swapVec[2] = make_pair( 6, 7 );
 966     swapVec[3] = make_pair( 10, 12 );
 967     break;
 968   case MED_PENTA15:
 969     swapVec.resize(4);
 970     swapVec[0] = make_pair( 1, 2 );
 971     swapVec[1] = make_pair( 4, 5 );
 972     swapVec[2] = make_pair( 6, 8 );
 973     swapVec[3] = make_pair( 9, 11 );
 974     break;
 975   case MED_HEXA20:
 976     swapVec.resize(7);
 977     swapVec[0] = make_pair( 1, 3 );
 978     swapVec[1] = make_pair( 5, 7 );
 979     swapVec[2] = make_pair( 8, 11 );
 980     swapVec[3] = make_pair( 9, 10 );
 981     swapVec[4] = make_pair( 12, 15 );
 982     swapVec[5] = make_pair( 13, 14 );
 983     swapVec[6] = make_pair( 17, 19 );
 984   default:;
 985   }
 986   END_OF("void getReverseVector()");
 987 }
 988
 989 //=======================================================================
 990 //function : orientElements
 991 //purpose  :
 992 //=======================================================================
 993
 994 static void orientElements( _intermediateMED& medi )
 995 {
 996   MESSAGE("orientElements()");
 997   set<_maille>::iterator elemIt = medi.maillage.begin();
 998
 999   if ( elemIt->sommets[0]->second.coord.size() == 2 ) { // space dimension
1000
1001     // --------------------------
1002     // Orient 2D faces clockwise
1003     // --------------------------
1004
1005     for ( ; elemIt != medi.maillage.end(); elemIt++ )
1006       if ( elemIt->dimension() == 2 )
1007       {
1008         // look for index of the most left node
1009         int iLeft = 0, iNode, nbNodes = elemIt->sommets.size();
1010         double minX = elemIt->sommets[0]->second.coord[0];
1011         for ( iNode = 1; iNode < nbNodes; ++iNode )
1012         {
1013           if ( minX > elemIt->sommets[ iNode ]->second.coord[ 0 ]) {
1014             minX = elemIt->sommets[ iNode ]->second.coord[ 0 ];
1015             iLeft = iNode;
1016           }
1017         }
1018         // indeces of the nodes neighboring the most left one
1019         int iPrev = ( iLeft - 1 < 0 ) ? nbNodes - 1 : iLeft - 1;
1020         int iNext = ( iLeft + 1 == nbNodes ) ? 0 : iLeft + 1;
1021         // find components of prev-left and left-next vectors
1022         double xP = elemIt->sommets[ iPrev ]->second.coord[ 0 ];
1023         double yP = elemIt->sommets[ iPrev ]->second.coord[ 1 ];
1024         double xN = elemIt->sommets[ iNext ]->second.coord[ 0 ];
1025         double yN = elemIt->sommets[ iNext ]->second.coord[ 1 ];
1026         double xL = elemIt->sommets[ iLeft ]->second.coord[ 0 ];
1027         double yL = elemIt->sommets[ iLeft ]->second.coord[ 1 ];
1028         double xPL = xL - xP, yPL = yL - yP; // components of prev-left vector
1029         double xLN = xN - xL, yLN = yN - yL; // components of left-next vector
1030         // normalise y of the vectors
1031         double modPL = sqrt ( xPL * xPL + yPL * yPL );
1032         double modLN = sqrt ( xLN * xLN + yLN * yLN );
1033         if ( modLN > DBL_MIN && modPL > DBL_MIN )
1034         {
1035           yPL /= modPL;
1036           yLN /= modLN;
1037           // summury direction of neighboring links must be positive
1038           bool clockwise = ( yPL + yLN > 0 );
1039           elemIt->reverse = ( !clockwise );
1040         }
1041       }
1042   }
1043   else {
1044
1045     int type = -100;
1046     vector< pair<int,int> > swapVec;
1047     for ( ; elemIt != medi.maillage.end(); elemIt++ ) {
1048       if ( elemIt->dimension() == 3 )
1049       {
1050         // ---------------------------------------------------
1051         // Orient volumes according to MED conventions:
1052         // normal of a bottom (first) face should be downward
1053         // ---------------------------------------------------
1054
1055         int nbBottomNodes = 0;
1056         switch ( elemIt->geometricType ) {
1057         case MED_TETRA4:
1058         case MED_TETRA10:
1059         case MED_PENTA6:
1060         case MED_PENTA15:
1061           nbBottomNodes = 3; break;
1062         case MED_PYRA5:
1063         case MED_PYRA13:
1064         case MED_HEXA8:
1065         case MED_HEXA20:
1066           nbBottomNodes = 4; break;
1067         default: continue;
1068         }
1069         // find a normal to the bottom face
1070         const _noeud* n[4] = {
1071           &elemIt->sommets[0]->second, // 3 bottom nodes
1072           &elemIt->sommets[1]->second,
1073           &elemIt->sommets[2]->second,
1074           &elemIt->sommets[nbBottomNodes]->second };// a top node
1075         double vec01 [3] = { // vector n[0]-n[1]
1076           n[1]->coord[0] - n[0]->coord[0],
1077           n[1]->coord[1] - n[0]->coord[1],
1078           n[1]->coord[2] - n[0]->coord[2], };
1079         double vec02 [3] = { // vector n[0]-n[2]
1080           n[2]->coord[0] - n[0]->coord[0],
1081           n[2]->coord[1] - n[0]->coord[1],
1082           n[2]->coord[2] - n[0]->coord[2] };
1083         double normal [3] = { // vec01 ^ vec02
1084           vec01[1] * vec02[2] - vec01[2] * vec02[1],
1085           vec01[2] * vec02[0] - vec01[0] * vec02[2],
1086           vec01[0] * vec02[1] - vec01[1] * vec02[0] };
1087         // check if the 102 angle is convex
1088         if ( nbBottomNodes > 3 ) {
1089           const _noeud* n3 = &elemIt->sommets[nbBottomNodes-1]->second;// last bottom node
1090           double vec03 [3] = { // vector n[0]-n3
1091             n3->coord[0] - n[0]->coord[0],
1092             n3->coord[1] - n[0]->coord[1],
1093             n3->coord[2] - n[0]->coord[2], };
1094           if ( fabs( normal[0]+normal[1]+normal[2] ) <= DBL_MIN ) { // vec01 || vec02
1095             normal[0] = vec01[1] * vec03[2] - vec01[2] * vec03[1]; // vec01 ^ vec03
1096             normal[1] = vec01[2] * vec03[0] - vec01[0] * vec03[2];
1097             normal[2] = vec01[0] * vec03[1] - vec01[1] * vec03[0];
1098           }
1099           else {
1100             double vec [3] = { // normal ^ vec01
1101               normal[1] * vec01[2] - normal[2] * vec01[1],
1102               normal[2] * vec01[0] - normal[0] * vec01[2],
1103               normal[0] * vec01[1] - normal[1] * vec01[0] };
1104             double dot2 = vec[0]*vec03[0] + vec[1]*vec03[1] + vec[2]*vec03[2]; // vec*vec03
1105             if ( dot2 < 0 ) { // concave -> reverse normal
1106               normal[0] *= -1;
1107               normal[1] *= -1;
1108               normal[2] *= -1;
1109             }
1110           }
1111         }
1112         // direction from top to bottom
1113         vector<double> tbDir(3);
1114         tbDir[0] = n[0]->coord[0] - n[3]->coord[0];
1115         tbDir[1] = n[0]->coord[1] - n[3]->coord[1];
1116         tbDir[2] = n[0]->coord[2] - n[3]->coord[2];
1117         // compare 2 directions: normal and top-bottom
1118         double dot = normal[0]*tbDir[0] + normal[1]*tbDir[1] + normal[2]*tbDir[2];
1119         bool reverse = ( dot < 0. );
1120         if ( reverse ) {
1121           if ( elemIt->geometricType != type ) {
1122             type = elemIt->geometricType;
1123             getReverseVector( type, swapVec );
1124 //             INFOS("vec01: " <<vec01[0] << " " <<vec01[1] << " " << vec01[2]);
1125 //             INFOS("vec02: " <<vec02[0] << " " <<vec02[1] << " " << vec02[2]);
1126 //             INFOS("normal: " <<normal[0] << " " <<normal[1] << " " << normal[2]);
1127 //             INFOS("tb: " << tbDir[0] << " " <<tbDir[1] << " " << tbDir[2]);
1128 //             INFOS( *elemIt );
1129 //             for ( vector< _maille::iter >::const_iterator si = elemIt->sommets.begin();
1130 //                  si != elemIt->sommets.end(); si++ )
1131 //               INFOS( (*si)->second );
1132           }
1133           _maille* ma = (_maille*) & (*elemIt);
1134           for ( int i = 0; i < swapVec.size(); ++i ) {
1135             _maille::iter tmp = ma->sommets[ swapVec[i].first ];
1136             ma->sommets[ swapVec[i].first ] = ma->sommets[ swapVec[i].second ];
1137             ma->sommets[ swapVec[i].second ] = tmp;
1138           }
1139         }
1140       } // dimension() == 3
1141     } // loop on maillage
1142
1143     // --------------------------------------
1144     // orient equally all connected 3D faces
1145     // --------------------------------------
1146
1147     // fill map of links and their faces
1148     set<const _maille*> faces;
1149     map<const _maille*, _groupe*> fgm;
1150     map<_link, list<const _maille*> > linkFacesMap;
1151     map<_link, list<const _maille*> >::iterator lfIt, lfIt2;
1152
1153     medi.treatGroupes(); // erase groupes that wont be converted
1154     for (unsigned int i=0; i!=medi.groupes.size(); ++i)
1155     {
1156       _groupe& grp = medi.groupes[i];
1157       _groupe::mailleIter maIt=grp.mailles.begin();
1158       if ( maIt==grp.mailles.end() || (*maIt)->dimension() != 2 )
1159         continue;
1160       for(; maIt!=grp.mailles.end(); ++maIt) {
1161         if ( faces.insert( &(**maIt )).second ) {
1162           for ( int j = 0; j < (*maIt)->sommets.size(); ++j )
1163             linkFacesMap[ (*maIt)->link( j ) ].push_back( &(**maIt) );
1164           fgm.insert( make_pair( &(**maIt), &grp ));
1165         }
1166       }
1167     }
1168     // dump linkFacesMap
1169 //     for ( lfIt = linkFacesMap.begin(); lfIt!=linkFacesMap.end(); lfIt++) {
1170 //       cout<< "LINK: " << lfIt->first.first << "-" << lfIt->first.second << endl;
1171 //       list<const _maille*> & fList = lfIt->second;
1172 //       list<const _maille*>::iterator fIt = fList.begin();
1173 //       for ( ; fIt != fList.end(); fIt++ )
1174 //         cout << "\t" << **fIt << fgm[*fIt]->nom << endl;
1175 //     }
1176
1177     // Each oriented link must appear in one face only, else a face is reversed.
1178
1179     queue<const _maille*> faceQueue; // the queue contains well oriented faces
1180     // whose neighbors orientation is to be checked
1181
1182     bool manifold = true;
1183     while ( !linkFacesMap.empty() )
1184     {
1185       if ( faceQueue.empty() ) {
1186         ASSERT( !linkFacesMap.begin()->second.empty() );
1187         faceQueue.push( linkFacesMap.begin()->second.front() );
1188       }
1189       while ( !faceQueue.empty() )
1190       {
1191         const _maille* face = faceQueue.front();
1192         faceQueue.pop();
1193
1194         // loop on links of <face>
1195         for ( int i = 0; i < face->sommets.size(); ++i ) {
1196           _link link = face->link( i );
1197           // find the neighbor faces
1198           lfIt = linkFacesMap.find( link );
1199           int nbFaceByLink = 0;
1200           list< const _maille* > ml;
1201           if ( lfIt != linkFacesMap.end() )
1202           {
1203             list<const _maille*> & fList = lfIt->second;
1204             list<const _maille*>::iterator fIt = fList.begin();
1205             ASSERT( fIt != fList.end() );
1206             for ( ; fIt != fList.end(); fIt++, nbFaceByLink++ ) {
1207               ml.push_back( *fIt );
1208               if ( *fIt != face ) // wrongly oriented neighbor face
1209               {
1210                 const _maille* badFace = *fIt;
1211                 // reverse and remove badFace from linkFacesMap
1212                 for ( int j = 0; j < badFace->sommets.size(); ++j ) {
1213                   _link badlink = badFace->link( j );
1214                   if ( badlink == link ) continue;
1215                   lfIt2 = linkFacesMap.find( badlink );
1216                   if ( lfIt2 != linkFacesMap.end() ) {
1217                     list<const _maille*> & ff = lfIt2->second;
1218                     ff.erase( find( ff.begin(), ff.end(), badFace ));
1219                         if ( ff.empty() )
1220                           linkFacesMap.erase( lfIt2 );
1221                   }
1222                 }
1223                 badFace->reverse = true; // reverse
1224                 //INFOS( "REVERSE " << *badFace );
1225                 faceQueue.push( badFace );
1226               }
1227             }
1228             linkFacesMap.erase( lfIt );
1229           }
1230           // add good neighbors to the queue
1231           _link revLink( link.second, link.first );
1232           lfIt = linkFacesMap.find( revLink );
1233           if ( lfIt != linkFacesMap.end() )
1234           {
1235             list<const _maille*> & fList = lfIt->second;
1236             list<const _maille*>::iterator fIt = fList.begin();
1237             for ( ; fIt != fList.end(); fIt++, nbFaceByLink++ ) {
1238               ml.push_back( *fIt );
1239               if ( *fIt != face )
1240                 faceQueue.push( *fIt );
1241             }
1242             linkFacesMap.erase( lfIt );
1243           }
1244           if ( nbFaceByLink > 2 ) {
1245             if ( manifold ) {
1246               list<const _maille*>::iterator i = ml.begin();
1247               INFOS(nbFaceByLink << " faces by 1 link:");
1248               for( ; i!= ml.end(); i++ ) {
1249                 INFOS("in object " << fgm[ *i ]->nom);
1250                 INFOS( **i );
1251               }
1252             }
1253             manifold = false;
1254           }
1255         } // loop on links of the being checked face
1256       } // loop on the face queue
1257     } // while ( !linkFacesMap.empty() )
1258
1259     if ( !manifold )
1260       INFOS(" -> Non manifold mesh, faces orientation may be incorrect");
1261
1262   } // space dimension == 3
1263 }
1264
1265 //=======================================================================
1266 //function : fillMesh
1267 //purpose  : load data from medi to mesh
1268 //=======================================================================
1269
1270 void GIBI_MESH_RDONLY_DRIVER::fillMesh(_intermediateMED* _ptrMedi)
1271 {
1272   const char * LOC = "GIBI_MESH_RDONLY_DRIVER::fillMesh(_intermediateMED* _ptrMedi) : " ;
1273   BEGIN_OF(LOC);
1274
1275   _ptrMesh->_name = _meshName;
1276
1277   if (_ptrMedi)
1278   {
1279     if (_ptrMedi->maillage.size()==0 ||
1280         _ptrMedi->groupes.size()==0 ||
1281         _ptrMedi->points.size()==0) {
1282       INFOS(" Error while reading file: the data read are not completed " ) ;
1283       return;
1284     }
1285     // fix element orientation
1286     orientElements( *_ptrMedi );
1287
1288     _ptrMesh->_spaceDimension = _ptrMedi->points.begin()->second.coord.size();
1289     _ptrMesh->_meshDimension = _ptrMedi->maillage.rbegin()->dimension();
1290     _ptrMesh->_numberOfNodes = _ptrMedi->points.size();
1291     _ptrMesh->_isAGrid = 0;
1292     _ptrMesh->_coordinate = _ptrMedi->getCoordinate();
1293
1294     //Construction des groupes
1295     _ptrMedi->getGroups(_ptrMesh->_groupCell,
1296                         _ptrMesh->_groupFace,
1297                         _ptrMesh->_groupEdge,
1298                         _ptrMesh->_groupNode, _ptrMesh);
1299
1300     _ptrMesh->_connectivity = _ptrMedi->getConnectivity();
1301
1302     // calcul de la connectivite d-1 complete, avec renumerotation des groupes
1303     //if (_ptrMesh->_spaceDimension==3)
1304     //    _ptrMesh->_connectivity->updateGroup(_ptrMesh->_groupFace) ;
1305     //else if (_ptrMesh->_spaceDimension==2)
1306     //    _ptrMesh->_connectivity->updateGroup(_ptrMesh->_groupEdge) ;
1307
1308     // Creation des familles à partir des groupes
1309     // NC : Cet appel pourra être différé quand la gestion de la cohérence famille/groupes sera assurée
1310     _ptrMesh->createFamilies();
1311     // TAKE CARE OF ELEMENTS ORDER IN GROUPS AFTER THEIR SPLITING INTO FAMILIES !!!!
1312     // _ptrMesh->createFamilies() breaks the order
1313 //     _ptrMedi->getFamilies(_ptrMesh->_familyCell,
1314 //                           _ptrMesh->_familyFace,
1315 //                           _ptrMesh->_familyEdge,
1316 //                           _ptrMesh->_familyNode, _ptrMesh);
1317
1318     // add attributes to families
1319     set<string> famNames;
1320     for (medEntityMesh entity=MED_CELL; entity<MED_ALL_ENTITIES; ++entity)
1321     {
1322       int i, nb = _ptrMesh->getNumberOfFamilies(entity);
1323       for ( i = 1; i <= nb; ++i ) {
1324         FAMILY* f = const_cast<FAMILY*>( _ptrMesh->getFamily( entity, i ));
1325         f->setNumberOfAttributes( 1 );
1326         int* attIDs = new int[1];
1327         attIDs[0] = 1;
1328         f->setAttributesIdentifiers( attIDs );
1329         int* attVals = new int[1];
1330         attVals[0] = 1;
1331         f->setAttributesValues( attVals );
1332         string* attDescr = new string[1];
1333         attDescr[0] = "med_family";
1334         f->setAttributesDescriptions( attDescr );
1335         // limit a name length
1336         if ( f->getName().length() > 31 ) {
1337           ostringstream name;
1338           name << "FAM" << f->getIdentifier();
1339           f->setName( name.str());
1340         }
1341         // check if family is on the whole mesh entity
1342         if (_ptrMesh->getNumberOfElements( entity, MED_ALL_ELEMENTS ) ==
1343             f->getNumberOfElements( MED_ALL_ELEMENTS ))
1344           f->setAll( true );
1345       }
1346       // setAll() for groups
1347       nb = _ptrMesh->getNumberOfGroups(entity);
1348       for ( i = 1; i <= nb; ++i ) {
1349         GROUP * g = const_cast<GROUP*>( _ptrMesh->getGroup( entity, i ));
1350         if (_ptrMesh->getNumberOfElements( entity, MED_ALL_ELEMENTS ) ==
1351             g->getNumberOfElements( MED_ALL_ELEMENTS ))
1352           g->setAll( true );
1353       }
1354     }
1355   }
1356   END_OF(LOC);
1357 }
1358
1359 void GIBI_MESH_RDONLY_DRIVER::write( void ) const
1360   throw (MEDEXCEPTION)
1361 {
1362   throw MEDEXCEPTION("GIBI_MESH_RDONLY_DRIVER::write : Can't write with a RDONLY driver !");
1363 }
1364
1365
1366 /*--------------------- WRONLY PART -------------------------------*/
1367
1368 GIBI_MESH_WRONLY_DRIVER::GIBI_MESH_WRONLY_DRIVER():GIBI_MESH_DRIVER()
1369 {
1370 }
1371 GIBI_MESH_WRONLY_DRIVER::GIBI_MESH_WRONLY_DRIVER(const string & fileName,
1372                                                  MESH * ptrMesh):
1373   GIBI_MESH_DRIVER(fileName,ptrMesh,MED_WRONLY)
1374 {
1375   MESSAGE("GIBI_MESH_WRONLY_DRIVER::GIBI_MESH_WRONLY_DRIVER(const string & fileName, MESH * ptrMesh) has been created");
1376 }
1377 GIBI_MESH_WRONLY_DRIVER::GIBI_MESH_WRONLY_DRIVER(const GIBI_MESH_WRONLY_DRIVER & driver):
1378   GIBI_MESH_DRIVER(driver)
1379 {
1380 }
1381 GIBI_MESH_WRONLY_DRIVER::~GIBI_MESH_WRONLY_DRIVER()
1382 {
1383   //MESSAGE("GIBI_MESH_WRONLY_DRIVER::GIBI_MESH_WRONLY_DRIVER(const string & fileName, MESH * ptrMesh) has been destroyed");
1384 }
1385 GENDRIVER * GIBI_MESH_WRONLY_DRIVER::copy(void) const
1386 {
1387   return new GIBI_MESH_WRONLY_DRIVER(*this);
1388 }
1389 void GIBI_MESH_WRONLY_DRIVER::read (void)
1390   throw (MEDEXCEPTION)
1391 {
1392   throw MEDEXCEPTION("GIBI_MESH_WRONLY_DRIVER::read : Can't read with a WRONLY driver !");
1393 }
1394
1395 //=======================================================================
1396 //function : open
1397 //purpose  :
1398 //=======================================================================
1399
1400 void GIBI_MESH_WRONLY_DRIVER::open()
1401   //     throw (MEDEXCEPTION)
1402 {
1403   const char * LOC = "GIBI_MESH_DRIVER::open()" ;
1404   BEGIN_OF(LOC);
1405
1406   MED_EN::med_mode_acces aMode = getAccessMode();
1407   switch (aMode) {
1408   case MED_EN::MED_REMP:
1409   case MED_EN::MED_ECRI: // should never append !!
1410     _gibi.open(_fileName.c_str(), ios::out);
1411     break;
1412   default:
1413     throw MEDEXCEPTION(LOCALIZED(STRING(LOC) << "Bad file mode access ! " << aMode));
1414   }
1415   if (_gibi &&
1416 #ifdef WNT
1417       _gibi.is_open()
1418 #else
1419       _gibi.rdbuf()->is_open()
1420 #endif
1421       )
1422     _status = MED_OPENED;
1423   else
1424   {
1425     _status = MED_CLOSED;
1426     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" Could not open file "<<_fileName));
1427   }
1428   END_OF(LOC);
1429 }
1430
1431 //=======================================================================
1432 //function : close
1433 //purpose  :
1434 //=======================================================================
1435
1436 void GIBI_MESH_WRONLY_DRIVER::close()
1437   //  throw (MEDEXCEPTION)
1438 {
1439     const char * LOC = "GIBI_MESH_DRIVER::close() " ;
1440     BEGIN_OF(LOC);
1441     if ( _status == MED_OPENED)
1442     {
1443         _gibi.close();
1444         _status = MED_CLOSED;
1445     }
1446     END_OF(LOC);
1447 }
1448
1449 //=======================================================================
1450 //function : write
1451 //purpose  :
1452 //=======================================================================
1453
1454 void GIBI_MESH_WRONLY_DRIVER::write(void) const
1455   throw (MEDEXCEPTION)
1456 {
1457   const char * LOC = "void GIBI_MESH_WRONLY_DRIVER::write(void) const : ";
1458   BEGIN_OF(LOC);
1459
1460   // we are going to modify the _gibi field
1461   GIBI_MESH_WRONLY_DRIVER * me = const_cast<GIBI_MESH_WRONLY_DRIVER *>(this);
1462 //  try {
1463     me->writeSupportsAndMesh();
1464     me->writeLastRecord();
1465 //   }
1466 //   catch (MEDEXCEPTION &ex)
1467 //   {
1468 //     INFOS( ex.what() );
1469 //   }
1470
1471   END_OF(LOC);
1472 }
1473
1474 //=======================================================================
1475 //function : getName
1476 //purpose  : return cleaned up support name
1477 //=======================================================================
1478
1479 static string cleanName( const string& theName )
1480 {
1481   string name = theName;
1482   if ( !name.empty() ) {
1483     // find a name string end
1484     int i, len = name.length();
1485     for ( i = 0; i < len; ++i ) {
1486       if ( name[i] == 0 )
1487         break;
1488     }
1489     // cut off trailing white spaces
1490     while ( i > 0 && name[i-1] == ' ' )
1491       i--;
1492     if ( i != len ) {
1493       name = name.substr( 0, i );
1494       len = i;
1495     }
1496   }
1497   return name;
1498 }
1499
1500 //=======================================================================
1501 //function : addSupport
1502 //purpose  :
1503 //=======================================================================
1504
1505 bool GIBI_MESH_WRONLY_DRIVER::addSupport( const SUPPORT * support )
1506 {
1507   if ( !support )
1508     return false;
1509   map<const SUPPORT*,supportData>::iterator su = _supports.find( support );
1510   if ( su != _supports.end() )
1511     return ( su->second.getNumberOfTypes() > 0 );
1512
1513   if ( support->getMesh() != _ptrMesh )
1514     throw MEDEXCEPTION(LOCALIZED(STRING("cant write support of other mesh" )));
1515
1516   // get sub-supports and define a support type name
1517   string supType;
1518   list<const SUPPORT*> sList;
1519   const GROUP* group = dynamic_cast< const GROUP* >(support);
1520   if ( group )
1521   {
1522     if ( group->getNumberOfTypes() > 0 || group->isOnAllElements() )
1523       sList.push_back( group );
1524     else {
1525       int iFam, nbFam = group->getNumberOfFamilies();
1526       for ( iFam = 1; iFam <= nbFam; ++iFam )
1527         sList.push_back( group->getFamily( iFam ));
1528     }
1529     supType = "group";
1530   }
1531   else
1532   {
1533     sList.push_back( support );
1534     supType = dynamic_cast< const FAMILY* >(support) ? "family" : "support";
1535   }
1536
1537   supportData & data = _supports[ support ];
1538   data._cleanName = cleanName( support->getName() );
1539
1540   // check if it is a writtable support, i.e.
1541   // nodal connectivity for a support entity exists
1542   medEntityMesh entity = support->getEntity();
1543   if ( entity != MED_NODE && !_ptrMesh->existConnectivity( MED_NODAL, entity )) {
1544     INFOS("Do not save " << supType << " of entity " << entity
1545           << " named <" << data._cleanName << "> nodal connectivity not defined");
1546     return false;
1547   }
1548
1549   // fill supportData
1550   list<const SUPPORT*>::iterator sIt = sList.begin();
1551   for ( ; sIt != sList.end(); sIt++ )
1552   {
1553     bool onAll = (*sIt)->isOnAllElements();
1554     int nbTypes = 0;
1555     if ( !onAll )
1556       nbTypes = (*sIt)->getNumberOfTypes();
1557     else
1558       nbTypes = _ptrMesh->getNumberOfTypes( entity );
1559     if ( nbTypes == 0 )
1560       continue;
1561     const medGeometryElement* types = 0;
1562     if ( !onAll )
1563       types = (*sIt)->getTypes();
1564     else if ( entity != MED_NODE )
1565       types = _ptrMesh->getTypes( entity );
1566     for ( int iType = 0; iType < nbTypes; ++iType )
1567     {
1568       medGeometryElement geomType = types ? types[ iType ] : MED_ALL_ELEMENTS;
1569       const int * ptrElemIDs = 0;
1570       int elemID1 = 0, nbElems = 0;
1571       if ( onAll ) {
1572         nbElems = _ptrMesh->getNumberOfElements( entity, geomType );
1573         elemID1 = (entity == MED_NODE) ? 1 : _ptrMesh->getGlobalNumberingIndex (entity)[ iType ];
1574       }
1575       else {
1576         nbElems = (*sIt)->getNumberOfElements( geomType );
1577         ptrElemIDs = (*sIt)->getNumber( geomType );
1578       }
1579       if ( geomType == 0 )
1580         geomType = MED_POINT1;
1581
1582       data.addTypeData( geomType, nbElems, ptrElemIDs, elemID1 );
1583     }
1584   }
1585
1586   if ( data.getNumberOfTypes() == 0 ) {
1587     INFOS("Do not save " << supType << " of entity " << entity
1588           << " named <" << data._cleanName << "> no geometric types");
1589     return false;
1590   }
1591
1592   return true;
1593 }
1594
1595 //=======================================================================
1596 //function : getSupportIndex
1597 //purpose  :
1598 //=======================================================================
1599
1600 int GIBI_MESH_WRONLY_DRIVER::getSubMeshIdAndSize(const SUPPORT *        support,
1601                                                  list<pair<int,int> > & idsAndSizes) const
1602 {
1603   idsAndSizes.clear();
1604   map<const SUPPORT*,supportData>::const_iterator su = _supports.find( support );
1605   if ( su == _supports.end() )
1606     return 0;
1607
1608   supportData * data = const_cast<supportData *>( & su->second );
1609   int id = data->_id;
1610   if ( data->getNumberObjects() > data->getNumberOfTypes() )
1611     id++;
1612   supportData::typeIterator tIt = data->_types.begin();
1613   for ( ; tIt != data->_types.end(); ++tIt )
1614   {
1615     int size = 0;
1616     list< typeData >& td = tIt->second;
1617     list< typeData >::iterator tdIt = td.begin();
1618     for ( ; tdIt != td.end(); ++tdIt )
1619       size += tdIt->_nbElems;
1620     idsAndSizes.push_back( make_pair( id++, size ));
1621   }
1622   return idsAndSizes.size();
1623 }
1624
1625 // ============================================================
1626 // the class writes endl to the file as soon as <limit> fields
1627 // have been written after the last endl
1628 // ============================================================
1629
1630 class TFieldCounter
1631 {
1632   fstream& _file;
1633   int _count, _limit;
1634  public:
1635   TFieldCounter(fstream& f, int limit=0): _file(f), _limit(limit) { init(); }
1636   void init(int limit=0) // init, is done by stop() as well
1637   { if (limit) _limit = limit; _count = 0; }
1638   void operator++(int) // next
1639   { if ( ++_count == _limit ) { _file << endl; init(); }}
1640   void stop() // init() and write endl if there was no endl after the last written field
1641   { if ( _count ) _file << endl; init(); }
1642 };
1643
1644 //=======================================================================
1645 //function : writeElements
1646 //purpose  : ptrElemIDs and elemID1 provide two alternative ways of giving
1647 //           elements to write.
1648 //           If elemSet != 0 then an element is
1649 //           ( addElemInSet ? <written and added to elemSet> : <ignored if id is in elemSet>)
1650 //=======================================================================
1651
1652 void GIBI_MESH_WRONLY_DRIVER::writeElements (medGeometryElement geomType,
1653                                              list< typeData >&  typeDataList,
1654                                              const int *        nodalConnect,
1655                                              const int *        nodalConnectIndex)
1656 {
1657   // ITYPEL : type de l'鬩ment 1=point, 2=segment ?eux noeuds...
1658   // NBSOUS : nombre de sous parties dans cet objet,
1659   //          une sous partie par type d'鬩ments le composant.
1660   // NBREF : nombre de sous r馩rences. Une r馩rence est par exemple le contour
1661   // NBNOEL : nombre de noeuds par 鬩ment
1662   // NBEL : nombre d'鬩ments
1663
1664   int castemType = GIBI_MESH_DRIVER::med2gibiGeom( geomType );
1665   char* zeroI8 = "       0"; // FORMAT(I8)
1666   int nbElemNodes = geomType % 100;
1667
1668   // count total nb of elements
1669   int nbElements = 0;
1670   list< typeData >::iterator td = typeDataList.begin();
1671   for ( ; td != typeDataList.end(); td++ )
1672     nbElements += td->_nbElems;
1673
1674   _gibi << setw(8) << castemType <<  // ITYPE
1675     zeroI8 <<                       // NBSOUS
1676       zeroI8 <<                     // NBREF
1677         setw(8) << nbElemNodes <<   // NBNOEL
1678           setw(8) << nbElements <<  // NBEL
1679             endl;
1680
1681   MESSAGE("writeElements(): geomType=" << geomType << " nbElements= " << nbElements)
1682
1683   // L 'enregistrement donnant le num? de la couleur des 鬩ments.
1684   // * 8000 FORMAT(10I8)
1685   TFieldCounter fcount( _gibi, 10 );
1686   int iElem = 0;
1687   for ( ; iElem < nbElements; ++iElem, fcount++ )
1688     _gibi << zeroI8;
1689   fcount.stop();
1690
1691   // Tableau des connectivit鳮 Description du premier 鬩ment puis du deuxi譥...
1692   // ATTENTION il ne s'agit pas de la num?tation vraie,
1693   // il faut la faire passer par le filtre du dernier tableau de la pile num? 32.
1694   //int nbSkipped = 0;
1695
1696   for ( td = typeDataList.begin(); td != typeDataList.end(); td++ )
1697   {
1698     for ( int i = 0; i < td->_nbElems; i++ )
1699     {
1700       iElem = td->_ptrElemIDs ? td->_ptrElemIDs[ i ] : td->_elemID1 + i;
1701       if ( geomType == MED_POINT1 )
1702       {
1703         _gibi << setw(8) << iElem;
1704         fcount++;
1705       }
1706       else
1707       {
1708         int nodeId = nodalConnectIndex[ iElem - 1 ] - 1;
1709         for ( int iNode = 0; iNode < nbElemNodes; ++iNode, fcount++ ) {
1710           _gibi << setw(8) << nodalConnect[ nodeId++ ];
1711         }
1712       }
1713     }
1714   }
1715
1716   fcount.stop();
1717 }
1718
1719 //=======================================================================
1720 //function : addName
1721 //purpose  : make name uppercase and shorter than 9, add it to nameNbMap,
1722 //           raise if not unique
1723 //=======================================================================
1724
1725 void GIBI_MESH_WRONLY_DRIVER::addName(map<string,int>& nameMap,
1726                                       string&          theName,
1727                                       int              index,
1728                                       string           prefix)
1729 {
1730   string name = cleanName( theName );
1731   if ( !name.empty() ) {
1732     int len = name.length();
1733     if ( len > 8 )
1734       throw MEDEXCEPTION(STRING("Can't write too long name: ") << name );
1735
1736     for ( int i = 0; i < len; ++i )
1737       name[i] = toupper( name[i] );
1738     if ( ! nameMap.insert( make_pair( name, index )).second )
1739       throw MEDEXCEPTION(STRING("Can't write not unique name: ") << name );
1740 //     bool ok = ( len <= 8 && len > 0 );
1741 //     if ( ok ) {
1742 //       for ( int i = 0; i < len; ++i )
1743 //         name[i] = toupper( name[i] );
1744 //       ok = nameMap.insert( make_pair( name, index )).second;
1745 //     }
1746 //     if ( !ok ) {
1747 //       char *str=new char[ prefix.size() + 13 ];
1748 //       int j = 1;
1749 //       do {
1750 //         sprintf( str, "%s_%d", prefix.c_str(), nameMap.size()+j );
1751 //         ok = nameMap.insert( make_pair( str, index )).second;
1752 //         j++;
1753 //       } while ( !ok );
1754 //       INFOS( "Save <" << name << "> as <" << str << ">");
1755 //       delete [] str;
1756 //     }
1757   }
1758 }
1759
1760 //=======================================================================
1761 //function : writeNames
1762 //purpose  :
1763 //=======================================================================
1764
1765 void GIBI_MESH_WRONLY_DRIVER::writeNames( map<string,int>& nameNbMap )
1766 {
1767   // La pile num? 1 est celle des objets de type maillage.
1768   // La ligne suivante donne le nom des objets maillages sauv鳮
1769   // * 8001       FORMAT(8(1X,A8))
1770   if ( !nameNbMap.empty() )
1771   {
1772     TFieldCounter fcount( _gibi, 8 );
1773     _gibi << left;
1774     map<string,int>::iterator nameNbIt = nameNbMap.begin();
1775     for ( ; nameNbIt != nameNbMap.end(); nameNbIt++, fcount++ ) {
1776       _gibi << " " << setw(8) << nameNbIt->first;
1777     }
1778     fcount.stop();
1779     _gibi << right;
1780     // La ligne suivante donne les num?s d'ordre, dans la pile,
1781     // des objets nomm?cit?pr飩demment.
1782     // *  8000 FORMAT(10I8)
1783     nameNbIt = nameNbMap.begin();
1784     for ( fcount.init(10); nameNbIt != nameNbMap.end(); nameNbIt++, fcount++ )
1785       _gibi << setw(8) << nameNbIt->second;
1786     fcount.stop();
1787   }
1788 }
1789
1790 //=======================================================================
1791 //function : writeSupportsAndMesh
1792 //purpose  :
1793 //=======================================================================
1794
1795 void GIBI_MESH_WRONLY_DRIVER::writeSupportsAndMesh()
1796 {
1797   const char * LOC = "void GIBI_MESH_WRONLY_DRIVER::writeSupportsAndMesh() ";
1798   BEGIN_OF(LOC);
1799
1800   if (_status!=MED_OPENED)
1801     throw MEDEXCEPTION(LOCALIZED(STRING(LOC) << "file " << _fileName<<  " is not opened." ));
1802   if (!_ptrMesh)
1803     throw MEDEXCEPTION(LOCALIZED(STRING(LOC) << "can't write a NULL mesh" ));
1804   if (!_ptrMesh->getConnectivityptr())
1805     throw MEDEXCEPTION(LOCALIZED(STRING(LOC) << "can't write a mesh with NULL connectivity" ));
1806
1807   // fill _supports with families and groups
1808   medEntityMesh entity;
1809   for (entity=MED_CELL; entity<MED_ALL_ENTITIES; ++entity)
1810   {
1811     int i, nb = _ptrMesh->getNumberOfGroups(entity);
1812     for ( i = 1; i <= nb; ++i )
1813       addSupport( _ptrMesh->getGroup( entity, i ));
1814 //     nb = _ptrMesh->getNumberOfFamilies(entity);
1815 //     for ( i = 1; i <= nb; ++i )
1816 //       addSupport( _ptrMesh->getFamily( entity, i ));
1817   }
1818
1819   // --------------------------------------------------------------------
1820   // Count total nb of objects: an object per an element type in support
1821   // plus an object per an element type not used in _supports.
1822   // Collect object names
1823   // --------------------------------------------------------------------
1824
1825   vector<int> nbSuppElemsByType(MED_HEXA20,0);
1826   map<string,int> nameNbMap;
1827   map<const SUPPORT*,supportData>::iterator supIt = _supports.begin();
1828   int i, nb_objects = 0;
1829   for ( ; supIt != _supports.end(); supIt++ )
1830   {
1831     supportData & data = supIt->second;
1832     int nbSupObj = data.getNumberObjects();
1833     if ( nbSupObj == 0 )
1834       continue;
1835     data._id = nb_objects + 1;
1836     nb_objects += nbSupObj;
1837
1838     addName( nameNbMap, data._cleanName, data._id, "C" );
1839     MESSAGE( "obj " << data._id << " " << data._cleanName);
1840
1841     // count elements: take into account supports on all elements and families only
1842     const SUPPORT* support = supIt->first;
1843     if ( support->isOnAllElements() || dynamic_cast< const FAMILY* >( support ))
1844     {
1845       supportData::typeIterator tIt = data._types.begin();
1846       for ( ; tIt != data._types.end(); ++tIt )
1847         if ( support->isOnAllElements() )
1848         {
1849           nbSuppElemsByType[ tIt->first ] = INT_MAX / 100;
1850         }
1851         else
1852         {
1853           list< typeData >& td = tIt->second;
1854           list< typeData >::iterator tdIt = td.begin();
1855           for ( ; tdIt != td.end(); ++tdIt )
1856             nbSuppElemsByType[ tIt->first] += tdIt->_nbElems;
1857         }
1858     }
1859   }
1860
1861   // count types of mesh elements that are not all in _supports
1862   int iType, nbTypes;
1863   entity = _ptrMesh->getConnectivityptr()->getEntity();
1864   for ( ; entity < MED_NODE; entity++ )
1865   {
1866     nbTypes = _ptrMesh->getNumberOfTypes( entity );
1867     if ( nbTypes == 0 || !_ptrMesh->existConnectivity( MED_NODAL, entity ))
1868       continue;
1869     const medGeometryElement* types = _ptrMesh->getTypes( entity );
1870     for ( iType = 0; iType < nbTypes; ++iType )
1871     {
1872       int nbElemInSups = nbSuppElemsByType[ types[ iType ]];
1873       int nbElemInMesh = _ptrMesh->getNumberOfElements(entity, types[ iType ]);
1874       if ( nbElemInSups < nbElemInMesh ) {
1875         nb_objects++;
1876         nbSuppElemsByType[ types[ iType ]] = -1; // to keep written elements of _supports
1877       }
1878     }
1879   }
1880
1881   // ------------
1882   // Write file
1883   // ------------
1884
1885   // Premier paquet dont le nombre de lignes ne varie pas.
1886   // On y trouve des indications g鮩rales.
1887   const int dim = _ptrMesh->getSpaceDimension();
1888   _gibi << " ENREGISTREMENT DE TYPE   4" << endl;
1889   _gibi << " NIVEAU  15 NIVEAU ERREUR   0 DIMENSION   " << dim <<endl;
1890   _gibi << " DENSITE  .00000E+00" << endl;
1891   _gibi << " ENREGISTREMENT DE TYPE   7" << endl;
1892   _gibi << " NOMBRE INFO CASTEM2000   8" <<endl;
1893   _gibi << " IFOUR  -1 NIFOUR   0 IFOMOD  -1 IECHO   1 IIMPI   0 IOSPI   0 ISOTYP   1" << endl;
1894   _gibi << " NSDPGE     0" << endl;
1895
1896   // Deuxi譥 paquet qui d馩nit toutes les piles
1897   // (une pile par type d'objet et certaines piles en plus).
1898   // Un enregistrement de type 2 pr鶩ent de l'飲iture d'une nouvelle pile,
1899   // celui de type 5 pr鶩ent de la fin.
1900   // * 800   FORMAT (' ENREGISTREMENT DE TYPE', I4)
1901   _gibi << " ENREGISTREMENT DE TYPE   2" << endl;
1902   // * 801     FORMAT(' PILE NUMERO',I4,'NBRE OBJETS NOMMES',I8,'NBRE OBJETS',I8)
1903   _gibi << " PILE NUMERO   1NBRE OBJETS NOMMES" << setw(8) << nameNbMap.size() <<
1904     "NBRE OBJETS" << setw(8) << nb_objects <<endl;
1905
1906   writeNames( nameNbMap );
1907
1908   // Passage ?a description des objets les uns apr賠les autres.
1909   // Le premier enregistrement de chaque objet est compos頤e 5 nombres repr鳥ntant :
1910   // ITYPEL : type de l'鬩ment 1=point, 2=segment ?eux noeuds...
1911   // NBSOUS : nombre de sous parties dans cet objet,
1912   //          une sous partie par type d'鬩ments le composant.
1913   // NBREF : nombre de sous r馩rences. Une r馩rence est par exemple le contour
1914   // NBNOEL : nombre de noeuds par 鬩ment
1915   // NBEL : nombre d'鬩ments
1916   // Si ITYPEL=0 alors NBSOUS diff鲥nt de z?. Dans ce cas on lira la liste des positions,
1917   // dans la pile des objets, des sous parties le composant.
1918   // Si NBSOUS=0, NBNOEL et NBEL sont diff鲥nts de z?, on trouve, au besoin,
1919   // la liste des r馩rences , les num?s des couleurs puis les connectivit鳮
1920
1921   TFieldCounter fcount( _gibi, 10 );
1922   char* zeroI8 = "       0"; // FORMAT(I8)
1923   for ( supIt = _supports.begin(); supIt != _supports.end(); supIt++ )
1924   {
1925     supportData & data = supIt->second;
1926     int nbSupObj = data.getNumberObjects();
1927     if ( nbSupObj == 0 )
1928       continue;
1929     MESSAGE("support " << data._id << "<" << data._cleanName << ">");
1930
1931     // write a compound object
1932     int nbTypes = data.getNumberOfTypes();
1933     if ( nbSupObj > nbTypes )
1934     {
1935       _gibi << zeroI8 << setw(8) << nbTypes << zeroI8 << zeroI8 << zeroI8 << endl;
1936       for ( int i_sub = 1; i_sub <= nbTypes; ++i_sub, fcount++ )
1937         _gibi << setw(8) << ( data._id + i_sub );
1938       fcount.stop();
1939     }
1940
1941     // write components
1942     entity = supIt->first->getEntity();
1943     const int * nodalConnect = 0, * nodalConnectIndex = 0;
1944     if ( entity != MED_NODE ) {
1945       nodalConnect = _ptrMesh->getConnectivity (MED_FULL_INTERLACE,
1946                                                 MED_NODAL,
1947                                                 entity,
1948                                                 MED_ALL_ELEMENTS);
1949       nodalConnectIndex = _ptrMesh->getConnectivityIndex (MED_NODAL,entity);
1950     }
1951     supportData::typeIterator tIt = data._types.begin();
1952     for ( ; tIt != data._types.end(); ++tIt )
1953     {
1954       writeElements (tIt->first,
1955                      tIt->second,
1956                      nodalConnect,
1957                      nodalConnectIndex);
1958     }
1959   }  // loop on _supports
1960
1961   // Write elements that are not in _supports
1962
1963   supportData data;
1964   entity = _ptrMesh->getConnectivityptr()->getEntity();
1965   for ( ; entity < MED_NODE; entity++ )
1966   {
1967     int nbTypes = _ptrMesh->getNumberOfTypes( entity );
1968     if ( nbTypes == 0 || !_ptrMesh->existConnectivity( MED_NODAL, entity ))
1969       continue;
1970     const medGeometryElement* types = _ptrMesh->getTypes( entity );
1971     const int * nbIndex = _ptrMesh->getGlobalNumberingIndex (entity);
1972     const int * nodalConnect = 0, * nodalConnectIndex = 0;
1973     nodalConnect = _ptrMesh->getConnectivity (MED_FULL_INTERLACE,
1974                                               MED_NODAL,
1975                                               entity,
1976                                               MED_ALL_ELEMENTS);
1977     nodalConnectIndex = _ptrMesh->getConnectivityIndex (MED_NODAL,entity);
1978
1979     for ( int iType = 1; iType <= nbTypes; ++iType )
1980     {
1981       int nbElements = nbIndex[ iType ] - nbIndex[ iType - 1 ];
1982       medGeometryElement geomType = types[ iType - 1 ];
1983       if ( nbSuppElemsByType[ geomType ] >= nbElements )
1984         continue; // all elements are written with _supports
1985
1986       int elemId1 = nbIndex[ iType - 1 ];
1987       data.addTypeData( geomType, nbElements, 0, elemId1 );
1988
1989       writeElements (geomType,
1990                      data._types[ geomType ],
1991                      nodalConnect,
1992                      nodalConnectIndex);
1993     }
1994   }
1995
1996   // D颵t de la pile 32 (celle des points)
1997
1998   int nbNodes = _ptrMesh->getNumberOfNodes();
1999   _gibi << " ENREGISTREMENT DE TYPE   2" << endl;
2000   _gibi << " PILE NUMERO  32NBRE OBJETS NOMMES       0" <<
2001     "NBRE OBJETS" << setw(8) << nbNodes << endl;
2002   // Liste des noms de points
2003   // * 8001       FORMAT(8(1X,A8))
2004   // No named nodes
2005   // suit le nombre de noeuds
2006   _gibi << setw(8) << nbNodes << endl;
2007   // Le tableau suivant donne le filtre pour avoir le vrai num? des noeuds
2008   // appartenant aux 鬩ments d飲its. Par exemple, si un 鬩ment, d飲it
2009   // dans la pile 1, fait r馩rence ?n num? de noeud 駡l ? il faut le
2010   // mettre 駡l ?2
2011   // * 8000 FORMAT(10I8)
2012   for ( i = 0; i < nbNodes; ++i, fcount++ )
2013     _gibi << setw(8) << i + 1;
2014   fcount.stop();
2015
2016   // D颵t de pile 33 (celle des configurations (coordonn?))
2017   _gibi << " ENREGISTREMENT DE TYPE   2" << endl;
2018   _gibi << " PILE NUMERO  33NBRE OBJETS NOMMES       0NBRE OBJETS       1" << endl;
2019   // Suit le nombre de points dont on donne les coordonn?
2020   int nbValues = nbNodes * ( dim + 1 );
2021   _gibi << setw(8) << nbValues << endl;
2022   // Les coordonn? sont donn? par noeuds. D'abord le premier puis le deuxi譥...
2023   // Pour chaque noeuds, on donne les 2 ou 3 coordonn? plus la densit頣ourante
2024   // au moment de sa cr顴ion.
2025   // * 8003   FORMAT(1P,3E22.14)
2026   _gibi.precision(14);
2027   _gibi.setf( ios_base::scientific, ios_base::floatfield );
2028   _gibi.setf( ios_base::uppercase );
2029   const double * coords = _ptrMesh->getCoordinates(MED_FULL_INTERLACE);
2030   int j = 0;
2031   for ( fcount.init(3),i = 0; i < nbNodes; ++i, j += dim )
2032   {
2033     for ( int iCoord = 0; iCoord < dim; ++iCoord, fcount++ )
2034       _gibi << setw(22) << coords[ j + iCoord ];
2035     _gibi << setw(22) << 0.0; // densite
2036     fcount++;
2037   }
2038   fcount.stop();
2039
2040   END_OF(LOC);
2041 }
2042
2043 //=======================================================================
2044 //function : writeLastRecord
2045 //purpose  :
2046 //=======================================================================
2047
2048 void GIBI_MESH_WRONLY_DRIVER::writeLastRecord()
2049 {
2050   _gibi << " ENREGISTREMENT DE TYPE   5" << endl;
2051   _gibi << "LABEL AUTOMATIQUE :   1" << endl;
2052 }
2053
2054 /*--------------------- RDWR PART -------------------------------*/
2055
2056 GIBI_MESH_RDWR_DRIVER::GIBI_MESH_RDWR_DRIVER():GIBI_MESH_DRIVER()
2057 {
2058 }
2059 GIBI_MESH_RDWR_DRIVER::GIBI_MESH_RDWR_DRIVER(const string & fileName,
2060                                              MESH * ptrMesh):
2061        GIBI_MESH_DRIVER(fileName,ptrMesh,MED_RDWR)
2062 {
2063   MESSAGE("GIBI_MESH_RDWR_DRIVER::GIBI_MESH_RDWR_DRIVER(const string & fileName, MESH * ptrMesh) has been created");
2064 }
2065 GIBI_MESH_RDWR_DRIVER::GIBI_MESH_RDWR_DRIVER(const GIBI_MESH_RDWR_DRIVER & driver):
2066   GIBI_MESH_RDONLY_DRIVER::GIBI_MESH_DRIVER(driver)
2067 {
2068   MESSAGE("GIBI_MESH_RDWR_DRIVER::GIBI_MESH_RDWR_DRIVER(driver) has been created");
2069 }
2070 GIBI_MESH_RDWR_DRIVER::~GIBI_MESH_RDWR_DRIVER() {
2071   MESSAGE("GIBI_MESH_RDWR_DRIVER::GIBI_MESH_RDWR_DRIVER(const string & fileName, MESH * ptrMesh) has been destroyed");
2072 }
2073 GENDRIVER * GIBI_MESH_RDWR_DRIVER::copy(void) const
2074 {
2075   BEGIN_OF( "GIBI_MESH_RDWR_DRIVER::copy()");
2076   GENDRIVER * driver = new GIBI_MESH_RDWR_DRIVER(*this);
2077   END_OF( "GIBI_MESH_RDWR_DRIVER::copy()");
2078   return driver;
2079 }
2080 void GIBI_MESH_RDWR_DRIVER::write(void) const
2081   throw (MEDEXCEPTION)
2082 {
2083   GIBI_MESH_RDWR_DRIVER * me = const_cast<GIBI_MESH_RDWR_DRIVER *>(this);
2084   me->GIBI_MESH_WRONLY_DRIVER::open();
2085   me->GIBI_MESH_WRONLY_DRIVER::write();
2086   me->GIBI_MESH_WRONLY_DRIVER::close();
2087 }
2088 void GIBI_MESH_RDWR_DRIVER::read (void)
2089   throw (MEDEXCEPTION)
2090 {
2091   BEGIN_OF( "GIBI_MESH_RDWR_DRIVER::read()");
2092   GIBI_MESH_RDONLY_DRIVER::open();
2093   GIBI_MESH_RDONLY_DRIVER::read();
2094   GIBI_MESH_RDONLY_DRIVER::close();
2095   END_OF( "GIBI_MESH_RDWR_DRIVER::read()");
2096 }
2097 void GIBI_MESH_RDWR_DRIVER::open()
2098   // throw (MEDEXCEPTION)
2099 {
2100 }
2101 void GIBI_MESH_RDWR_DRIVER::close()
2102   // throw (MEDEXCEPTION)
2103 {
2104 }
2105
2106 //============================== ====================================================
2107 //============================== FIELD Reading Driver ==============================
2108 //============================== ====================================================
2109
2110 GIBI_MED_RDONLY_DRIVER::GIBI_MED_RDONLY_DRIVER():GIBI_MESH_RDONLY_DRIVER()
2111 {
2112 }
2113 GIBI_MED_RDONLY_DRIVER::GIBI_MED_RDONLY_DRIVER(const string & fileName, MED * ptrMed):
2114        GIBI_MESH_RDONLY_DRIVER(fileName,NULL), _med( ptrMed )
2115 {
2116   MESSAGE("GIBI_MED_RDONLY_DRIVER(const string & fileName, MED * ptrMed) has been created");
2117   _fileName = fileName;
2118   _accessMode = MED_RDONLY;
2119 }
2120 GIBI_MED_RDONLY_DRIVER::GIBI_MED_RDONLY_DRIVER(const GIBI_MED_RDONLY_DRIVER & driver)
2121 {
2122 }
2123  GIBI_MED_RDONLY_DRIVER::~GIBI_MED_RDONLY_DRIVER()
2124 {
2125 }
2126 GENDRIVER * GIBI_MED_RDONLY_DRIVER::copy ( void ) const
2127 {
2128   return new GIBI_MED_RDONLY_DRIVER(*this);
2129 }
2130
2131 //=======================================================================
2132 //function : read
2133 //purpose  :
2134 //=======================================================================
2135
2136 void GIBI_MED_RDONLY_DRIVER::read ( void ) throw (MEDEXCEPTION)
2137 {
2138   const char * LOC = "GIBI_MED_RDONLY_DRIVER::read() : " ;
2139   BEGIN_OF(LOC);
2140
2141   if (_status!=MED_OPENED)
2142     throw MEDEXCEPTION(LOCALIZED(STRING(LOC) << "file " << _fileName<<" is not opened." ));
2143
2144   _ptrMesh = new MESH();
2145
2146   _intermediateMED medi;
2147   try {
2148     if ( !readFile( &medi, true ) )
2149       return;
2150
2151     // set name of field and its support if it is empty
2152     set<string> fnames;
2153     list< _fieldBase* >::iterator fIt = medi.fields.begin();
2154     for ( ; fIt != medi.fields.end(); fIt++ )
2155       fnames.insert( (*fIt)->name );
2156     int i = 0;
2157     for (fIt = medi.fields.begin(); fIt != medi.fields.end(); fIt++ ) {
2158       _fieldBase* f = *fIt;
2159       if ( f->name.empty() ) {
2160         do {
2161           ostringstream name;
2162           name << "F_" << ++i;
2163           f->name = name.str();
2164         } while ( !fnames.insert( f->name ).second );
2165       }
2166       _groupe& g = medi.groupes[ f->group_id ];
2167       if ( g.nom.empty() )
2168         g.nom = "GRP_" + f->name;
2169     }
2170     //MESSAGE(LOC <<  medi );
2171     fillMesh( &medi );
2172     MESSAGE(LOC << "GIBI_MED_RDONLY_DRIVER::read : RESULTATS STRUCTURE INTERMEDIAIRES : ");
2173     MESSAGE(LOC <<  medi );
2174
2175     list< FIELD_* > fields;
2176     medi.getFields( fields );
2177     MESSAGE( "nb fields: " << fields.size() );
2178
2179     if ( _ptrMesh->getName().empty() )
2180       _ptrMesh->setName( "MESH" );
2181
2182     _med->addMesh( _ptrMesh );
2183
2184     list< FIELD_* >::iterator it = fields.begin();
2185     for ( ; it != fields.end(); it++ )
2186       _med->addField( *it );
2187   }
2188   catch (MEDEXCEPTION &ex)
2189   {
2190     INFOS( ex.what() );
2191   }
2192
2193   END_OF(LOC);
2194 }
2195
2196 //============================== ====================================================
2197 //============================== FIELD Writting Driver ==============================
2198 //============================== ====================================================
2199
2200 GIBI_MED_WRONLY_DRIVER::GIBI_MED_WRONLY_DRIVER():GIBI_MESH_WRONLY_DRIVER()
2201 {
2202 }
2203 GIBI_MED_WRONLY_DRIVER::GIBI_MED_WRONLY_DRIVER(const string & fileName,
2204                                                    MED * ptrMed,
2205                                                    MESH * ptrMesh)
2206      :GIBI_MESH_WRONLY_DRIVER(fileName,ptrMesh), _med( ptrMed )
2207 {
2208   const char * LOC =
2209     "GIBI_MED_WRONLY_DRIVER(const string & fileName, MED * ptrMed, MESH * ptrMesh)" ;
2210   BEGIN_OF(LOC);
2211
2212   _fileName = fileName;
2213   _accessMode = MED_WRONLY;
2214   _ptrMesh = ptrMesh;
2215   if ( !_med || !_ptrMesh )
2216     throw MEDEXCEPTION(LOCALIZED(STRING(LOC) << " Bad params " << ptrMed << " " << ptrMesh ));
2217 }
2218 GIBI_MED_WRONLY_DRIVER::GIBI_MED_WRONLY_DRIVER(const GIBI_MED_WRONLY_DRIVER & driver)
2219 {
2220 }
2221 GIBI_MED_WRONLY_DRIVER::~GIBI_MED_WRONLY_DRIVER()
2222 {
2223 }
2224 GENDRIVER * GIBI_MED_WRONLY_DRIVER::copy ( void ) const
2225 {
2226   return new GIBI_MED_WRONLY_DRIVER(*this);
2227 }
2228
2229 //=======================================================================
2230 //function : writeDataSection
2231 //purpose  :
2232 //=======================================================================
2233
2234 template< class T >
2235   static void writeDataSection (fstream& file,
2236                                 FIELD_*  field,
2237                                 int      id1,
2238                                 int      id2)
2239 {
2240   FIELD<T>* f = dynamic_cast<FIELD<T>*>( field );
2241   if (!f) return;
2242   MEDARRAY<T>* array = f->getvalue();
2243   int ld = array->getLeadingValue();
2244   //SCRUTE( array->getLengthValue() );
2245   for ( int iComp = 0; iComp < ld; ++iComp )
2246   {
2247     file << setw(8) << 1          // nb scalar values by element
2248       << setw(8) << ( id2 - id1 ) // total nb of scalar values
2249         << setw(8) << 0
2250           << setw(8) << 0 << endl;
2251     // * 8003   FORMAT(1P,3E22.14)
2252     int id = id1;
2253     while ( id < id2 )
2254     {
2255       for ( int i = 0; id < id2 && i < 3; ++i )
2256         file << setw(22) << array->getIJ( id++, iComp + 1);
2257       file << endl;
2258     }
2259   }
2260 }
2261
2262 //=======================================================================
2263 //function : getSupports
2264 //purpose  :
2265 //=======================================================================
2266
2267 // int GIBI_MED_WRONLY_DRIVER::getSupports(const FIELD_*         field,
2268 //                                         list<const SUPPORT*>& supList) const
2269 // {
2270 //   int nbSup = 0;
2271 //   const SUPPORT* sup = field->getSupport();
2272 //   const GROUP* group = dynamic_cast<const GROUP*>( sup );
2273 //   if ( group ) {
2274 //     int nbFam = group->getNumberOfFamilies();
2275 //     for ( int iFam = 1; iFam <= nbFam; ++iFam ) {
2276 //       const FAMILY* fam = group->getFamily( iFam );
2277 //       if ( getSupportIndex( fam ) >= 0 ) {
2278 //         supList.push_back( fam );
2279 //         nbSup++;
2280 //       }
2281 //     }
2282 //   }
2283 //   else {
2284 //     supList.push_back( sup );
2285 //     nbSup++;
2286 //   }
2287 //   return nbSup;
2288 // }
2289
2290 //=======================================================================
2291 //function : write
2292 //purpose  :
2293 //=======================================================================
2294
2295 void GIBI_MED_WRONLY_DRIVER::write( void ) const throw (MEDEXCEPTION)
2296 {
2297   const char * LOC = "void GIBI_MED_WRONLY_DRIVER::write(void) const : ";
2298   BEGIN_OF(LOC);
2299
2300   // we are going to modify the _gibi field
2301   GIBI_MED_WRONLY_DRIVER * me = const_cast<GIBI_MED_WRONLY_DRIVER *>(this);
2302
2303   // get all fields on _ptrMesh and add their support to be written
2304   list<FIELD_*> fields;
2305   int iField, nbFileds = _med->getNumberOfFields();
2306   int nb_obj = 0;
2307   list<int> nb_sub_list;
2308   map<string,int> nameNbMap;
2309
2310   list<pair<int,int> >           subIdSizeList; // pair( <submesh id>, <submesh size> );
2311   list<pair<int,int> >::iterator idsize;
2312
2313   string *names=new string[ nbFileds ];
2314   _med->getFieldNames( names );
2315   for ( iField = 0; iField < nbFileds; ++iField )
2316   {
2317     int nb_sub = 0;
2318     deque<DT_IT_> dtit = _med->getFieldIteration( names[ iField ]);
2319     deque<DT_IT_>::iterator fIt = dtit.begin();
2320     for ( ; fIt != dtit.end(); fIt++ )
2321     {
2322       FIELD_ * f = _med->getField( names[ iField ], fIt->dt, fIt->it );
2323       if ( !dynamic_cast< FIELD<double >* >( f ))
2324       {
2325         MESSAGE("GIBI_MED_WRONLY_DRIVER::write( FIELD< int > ) not implemented");
2326         continue;
2327       }
2328       const SUPPORT * sup = f->getSupport();
2329       if ( me->addSupport( sup ) ) {
2330         fields.push_back( f );
2331         nb_sub += getSubMeshIdAndSize( sup, subIdSizeList );
2332       }
2333     }
2334     if ( nb_sub ) {
2335       addName( nameNbMap, names[ iField ], ++nb_obj, "F" );
2336       nb_sub_list.push_back( nb_sub );
2337     }
2338   }
2339
2340   // write mesh
2341
2342   //try {
2343     me->writeSupportsAndMesh();
2344 //   }
2345 //   catch (MEDEXCEPTION &ex)
2346 //   {
2347 //     INFOS( ex.what() );
2348 //     return;
2349 //   }
2350
2351   // write fields
2352
2353   if ( !fields.empty() ) {
2354
2355     fstream & gibi = me->_gibi;
2356
2357     TFieldCounter fcount( gibi, 10 );
2358
2359     gibi << " ENREGISTREMENT DE TYPE   2" << endl;
2360     gibi << " PILE NUMERO  39NBRE OBJETS NOMMES" << setw(8) << nameNbMap.size()
2361       << "NBRE OBJETS" << setw(8) << nb_obj << endl;
2362
2363     me->writeNames( nameNbMap );
2364
2365     list<FIELD_*>::iterator itF = fields.begin();
2366     list<int>::iterator itNbSub = nb_sub_list.begin();
2367     int nb_sub = 0, cur_nb_sub = 0;
2368     for ( ; itF != fields.end(); itF++ )
2369     {
2370       if ( cur_nb_sub == nb_sub && itNbSub != nb_sub_list.end() ) {
2371         // start the next field writting
2372         nb_sub = *(itNbSub++);
2373         gibi << setw(8) << nb_sub << "      -1       6      72" << endl;
2374         gibi << left;
2375         gibi << setw(72) << " Field" << endl;
2376         gibi << right;
2377         gibi << setw(72) << " " << endl;
2378
2379         // Sub Components section
2380         list<FIELD_*>::iterator itF2 = itF;
2381         vector<int> vals( 9, 0 );
2382         vals[8] = 2;
2383         fcount.init(10);
2384         cur_nb_sub = 0;
2385         while ( itF2 != fields.end() && cur_nb_sub < nb_sub )
2386         {
2387           FIELD_* f = *itF2++;
2388           vals[2] = f->getNumberOfComponents();
2389           getSubMeshIdAndSize( f->getSupport(), subIdSizeList );
2390           for ( idsize = subIdSizeList.begin(); idsize != subIdSizeList.end(); idsize++ )
2391           {
2392             ++cur_nb_sub;
2393             vals[0] = -idsize->first; // support id
2394             for ( int i = 0; i < vals.size(); ++i, fcount++ )
2395               gibi << setw(8) << vals[ i ];
2396           }
2397         }
2398         fcount.stop();
2399         cur_nb_sub = 0;
2400         // dummy strings
2401         int i_sub;
2402         for ( fcount.init(4), i_sub = 0; i_sub < nb_sub; ++i_sub, fcount++ )
2403           gibi << "                  ";
2404         fcount.stop();
2405         for ( fcount.init(8), i_sub = 0; i_sub < nb_sub; ++i_sub, fcount++ )
2406           gibi << "         ";
2407         fcount.stop();
2408       }
2409
2410       FIELD_* f = *itF;
2411       int id1 = 1;
2412       int iComp, nbComp = f->getNumberOfComponents();
2413       // loop on sub-components
2414       getSubMeshIdAndSize( f->getSupport(), subIdSizeList );
2415       for ( idsize = subIdSizeList.begin(); idsize != subIdSizeList.end(); idsize++ )
2416       {
2417         cur_nb_sub++;
2418         // component addresses
2419         for ( fcount.init(10), iComp = 0; iComp < nbComp; ++iComp, fcount++ )
2420           gibi << setw(8) << 777; // a good number
2421         fcount.stop();
2422         // component names
2423         gibi << left;
2424         for ( fcount.init(8), iComp = 0; iComp < nbComp; ++iComp, fcount++ )
2425           gibi << " "  << setw(8) << f->getComponentName( iComp + 1 );
2426         fcount.stop();
2427         // component types
2428         for ( fcount.init(4), iComp = 0; iComp < nbComp; ++iComp, fcount++ )
2429           gibi << " "  << setw(17) << "REAL*8";
2430         fcount.stop();
2431         gibi << right;
2432
2433         // Data section
2434         int id2 = id1 + idsize->second;
2435         writeDataSection<double>( gibi, f, id1, id2 );
2436         id1 = id2;
2437       }
2438     } // loop on fields
2439   }
2440   me->writeLastRecord();
2441   delete [] names;
2442   END_OF(LOC);
2443 }