src/MEDMEM/MEDMEM_GibiMeshDriver.cxx

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