src/MEDMEM/MEDMEM_GibiMeshDriver.cxx

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