src/MEDMEM/MEDMEM_EnsightFieldDriver.cxx

   1 //  Copyright (C) 2007-2008  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 //  Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 //  CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 //  This library is free software; you can redistribute it and/or
   7 //  modify it under the terms of the GNU Lesser General Public
   8 //  License as published by the Free Software Foundation; either
   9 //  version 2.1 of the License.
  10 //
  11 //  This library is distributed in the hope that it will be useful,
  12 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 //  Lesser General Public License for more details.
  15 //
  16 //  You should have received a copy of the GNU Lesser General Public
  17 //  License along with this library; if not, write to the Free Software
  18 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 //  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22 #include <fstream>
  23 #include <sstream>
  24 #include <iomanip>
  25
  26 #include "MEDMEM_Utilities.hxx"
  27
  28 #include "MEDMEM_Field.hxx"
  29 #include "MEDMEM_EnsightUtils.hxx"
  30
  31 #define TStrTool _ASCIIFileReader
  32
  33 using namespace MED_EN;
  34 using namespace MEDMEM_ENSIGHT;
  35
  36 namespace {
  37
  38   // ==============================================================================
  39   /*!
  40    * \brief Iterator on specified values of a field component
  41    */
  42   template <typename T> class _SelectedValueIterator
  43   {
  44     int        myDelta;
  45     const T*   myPtr;
  46     int        myIndexDelta;
  47     const int* myIndex;
  48   public:
  49     _SelectedValueIterator() // by default next() returns zero
  50       : myDelta( 0 ),      myPtr  (_ValueIterator<T>::zeroPtr()),
  51         myIndexDelta( 0 ), myIndex(_ValueIterator<int>::zeroPtr()) {}
  52
  53     _SelectedValueIterator(const T* values, int delta, const int* index)
  54       : myDelta( delta ), myPtr(values), myIndexDelta(1), myIndex(index-1) {}
  55
  56     const T & next() { myIndex += myIndexDelta; return myPtr[ (*myIndex) * myDelta ]; }
  57   };
  58
  59   // ==============================================================================
  60   /*!
  61    * \brief Container of values of any type to read/write from/to EnSight variable.
  62    *  It stores values relating to one element type within a part
  63    */
  64   struct _SubPartValues
  65   {
  66     _SubPart      mySubPart;
  67     const char *  myValues;
  68     medModeSwitch myInterlace;
  69     bool          myOwnValues; // to delete or not myValues
  70     string        myConstValue; // if myConstValue is set, other members are meaningless
  71
  72     // for writing values of a field-on-all-entities for a group
  73     const int *   myNumbers;
  74
  75     string        myUndefValue;
  76     set<int>      myUndefIndices;
  77     vector<int>   myPartialIndices;
  78
  79     _SubPartValues(const _SubPart& sub=_SubPart())
  80       : mySubPart(sub), myValues(0), myInterlace(MED_NO_INTERLACE), myOwnValues(true) {};
  81
  82     ~_SubPartValues() { clear(); }
  83     void clear() { if ( myOwnValues ) delete [] myValues; myValues = 0; }
  84
  85     // -------------------------------------------------------------------
  86     //!< copying removes data from the source object
  87     _SubPartValues(const _SubPartValues& other)
  88       :mySubPart(other.mySubPart), myValues(other.myValues),
  89        myInterlace(other.myInterlace), myOwnValues(true)
  90     { ((_SubPartValues&)other).myValues=0; }
  91
  92     // -------------------------------------------------------------------
  93     //!< Return iterator on values of a component counted from 1
  94     //   nbElements is a number of values of a component located one by one in
  95     //   not full interlace
  96     template <typename T> _ValueIterator<T>
  97     getValues( int component, int nbElements, int nbComponents ) const
  98     {
  99       const T* values = (T*) myValues;
 100       if ( myInterlace == MED_FULL_INTERLACE ) {
 101         values += (component - 1);
 102         return _ValueIterator<T>( values, nbComponents );
 103       }
 104       else {
 105         values += (component - 1) * nbElements;
 106         return _ValueIterator<T>( values, 1 );
 107       }
 108     }
 109     // -------------------------------------------------------------------
 110     //!< Return iterator on specified values of a component counted from 1.
 111     //   nbElements is a number of values of a component located one by one in
 112     //   not full interlace
 113     template <typename T> _SelectedValueIterator<T>
 114     getSelectedValues( int component, int nbElements, int nbComponents ) const
 115     {
 116       // Values are specified by myNumbers. myValues points to the value
 117       // relating to myNumbers[0] element.
 118       const T* values = (T*) myValues;
 119       if ( myInterlace == MED_FULL_INTERLACE ) {
 120         values -= myNumbers[0] * nbComponents;
 121         values += (component - 1);
 122         return _SelectedValueIterator<T>( values, nbComponents, myNumbers );
 123       }
 124       else {
 125         values -= myNumbers[0];
 126         values += (component - 1) * nbElements;
 127         return _SelectedValueIterator<T>( values, 1, myNumbers );
 128       }
 129     }
 130     // -------------------------------------------------------------------
 131     //!< Return iterators on all values of all components
 132     template <typename T> vector< _ValueIterator<T> >
 133     getAllValues( int nbElements, int nbComponents, bool add3dComp ) const
 134     {
 135       // add the 3-d component to a vector in 2D space
 136       int addComp = ( add3dComp && nbComponents == 2 ) ? 1 : 0;
 137       vector< _ValueIterator<T> > compValIt( nbComponents + addComp);
 138       for ( int j = 1; j <= nbComponents; ++j )
 139         compValIt[ j-1 ] = getValues<T>( j, nbElements, nbComponents );
 140       return compValIt;
 141     }
 142     // -------------------------------------------------------------------
 143     //!< Return iterators on specified values for all components
 144     template <typename T> vector< _SelectedValueIterator<T> >
 145     getAllSelectedValues( int nbElements, int nbComponents, bool add3dComp ) const
 146     {
 147       // add the 3-d component to a vector in 2D space
 148       int addComp = ( add3dComp && nbComponents == 2 ) ? 1 : 0;
 149       vector< _SelectedValueIterator<T> > compValIt( nbComponents + addComp);
 150       for ( int j = 1; j <= nbComponents; ++j )
 151         compValIt[ j-1 ] = getSelectedValues<T>( j, nbElements, nbComponents );
 152       return compValIt;
 153     }
 154   };
 155
 156   //=======================================================================
 157   /*!
 158     use file reader to get values of a _SubPart
 159   */
 160   template <class FileReader>
 161   void readSubPartValues( FileReader&      reader,
 162                           const FIELD_*    field,
 163                           _SubPartValues & subValues)
 164   {
 165     medEntityMesh entity = field->getSupport()->getEntity();
 166     int     nbComponents = field->getNumberOfComponents();
 167
 168     int nbElements = 0;
 169     if ( entity == MED_NODE )
 170       nbElements = subValues.mySubPart.myNbNodes;
 171     else
 172       nbElements = subValues.mySubPart.myNbCells;
 173
 174     int nbValues = nbElements * nbComponents;
 175     if ( nbValues < 1 )
 176       return;
 177
 178     const char* undefValue = 0;
 179     set<int>* undefIndices = &subValues.myUndefIndices;
 180     vector<int>* partial   = 0;
 181     if ( !subValues.myUndefValue.empty() )
 182       undefValue = subValues.myUndefValue.c_str();
 183     if ( !subValues.myPartialIndices.empty() )
 184       partial = &subValues.myPartialIndices;
 185
 186     if ( field->getValueType() == MED_REEL64 )
 187       subValues.myValues = reader.template convertReals<double>( nbValues,
 188                                                                  undefValue, undefIndices,
 189                                                                  partial, nbComponents );
 190     else if ( field->getValueType() == MED_INT32 )
 191       subValues.myValues = reader.template convertReals<int>   ( nbValues,
 192                                                                  undefValue, undefIndices,
 193                                                                  partial, nbComponents );
 194     else
 195       subValues.myValues = reader.template convertReals<long>  ( nbValues,
 196                                                                  undefValue, undefIndices,
 197                                                                  partial, nbComponents );
 198     // convert partial indices into undef indices
 199     if ( !subValues.myPartialIndices.empty() )
 200     {
 201       // sort partial
 202       set<int> sortedPartial;
 203       vector<int>::iterator p = partial->begin(), pEnd = partial->end();
 204       while ( p != pEnd )
 205         sortedPartial.insert( sortedPartial.end(), *p++ );
 206       partial->clear();
 207
 208       // fill undef indices
 209       set<int> & undef = subValues.myUndefIndices;
 210       set<int>::iterator sp = sortedPartial.begin();
 211       int i = 1;
 212       for ( ; sp != sortedPartial.end() && i <= nbElements; ++i, ++sp ) {
 213         while ( i < *sp )
 214           undef.insert( undef.end(), i++ );
 215       }
 216       while ( i <= nbElements )
 217         undef.insert( undef.end(), i++ );
 218     }
 219   }
 220   //=======================================================================
 221   /*!
 222     set values to an array
 223   */
 224   template <class T, class TMailIter, class TArray >
 225   void _setValuesToArray( vector< _ValueIterator<T> > & compValIt,
 226                           const int                     nbValues,
 227                           TArray*                       array,
 228                           TMailIter&                    cellIt,
 229                           _Support*                     interSupport,
 230                           const set<int>&               localUndefIndices,
 231                           set<int> &                    globalUndefIndices)
 232   {
 233     int nbComponents = array->getDim();
 234     if ( localUndefIndices.empty() )
 235     {
 236       for ( int c = 0; c < nbValues; ++c, ++cellIt )
 237       {
 238         int i = interSupport->getIndex( *cellIt );
 239         for ( int j = 1; j <= nbComponents; ++j )
 240           array->setIJ( i, j, compValIt[ j ].next() );
 241       }
 242     }
 243     else
 244     {
 245       // set values from the first till undefined end
 246       set<int>::const_iterator undef = localUndefIndices.begin();
 247       int c = 1, last = min( nbValues, *localUndefIndices.rbegin() );
 248       for ( ; c <= last; ++c, ++cellIt )
 249       {
 250         int i = interSupport->getIndex( *cellIt );
 251         if ( c == *undef ) {
 252           globalUndefIndices.insert( i );
 253           ++undef;
 254         }
 255         for ( int j = 1; j <= nbComponents; ++j )
 256           array->setIJ( i, j, compValIt[ j ].next() );
 257       }
 258       // set all the rest values
 259       for ( ; c <= nbValues; ++c, ++cellIt )
 260       {
 261         int i = interSupport->getIndex( *cellIt );
 262         for ( int j = 1; j <= nbComponents; ++j )
 263           array->setIJ( i, j, compValIt[ j ].next() );
 264       }
 265     }
 266   }
 267
 268   //================================================================================
 269   /*!
 270    * \brief creates a new support excluding undef indices
 271    */
 272   //================================================================================
 273
 274   SUPPORT* makeShiftedSupport(const SUPPORT* support,
 275                               const set<int> undefIndices)
 276   {
 277     int nbElements = support->getNumberOfElements(MED_ALL_ELEMENTS);
 278     int nbTypes    = support->getNumberOfTypes();
 279
 280     int i, shitf = 1;
 281     set<int>::const_iterator undef;
 282
 283     // shift support numbers
 284     int * index  = new int[ nbTypes + 1 ];
 285     int * number = new int[ nbElements - undefIndices.size() ];
 286     if ( support->isOnAllElements() ) {
 287       // make shifted number
 288       int * numberPtr = number;
 289       for ( i = 0, undef = undefIndices.begin(); undef != undefIndices.end(); ++undef )
 290         while ( ++i < *undef )
 291           *numberPtr++ = i;
 292       while ( ++i <= nbElements )
 293         *numberPtr++ = i;
 294       // fill index but without shift
 295       const int * oldNbElemByType = support->getNumberOfElements();
 296       index[0] = 1;
 297       for ( int t = 0; t < nbTypes; ++t )
 298         index[ t+1 ] = index[ t ] + oldNbElemByType[ t ];
 299     }
 300     else {
 301       // shift existing number
 302       shitf = 1;
 303       const int * oldNumber = support->getNumber( MED_ALL_ELEMENTS );
 304       std::copy( oldNumber, oldNumber + *undefIndices.begin()-1, number ); // copy [0:firstUndef]
 305       for ( undef = undefIndices.begin(); undef != undefIndices.end(); ) {
 306         i = *undef++;
 307         int nextUndef = ( undef != undefIndices.end() ) ? *undef : nbElements + 1;
 308         while ( ++i < nextUndef )
 309           number[ i-1 - shitf ] = oldNumber[ i-1 ]; // shift number
 310         ++shitf;
 311       }
 312       // copy index
 313       const int * oldIndex  = support->getNumberIndex();
 314       std::copy( oldIndex, oldIndex + nbTypes + 1, index );
 315     }
 316     // update index
 317     {
 318       set<int>::const_reverse_iterator undef = undefIndices.rbegin();
 319       for ( ; undef != undefIndices.rend(); ++undef ) {
 320         for ( int t = nbTypes; t ; --t )
 321           if ( *undef < index[ t ] )
 322             --index[ t ];
 323           else
 324             break;
 325       }
 326     }
 327     // count number of elements by type
 328     vector<int> nbElem( nbTypes );
 329     for ( int t = 0; t < nbTypes; ++t )
 330       nbElem[ t ] = index[ t+1 ] - index[ t ];
 331
 332     SUPPORT* newSup = new SUPPORT();
 333     newSup->setMesh                 ( support->getMesh() );
 334     newSup->setNumberOfGeometricType( nbTypes );
 335     newSup->setGeometricType        ( support->getTypes() );
 336     newSup->setNumberOfElements     ( &nbElem[0] );
 337     newSup->setNumber               ( index, number, /*shallowCopy =*/ true );
 338     newSup->setEntity               ( support->getEntity() );
 339     newSup->setAll                  ( false );
 340
 341     return newSup;
 342   }
 343
 344   //=======================================================================
 345   /*!
 346     set values to a field
 347   */
 348   template <class T, class INTERLACE>
 349   void _setValuesToField( FIELD<T,INTERLACE>*    field,
 350                           _Support*              interSupport,
 351                           list<_SubPartValues> & subPartValues )
 352   {
 353     medEntityMesh entity = field->getSupport()->getEntity();
 354     SUPPORT* support = interSupport->medSupport( entity );
 355     field->setSupport( support );
 356
 357     int j, nbComponents = field->getNumberOfComponents();
 358     int      nbElements = field->getSupport()->getNumberOfElements(MED_ALL_ELEMENTS);
 359
 360     typedef typename MEDMEM_ArrayInterface<T,INTERLACE,NoGauss>::Array TArray;
 361
 362     const string& constValue = subPartValues.front().myConstValue;
 363     if ( !constValue.empty() )
 364     {
 365       vector<T> values(nbComponents * nbElements, (T)atof(constValue.c_str()));
 366       field->setArray( new TArray( &values[0], nbComponents, nbElements ));
 367       field->setNumberOfValues( nbElements );
 368       return;
 369     }
 370
 371     field->allocValue( nbComponents, nbElements );
 372     TArray * array = field->getArrayNoGauss();
 373
 374     set<int> undefIndices; // indices of undefined values
 375
 376     // loop on subParts conatining ensight field values
 377     list<_SubPartValues>::iterator subValue = subPartValues.begin();
 378     for ( ; subValue != subPartValues.end(); ++subValue )
 379     {
 380       _SubPart & subPart = subValue->mySubPart;
 381       int nbValues = entity==MED_NODE ? subPart.myNbNodes : subPart.myNbCells;
 382       if ( nbValues == 0 )
 383         continue;
 384       // value iterator for each component
 385       vector< _ValueIterator<T> > compValIt( nbComponents+1 );
 386       for ( j = 1; j <= nbComponents; ++j )
 387         compValIt[ j ] = subValue->template getValues<T>( j, nbValues, nbComponents );
 388
 389       // Set values
 390
 391       if ( entity == MED_NODE ) {
 392         map< int, _noeud >::iterator inode = subPart.myFirstNode;
 393         _setValuesToArray( compValIt, nbValues, array, inode, interSupport,
 394                            subValue->myUndefIndices, undefIndices );
 395       }
 396       else {
 397         _groupe::TMailleIter cell = subPart.myFirstCell;
 398         _setValuesToArray( compValIt, nbValues, array, cell, interSupport,
 399                            subValue->myUndefIndices, undefIndices );
 400       }
 401     }
 402
 403     if ( !undefIndices.empty() )
 404     {
 405       // Some values are undefined; it's necessary to compact values in the array
 406       // and create a new support of defined numbers only.
 407
 408       // shift values
 409       int i, shitf = 1;
 410       set<int>::iterator undef = undefIndices.begin();
 411       while ( undef != undefIndices.end() ) {
 412         i = *undef++;
 413         int nextUndef = ( undef != undefIndices.end() ) ? *undef : nbElements + 1;
 414         while ( ++i < nextUndef ) {
 415           for ( int j = 1; j <= nbComponents; ++j )
 416             array->setIJ( i - shitf, j, array->getIJ( i, j )); // shift value
 417           //number[ i-1 - shitf ] = oldNumber[ i-1 ]; // shift number
 418         }
 419         ++shitf;
 420       }
 421       array->_nbelem -= undefIndices.size();
 422       array->_arraySize = array->_nbelem * nbComponents;
 423
 424       support = makeShiftedSupport( support, undefIndices );
 425       support->setName( STRING("Partial for ") << field->getName() );
 426       field->setSupport( support );
 427       field->setNumberOfValues( array->_nbelem );
 428     }
 429     if ( support->getName().empty() )
 430       support->setName( STRING("Support for ") << field->getName() );
 431   }
 432
 433   //=======================================================================
 434   /*!
 435     set values to a field
 436   */
 437   void setValuesToField( FIELD_*                field,
 438                          _Support*              interSupport,
 439                          list<_SubPartValues> & subPartValues )
 440   {
 441     switch ( field->getInterlacingType() ) {
 442     case MED_FULL_INTERLACE:
 443       if ( field->getValueType() == MED_REEL64 )
 444         _setValuesToField( static_cast< FIELD<double, FullInterlace>* >( field ),
 445                            interSupport, subPartValues );
 446       else
 447         _setValuesToField( static_cast< FIELD<int, FullInterlace>* >( field ),
 448                            interSupport, subPartValues );
 449       break;
 450     case MED_NO_INTERLACE:
 451       if ( field->getValueType() == MED_REEL64 )
 452         _setValuesToField( static_cast< FIELD<double, NoInterlace>* >( field ),
 453                            interSupport, subPartValues );
 454       else
 455         _setValuesToField( static_cast< FIELD<int, NoInterlace>* >( field ),
 456                            interSupport, subPartValues );
 457       break;
 458     case MED_NO_INTERLACE_BY_TYPE:
 459       if ( field->getValueType() == MED_REEL64 )
 460         _setValuesToField( static_cast< FIELD<double, NoInterlaceByType>* >( field ),
 461                            interSupport, subPartValues );
 462       else
 463         _setValuesToField( static_cast< FIELD<int, NoInterlaceByType>* >( field ),
 464                            interSupport, subPartValues );
 465       break;
 466     default:;
 467     }
 468   }
 469
 470
 471   //================================================================================
 472   /*!
 473    * \brief Return node numbers of a non-nodal support
 474    */
 475   //================================================================================
 476
 477   int* getNodeNumbers( const SUPPORT* support, int & nbNodes )
 478   {
 479     map<int, int> nodeNumbers;
 480     _CaseFileDriver_User::getSupportNodes( support, nodeNumbers );
 481     nbNodes = nodeNumbers.size();
 482     int* nNumbers = new int[ nbNodes ];
 483     int* nPtr = nNumbers;
 484     map<int, int>::iterator n = nodeNumbers.begin(), nEnd = nodeNumbers.end();
 485     while ( n != nEnd )
 486       *nPtr++ = n->first, n++;
 487
 488     return nNumbers;
 489   }
 490
 491 } // namespace
 492
 493 //=======================================================================
 494 /*!
 495   Constructor.
 496 */
 497 ENSIGHT_FIELD_DRIVER::ENSIGHT_FIELD_DRIVER():
 498   _CaseFileDriver_User(),
 499   _ptrField((FIELD_ *) 0)
 500 {}
 501 //=======================================================================
 502 /*!
 503   Constructor.
 504 */
 505 ENSIGHT_FIELD_DRIVER::ENSIGHT_FIELD_DRIVER(const string &  fileName,
 506                                            FIELD_       *  ptrField,
 507                                            med_mode_acces  accessMode):
 508   _CaseFileDriver_User(fileName,accessMode),
 509   _ptrField((FIELD_ *) ptrField)
 510 {
 511 }
 512 //=======================================================================
 513 /*!
 514   Copy constructor.
 515 */
 516 ENSIGHT_FIELD_DRIVER::ENSIGHT_FIELD_DRIVER(const ENSIGHT_FIELD_DRIVER & fieldDriver):
 517   _CaseFileDriver_User(fieldDriver),
 518   _ptrField(fieldDriver._ptrField),
 519   _fieldName(fieldDriver._fieldName)
 520 {
 521 }
 522 //=======================================================================
 523 /*!
 524   Take missing data from other driver.
 525 */
 526 void ENSIGHT_FIELD_DRIVER::merge ( const GENDRIVER& driver )
 527 {
 528   _CaseFileDriver_User::merge( driver );
 529
 530   const ENSIGHT_FIELD_DRIVER* other = dynamic_cast< const ENSIGHT_FIELD_DRIVER* >( &driver );
 531   if ( other ) {
 532     if ( !_ptrField )
 533       _ptrField = other->_ptrField;
 534     if ( _constantValue.empty() )
 535       _constantValue = other->_constantValue;
 536     // _fieldName is copied by GENDRIVER
 537   }
 538 }
 539 //=======================================================================
 540 /*!
 541   Destructor.
 542 */
 543 ENSIGHT_FIELD_DRIVER::~ENSIGHT_FIELD_DRIVER()
 544 {
 545 }
 546 //=======================================================================
 547 /*!
 548   Set the name of the FIELD in EnSight file
 549 */
 550 void ENSIGHT_FIELD_DRIVER::setFieldName(const string & fieldName) throw (MEDEXCEPTION)
 551 {
 552   const char* LOC = "ENSIGHT_FIELD_DRIVER::setFieldName(): ";
 553   if ( fieldName.size() > MAX_FIELD_NAME_LENGTH )
 554     throw MEDEXCEPTION( compatibilityPb(LOC) << "too long name (> " <<
 555                         MAX_FIELD_NAME_LENGTH << "): " << fieldName);
 556
 557   // look for illegal characters
 558   string::size_type pos = fieldName.find_first_of( ILLEGAL_FIELD_NAME_CHARACTERS );
 559   if ( pos != fieldName.npos )
 560     throw MEDEXCEPTION( compatibilityPb(LOC) << "Character " << pos << " (" << fieldName[pos] <<
 561                         ") "<< " in " << fieldName << " is not allowed in field name in EnSight\n"
 562                         "The illegal characters are: '" << ILLEGAL_FIELD_NAME_CHARACTERS << "'");
 563   _fieldName = fieldName;
 564 }
 565 //=======================================================================
 566 /*!
 567   Check possibility to open a case file or a data file
 568 */
 569 void ENSIGHT_FIELD_DRIVER::openConst(bool checkDataFile) const throw (MEDEXCEPTION)
 570 {
 571   const char * LOC ="ENSIGHT_FIELD_DRIVER::open() : ";
 572   BEGIN_OF_MED(LOC);
 573
 574   if ( checkDataFile )
 575   {
 576     if ( !getConstantValue().empty() )
 577       return; // constant value is either stored in case file or is read by _CaseFileDriver
 578
 579     if ( getDataFileName().empty() )
 580       throw MED_EXCEPTION
 581         ( LOCALIZED( STRING(LOC) << "Internal error, variable file name is empty"));
 582
 583     if (!canOpenFile( getDataFileName(), getAccessMode() ))
 584       throw MED_EXCEPTION
 585         ( LOCALIZED( STRING(LOC) << "Can not open Ensight Variable file " << getDataFileName()
 586                      << " in access mode " << getAccessMode()));
 587   }
 588   else
 589   {
 590     if ( getCaseFileName().empty() )
 591       throw MED_EXCEPTION
 592         ( LOCALIZED( STRING(LOC) << "Case file name is empty, "
 593                      "please set a correct fileName before calling open()"));
 594
 595     if ( !canOpenFile( getCaseFileName(), getAccessMode() ))
 596       throw MED_EXCEPTION
 597         ( LOCALIZED( STRING(LOC) << "Can not open Ensight Case file " << getCaseFileName()
 598                      << " in access mode " << getAccessMode()));
 599   }
 600
 601   END_OF_MED(LOC);
 602 }
 603 // ==========================================================================================
 604 //function : ENSIGHT_FIELD_RDONLY_DRIVER
 605
 606 ENSIGHT_FIELD_RDONLY_DRIVER::ENSIGHT_FIELD_RDONLY_DRIVER()
 607   :ENSIGHT_FIELD_DRIVER(), _fieldStep(1)
 608 {
 609 }
 610
 611 //=======================================================================
 612 /*!
 613   Constructor to read a field of index-th time step
 614 */
 615 //=======================================================================
 616
 617 ENSIGHT_FIELD_RDONLY_DRIVER::ENSIGHT_FIELD_RDONLY_DRIVER(const string & fileName,
 618                                                          FIELD_ *       ptrField,
 619                                                          int            step)
 620   :ENSIGHT_FIELD_DRIVER(fileName,ptrField,MED_EN::RDONLY),_fieldStep(step)
 621 {
 622 }
 623
 624 //=======================================================================
 625 //function : ENSIGHT_FIELD_RDONLY_DRIVER
 626 //=======================================================================
 627
 628 ENSIGHT_FIELD_RDONLY_DRIVER::ENSIGHT_FIELD_RDONLY_DRIVER(const ENSIGHT_FIELD_RDONLY_DRIVER & drv)
 629   :ENSIGHT_FIELD_DRIVER(drv),_fieldStep(drv._fieldStep)
 630 {
 631 }
 632
 633 //=======================================================================
 634 //function : ~ENSIGHT_FIELD_RDONLY_DRIVER
 635 //=======================================================================
 636
 637 ENSIGHT_FIELD_RDONLY_DRIVER::~ENSIGHT_FIELD_RDONLY_DRIVER()
 638 {
 639   close();
 640 }
 641
 642 //=======================================================================
 643 //function : copy
 644 //=======================================================================
 645
 646 GENDRIVER * ENSIGHT_FIELD_RDONLY_DRIVER::copy(void) const
 647 {
 648   ENSIGHT_FIELD_RDONLY_DRIVER * myDriver = new ENSIGHT_FIELD_RDONLY_DRIVER(*this);
 649
 650   return myDriver ;
 651 }
 652 //=======================================================================
 653 //  Take missing data from other driver.
 654 //=======================================================================
 655
 656 void ENSIGHT_FIELD_RDONLY_DRIVER::merge ( const GENDRIVER& driver )
 657 {
 658   ENSIGHT_FIELD_DRIVER::merge( driver );
 659
 660   const ENSIGHT_FIELD_RDONLY_DRIVER* other =
 661     dynamic_cast< const ENSIGHT_FIELD_RDONLY_DRIVER* >( &driver );
 662   if ( other ) {
 663     if ( _fieldStep < other->_fieldStep )
 664       _fieldStep = other->_fieldStep;
 665   }
 666 }
 667
 668 //=======================================================================
 669 //function : read
 670 //=======================================================================
 671
 672 void ENSIGHT_FIELD_RDONLY_DRIVER::read (void)
 673   throw (MEDEXCEPTION)
 674 {
 675   const char * LOC = "ENSIGHT_FIELD_RDONLY_DRIVER::read() : " ;
 676   BEGIN_OF_MED(LOC);
 677
 678   openConst(false); // check if can read case file
 679
 680   _CaseFileDriver caseFile( getCaseFileName(), this);
 681   caseFile.read();
 682
 683   // find out index of variable to read
 684   int variableIndex = caseFile.getVariableIndex( _fieldName );
 685   if ( !variableIndex )
 686     variableIndex = caseFile.getVariableIndex( _ptrField->getName() );
 687   if ( !variableIndex ) {
 688     if ( !_fieldName.empty() )
 689       throw MEDEXCEPTION
 690         (LOCALIZED(STRING(LOC) << "no field found by name |" << _fieldName << "|"));
 691     else
 692       throw MEDEXCEPTION
 693         (LOCALIZED(STRING(LOC) << "no field found by name |" << _ptrField->getName() << "|"));
 694   }
 695
 696   //  here data from Case File is passed:
 697   // * field name
 698   // * number of components
 699   // * etc.
 700   caseFile.setDataFileName( variableIndex, _fieldStep, this );
 701
 702   openConst(true); // check if can read data file
 703
 704   getInterData(); // get data on nb of entities by type and on their relocation
 705
 706   cout << "-> Entering into the field file " << getDataFileName() << endl  ;
 707
 708   if ( !getConstantValue().empty() )
 709   {
 710     // Create a constant value field
 711
 712     medEntityMesh entity = MED_CELL;
 713     GROUP* support = new GROUP;
 714     support->setName( string("SupportOnAll_")+entNames[entity]);
 715     support->setMesh( getInterData()->_medMesh );
 716     support->setAll( true );
 717     support->setEntity( entity );
 718     support->update();
 719
 720     _groupe interGroup;
 721     interGroup.medGroup = support;
 722     _Support interSupport;
 723     interSupport.setGroup( &interGroup );
 724
 725     list<_SubPartValues> subPartValues( 1 );
 726     subPartValues.back().myConstValue = getConstantValue();
 727
 728     setValuesToField( _ptrField, &interSupport, subPartValues );
 729   }
 730   else
 731   {
 732     // Read values
 733
 734     if ( isBinaryDataFile( getDataFileName() ) ) // binary
 735     {
 736       if ( isGoldFormat() ) // Gold
 737       {
 738         readGoldBinary();
 739       }
 740       else // EnSight6
 741       {
 742         read6Binary();
 743       }
 744     }
 745     else // ASCII
 746     {
 747       if ( isGoldFormat() ) // Gold
 748       {
 749         readGoldASCII();
 750       }
 751       else // EnSight6
 752       {
 753         read6ASCII();
 754       }
 755     }
 756   }
 757 }
 758
 759 //=======================================================================
 760 /*!
 761   Read Ensight Gold ASCII
 762 */
 763 //=======================================================================
 764
 765 void ENSIGHT_FIELD_RDONLY_DRIVER::readGoldASCII ()
 766 {
 767
 768   // data that was set by CaseFile Driver
 769   medEntityMesh entity = _ptrField->getSupport()->getEntity();
 770
 771   _SupportDesc         supportDescriptor;
 772   list<_SubPartValues> subPartValues;
 773
 774   _ASCIIFileReader varFile( getDataFileName() );
 775
 776   if ( isSingleFileMode() ) {
 777     int curTimeStep = 1;
 778     while ( curTimeStep < getIndexInDataFile() ) {
 779       while ( !varFile.isTimeStepEnd() )
 780         varFile.getLine();
 781       curTimeStep++;
 782     }
 783     while ( !varFile.isTimeStepBeginning() )
 784       varFile.getLine();
 785   }
 786   string description = varFile.getLine(); // description line
 787   _ptrField->setDescription( description );
 788
 789   int partNumber = 0;
 790   _SubPartValues subValuesSample; // to keep interlace
 791   subValuesSample.myInterlace = MED_NO_INTERLACE;
 792   while ( !varFile.isTimeStepEnd() )
 793   {
 794     string word, restLine, line = varFile.getLine();
 795     TStrTool::split( line, word, restLine );
 796     if ( word == "part" ) {
 797       partNumber = varFile.getInt(); // part #
 798       continue;
 799     }
 800     subPartValues.push_back( subValuesSample );
 801     _SubPartValues & subValues = subPartValues.back();
 802
 803     if ( restLine == "undef" )
 804       subValues.myUndefValue = varFile.getWord(); // undef
 805
 806     if ( restLine == "partial" ) {
 807       int nbPartial = varFile.getInt(); // ne
 808       subValues.myPartialIndices.reserve( nbPartial );
 809       while ( nbPartial-- )
 810         subValues.myPartialIndices.push_back( varFile.getInt() ); // i_ne
 811     }
 812     _SubPartDesc desc( partNumber, word );
 813     subValues.mySubPart = *getSubPart( desc );
 814     readSubPartValues( varFile, _ptrField, subValues );
 815     supportDescriptor.insert( desc );
 816   }
 817
 818   _Support* interSupport = getSupport( supportDescriptor, entity );
 819   setValuesToField( _ptrField, interSupport, subPartValues );
 820
 821 }
 822
 823 //=======================================================================
 824 /*!
 825   Read Ensight GOLD binary
 826 */
 827 //=======================================================================
 828
 829 void ENSIGHT_FIELD_RDONLY_DRIVER::readGoldBinary()
 830 {
 831
 832   // data that was set by CaseFile Driver
 833   medEntityMesh entity = _ptrField->getSupport()->getEntity();
 834
 835   _SupportDesc         supportDescriptor;
 836   list<_SubPartValues> subPartValues;
 837
 838   _BinaryFileReader varFile( getDataFileName() );
 839
 840   // check if swapping bytes needed
 841   try {
 842     skipTimeStamp( varFile );
 843   }
 844   catch (...) {
 845     varFile.swapBytes();
 846     varFile.rewind();
 847   }
 848   if ( getIndexInDataFile() <= 1 )
 849     varFile.rewind();
 850
 851   if ( isSingleFileMode() ) {
 852     // one time step may be skipped by skipTimeStamp()
 853     int curTimeStep = varFile.getPosition() ? 2 : 1 ;
 854     while ( curTimeStep < getIndexInDataFile() ) {
 855       skipTimeStamp( varFile );
 856       curTimeStep++;
 857     }
 858     varFile.skipTimeStepBeginning();
 859   }
 860   TStrOwner description( varFile.getLine() ); // description line
 861   _ptrField->setDescription( description.myValues );
 862
 863   int partNumber = 0;
 864   _SubPartValues subValuesSample; // to keep interlace
 865   subValuesSample.myInterlace = MED_NO_INTERLACE;
 866   while ( !varFile.eof() )
 867   {
 868     TStrOwner line( varFile.getLine() );
 869     if ( isTimeStepEnd( line.myValues ))
 870       break;
 871     string word, restLine;
 872     TStrTool::split( line.myValues, word, restLine );
 873     if ( word == "part" ) {
 874       partNumber = *TIntOwner( varFile.getInt(1)); // part #
 875       continue;
 876     }
 877     subPartValues.push_back( subValuesSample );
 878     _SubPartValues & subValues = subPartValues.back();
 879
 880     if ( restLine == "undef" )
 881       subValues.myUndefValue = (STRING( *TFltOwner( varFile.getFlt(1) ))); // undef
 882
 883     if ( restLine == "partial" ) {
 884       int nbPartial = *TIntOwner( varFile.getInt(1) ); // ne
 885       TIntOwner partial( varFile.getInt( nbPartial ));
 886       subValues.myPartialIndices.assign( partial.myValues,
 887                                          partial.myValues+nbPartial ); // i_ne
 888     }
 889     _SubPartDesc desc( partNumber, word );
 890     subValues.mySubPart = *getSubPart( desc );
 891     readSubPartValues( varFile, _ptrField, subValues );
 892     supportDescriptor.insert( desc );
 893   }
 894
 895   _Support* interSupport = getSupport( supportDescriptor, entity );
 896   setValuesToField( _ptrField, interSupport, subPartValues );
 897
 898 }
 899
 900 //=======================================================================
 901 /*!
 902   Read Ensight6 ASCII
 903 */
 904 //=======================================================================
 905
 906 void ENSIGHT_FIELD_RDONLY_DRIVER::read6ASCII()
 907 {
 908
 909   // data that was set by CaseFile Driver
 910   medEntityMesh entity = _ptrField->getSupport()->getEntity();
 911
 912   _SupportDesc         supportDescriptor;
 913   list<_SubPartValues> subPartValues;
 914
 915   _ASCIIFileReader varFile( getDataFileName() );
 916
 917   if ( isSingleFileMode() ) {
 918     int curTimeStep = 1;
 919     while ( curTimeStep < getIndexInDataFile() ) {
 920       while ( !varFile.isTimeStepEnd() )
 921         varFile.getLine();
 922       curTimeStep++;
 923     }
 924     while ( !varFile.isTimeStepBeginning() )
 925       varFile.getLine();
 926   }
 927   string description = varFile.getLine(); // description line
 928   _ptrField->setDescription( description );
 929
 930   const _SubPart* subPart = 0;
 931   if ( entity == MED_NODE ) // variable per node
 932   {
 933     _SubPartDesc desc = _SubPartDesc::globalCoordDesc();
 934     subPart = getSubPart( desc );
 935     if ( subPart ) {
 936       supportDescriptor.insert( desc );
 937       _SubPartValues subValues( *subPart );
 938       subValues.myInterlace = MED_FULL_INTERLACE;
 939       readSubPartValues( varFile, _ptrField, subValues );
 940       subPartValues.push_back( subValues );
 941     }
 942   }
 943   int partNumber = 0;
 944   while ( !varFile.isTimeStepEnd() )
 945   {
 946     string word = varFile.getWord();
 947     if ( word == "part" ) {
 948       partNumber = varFile.getInt();
 949       continue;
 950     }
 951     _SubPartDesc desc( partNumber, word );
 952     supportDescriptor.insert( desc );
 953     subPart = getSubPart( desc );
 954     _SubPartValues subValues( *subPart );
 955     if ( desc.typeName() == "block" )
 956       subValues.myInterlace = MED_NO_INTERLACE;
 957     else
 958       subValues.myInterlace = MED_FULL_INTERLACE;
 959     readSubPartValues( varFile, _ptrField, subValues );
 960     subPartValues.push_back( subValues );
 961   }
 962
 963   _Support* interSupport = getSupport( supportDescriptor, entity );
 964   setValuesToField( _ptrField, interSupport, subPartValues );
 965
 966 }
 967
 968 //=======================================================================
 969 /*!
 970   Read Ensight6 binary
 971 */
 972 //=======================================================================
 973
 974 void ENSIGHT_FIELD_RDONLY_DRIVER::read6Binary()
 975 {
 976
 977   // data that was set by CaseFile Driver
 978   medEntityMesh entity = _ptrField->getSupport()->getEntity();
 979
 980   _SupportDesc         supportDescriptor;
 981   list<_SubPartValues> subPartValues;
 982   const _SubPart*      subPart = 0;
 983   int                  partNumber = 0;
 984
 985   _BinaryFileReader varFile( getDataFileName() );
 986
 987   // check if swapping bytes needed
 988   try {
 989     skipTimeStamp( varFile );
 990   }
 991   catch (...) {
 992     varFile.swapBytes();
 993     varFile.rewind();
 994   }
 995   if ( getIndexInDataFile() <= 1 )
 996     varFile.rewind();
 997
 998   if ( isSingleFileMode() ) {
 999     // one time step may be skipped by skipTimeStamp()
1000     int curTimeStep = varFile.getPosition() ? 2 : 1 ;
1001     while ( curTimeStep < getIndexInDataFile() ) {
1002       skipTimeStamp( varFile );
1003       curTimeStep++;
1004     }
1005     varFile.skipTimeStepBeginning();
1006   }
1007
1008   TStrOwner description( varFile.getLine() ); // description line
1009   _ptrField->setDescription( description.myValues );
1010
1011   if ( entity == MED_NODE ) // global unstructured node values
1012   {
1013     _SubPartDesc desc = _SubPartDesc::globalCoordDesc();
1014     subPart = getSubPart( desc );
1015     if ( subPart ) {
1016       supportDescriptor.insert( desc );
1017       _SubPartValues subValues( *subPart );
1018       subValues.myInterlace = MED_FULL_INTERLACE;
1019       readSubPartValues( varFile, _ptrField, subValues );
1020       subPartValues.push_back( subValues );
1021     }
1022   }
1023   while ( !varFile.eof() )
1024   {
1025     TStrOwner line( varFile.getLine() );
1026     if ( isTimeStepEnd( line ))
1027       break;
1028     string word, restLine;
1029     TStrTool::split( line.myValues, word, restLine );
1030     if ( word == "part" ) {
1031       partNumber = atoi( restLine.c_str() );
1032       continue;
1033     }
1034     _SubPartDesc desc( partNumber, word );
1035     supportDescriptor.insert( desc );
1036     subPart = getSubPart( desc );
1037     _SubPartValues subValues( *subPart );
1038     if ( desc.typeName() == "block" )
1039       subValues.myInterlace = MED_NO_INTERLACE;
1040     else
1041       subValues.myInterlace = MED_FULL_INTERLACE;
1042     readSubPartValues( varFile, _ptrField, subValues );
1043     subPartValues.push_back( subValues );
1044   }
1045
1046   _Support* interSupport = getSupport( supportDescriptor, entity );
1047   setValuesToField( _ptrField, interSupport, subPartValues );
1048
1049 }
1050
1051 //================================================================================
1052 /*!
1053  * \brief Looks for beginning of the next time stamp
1054  */
1055 //================================================================================
1056
1057 void ENSIGHT_FIELD_RDONLY_DRIVER::skipTimeStamp(_BinaryFileReader& varFile)
1058 {
1059   medEntityMesh entity = _ptrField->getSupport()->getEntity();
1060   int     nbComponents = _ptrField->getNumberOfComponents();
1061
1062   if ( isSingleFileMode() )
1063     varFile.skipTimeStepBeginning();
1064
1065   varFile.skip( 80 ); // description
1066
1067   _SubPart* subPart;
1068   if ( entity == MED_NODE && !isGoldFormat() ) { // skip values on global nodes
1069     subPart = getSubPart( _SubPartDesc::globalCoordDesc() );
1070     if ( subPart )
1071       varFile.skip( subPart->myNbNodes * nbComponents * sizeof(float) );
1072   }
1073   int partNumber;
1074   while ( !varFile.eof() )
1075   {
1076     TStrOwner line( varFile.getLine() );
1077     if ( isTimeStepEnd( line ))
1078       break;
1079
1080     string word, restLine;
1081     TStrTool::split( line.myValues, word, restLine );
1082
1083     if ( word == "part" ) {
1084       if ( isGoldFormat() )
1085         partNumber = *TIntOwner( varFile.getInt(1));
1086       else
1087         partNumber = atoi( restLine.c_str() );
1088       continue;
1089     }
1090     if ( restLine == "undef" ) {
1091       varFile.skip( sizeof(int) );
1092     }
1093     if ( restLine == "partial" ) {
1094       int nbPartial = *TIntOwner( varFile.getInt(1)); // ne
1095       varFile.skip( nbPartial * sizeof(int));
1096     }
1097
1098     _SubPartDesc desc( partNumber, word );
1099     subPart = getSubPart( desc );
1100     int nbElems = ( entity == MED_NODE ) ? subPart->myNbNodes : subPart->myNbCells;
1101     varFile.skip( nbElems * nbComponents * sizeof(float) );
1102   }
1103 }
1104
1105 //=======================================================================
1106 //function : write
1107 //=======================================================================
1108
1109 void ENSIGHT_FIELD_RDONLY_DRIVER::write(void) const
1110   throw (MEDEXCEPTION)
1111 {
1112   throw MEDEXCEPTION("ENSIGHT_FIELD_RDONLY_DRIVER::write : Can't write with a RDONLY driver !");
1113 }
1114
1115 //=======================================================================
1116 //function : ENSIGHT_FIELD_WRONLY_DRIVER
1117 //=======================================================================
1118
1119 ENSIGHT_FIELD_WRONLY_DRIVER::ENSIGHT_FIELD_WRONLY_DRIVER()
1120   :ENSIGHT_FIELD_DRIVER()
1121 {
1122 }
1123
1124 //=======================================================================
1125 //function : ENSIGHT_FIELD_WRONLY_DRIVER
1126 //=======================================================================
1127
1128 ENSIGHT_FIELD_WRONLY_DRIVER::ENSIGHT_FIELD_WRONLY_DRIVER(const string & fileName,
1129                                                          FIELD_ *       ptrField)
1130   :ENSIGHT_FIELD_DRIVER(fileName,ptrField,MED_EN::WRONLY)
1131 {
1132 }
1133
1134 //=======================================================================
1135 //function : ENSIGHT_FIELD_WRONLY_DRIVER
1136 //=======================================================================
1137
1138
1139 ENSIGHT_FIELD_WRONLY_DRIVER::ENSIGHT_FIELD_WRONLY_DRIVER(const ENSIGHT_FIELD_WRONLY_DRIVER & drv)
1140   :ENSIGHT_FIELD_DRIVER(drv)
1141 {
1142 }
1143
1144 //=======================================================================
1145 //function : ~ENSIGHT_FIELD_WRONLY_DRIVER
1146 //=======================================================================
1147
1148 ENSIGHT_FIELD_WRONLY_DRIVER::~ENSIGHT_FIELD_WRONLY_DRIVER()
1149 {
1150   close();
1151 }
1152
1153
1154 //=======================================================================
1155 //function : copy
1156 //=======================================================================
1157
1158 GENDRIVER * ENSIGHT_FIELD_WRONLY_DRIVER::copy(void) const
1159 {
1160   ENSIGHT_FIELD_WRONLY_DRIVER * myDriver = new ENSIGHT_FIELD_WRONLY_DRIVER(*this);
1161
1162   return myDriver ;
1163 }
1164
1165 //=======================================================================
1166 //function : read
1167 //=======================================================================
1168
1169 void ENSIGHT_FIELD_WRONLY_DRIVER::read (void)
1170   throw (MEDEXCEPTION)
1171 {
1172   throw MEDEXCEPTION("ENSIGHT_FIELD_WRONLY_DRIVER::read : Can't read with a write only!");
1173 }
1174
1175 //================================================================================
1176 /*!
1177  * \brief Return pointer to the value of the i-th element
1178  */
1179 //================================================================================
1180
1181 const char* getValuePointer( int i, const FIELD_* field )
1182 {
1183   switch ( field->getInterlacingType() ) {
1184   case MED_FULL_INTERLACE:
1185     if ( field->getValueType() == MED_REEL64 )
1186       return (const char*) & static_cast<const FIELD<double, FullInterlace>* >
1187         ( field )->getArrayNoGauss()->getIJ( i, 1 );
1188     else
1189       return (const char*) & static_cast<const FIELD<int, FullInterlace>* >
1190         ( field )->getArrayNoGauss()->getIJ( i, 1 );
1191   case MED_NO_INTERLACE:
1192     if ( field->getValueType() == MED_REEL64 )
1193       return (const char*) & static_cast<const FIELD<double, NoInterlace>* >
1194         ( field )->getArrayNoGauss()->getIJ( i, 1 );
1195     else
1196       return (const char*) & static_cast<const FIELD<int, NoInterlace>* >
1197         ( field )->getArrayNoGauss()->getIJ( i, 1 );
1198   case MED_NO_INTERLACE_BY_TYPE:
1199     if ( field->getValueType() == MED_REEL64 )
1200       return (const char*) & static_cast<const FIELD<double, NoInterlaceByType>* >
1201         ( field )->getArrayNoGauss()->getIJ( i, 1 );
1202     else
1203       return (const char*) & static_cast<const FIELD<int, NoInterlaceByType>* >
1204         ( field )->getArrayNoGauss()->getIJ( i, 1 );
1205   default:;
1206   }
1207   return 0;
1208 }
1209
1210 //=======================================================================
1211 /*!
1212  * \brief Write field values in EnSight6 ASCII format
1213 */
1214 //=======================================================================
1215
1216 template< typename T >
1217 void write6ASCIIValues( const _SubPartValues& subPartValues,
1218                         int                   nbValues,
1219                         int                   nbComponents,
1220                         int                   componentShift,
1221                         ofstream &            ensightDataFile )
1222 {
1223   // add the 3-d component to a vector in 2D space
1224   const bool add3dComp = true;
1225   int mypoint=0;
1226   if ( subPartValues.myNumbers )
1227   {
1228     vector< _SelectedValueIterator<T> > compValIt =
1229       subPartValues.template getAllSelectedValues<T>( componentShift, nbComponents, add3dComp );
1230     nbComponents = compValIt.size();
1231     // in full interlace
1232     for (int i=0; i<nbValues; i++)
1233       for(int j=0; j<nbComponents; j++) {
1234         ensightDataFile << setw(12) << _toFloat( compValIt[j].next() );
1235         if (++mypoint == 6) {
1236           ensightDataFile << endl ;
1237           mypoint=0;
1238         }
1239       }
1240   }
1241   else {
1242     vector< _ValueIterator<T> > compValIt =
1243       subPartValues.template getAllValues<T>( componentShift, nbComponents, add3dComp );
1244     nbComponents = compValIt.size();
1245     // in full interlace
1246     for (int i=0; i<nbValues; i++)
1247       for(int j=0; j<nbComponents; j++) {
1248         ensightDataFile << setw(12) << _toFloat( compValIt[j].next() );
1249         if (++mypoint == 6) {
1250           ensightDataFile << endl ;
1251           mypoint=0;
1252         }
1253       }
1254   }
1255   if ( nbValues * nbComponents % 6 )
1256     ensightDataFile << endl;
1257 }
1258
1259 //=======================================================================
1260 /*!
1261  * \brief Write field values in EnSight GOLD ASCII format
1262 */
1263 //=======================================================================
1264
1265 template< typename T >
1266 void writeGoldASCIIValues( const _SubPartValues& subPartValues,
1267                            int                   nbValues,
1268                            int                   nbComponents,
1269                            int                   componentShift,
1270                            ofstream &            ensightDataFile )
1271 {
1272   // in no-interlace
1273   for(int j=1; j<=nbComponents; j++) {
1274     if ( subPartValues.myNumbers )
1275     {
1276       _SelectedValueIterator<T> values
1277         = subPartValues.template getSelectedValues<T>( j, componentShift, nbComponents );
1278       for (int i=0; i<nbValues; i++)
1279         ensightDataFile << setw(12) << _toFloat( values.next() ) << endl;
1280     }
1281     else {
1282       _ValueIterator<T> values
1283         = subPartValues.template getValues<T>( j, componentShift, nbComponents );
1284       for (int i=0; i<nbValues; i++)
1285         ensightDataFile << setw(12) << _toFloat( values.next() )<< endl;
1286     }
1287   }
1288   // add the 3-d component to a vector in 2D space
1289   if ( nbComponents == 2 ) {
1290     for (int i=0; i<nbValues; i++)
1291       ensightDataFile << " 0.00000e+00" << endl;
1292   }
1293 }
1294
1295 //=======================================================================
1296 /*!
1297  * \brief Write field values in EnSight6 Binary format
1298 */
1299 //=======================================================================
1300
1301 template< typename T >
1302 void write6BinaryValues( const _SubPartValues& subPartValues,
1303                          int                   nbValues,
1304                          int                   nbComponents,
1305                          int                   componentShift,
1306                          _BinaryFileWriter &   ensightDataFile )
1307 {
1308   const bool add3dComp = true;
1309   if ( subPartValues.myNumbers )
1310   {
1311     vector< _SelectedValueIterator<T> > compValIt =
1312       subPartValues.template getAllSelectedValues<T>( componentShift, nbComponents, add3dComp );
1313     ensightDataFile.addReal( compValIt, nbValues, MED_FULL_INTERLACE );
1314   }
1315   else
1316   {
1317     vector< _ValueIterator<T> > compValIt =
1318       subPartValues.template getAllValues<T>( componentShift, nbComponents, add3dComp );
1319     ensightDataFile.addReal( compValIt, nbValues, MED_FULL_INTERLACE );
1320   }
1321 }
1322
1323 //=======================================================================
1324 /*!
1325  * \brief Write field values in EnSight GOLD Binary format
1326 */
1327 //=======================================================================
1328
1329 template< typename T >
1330 void writeGoldBinaryValues( const _SubPartValues& subPartValues,
1331                             int                   nbValues,
1332                             int                   nbComponents,
1333                             int                   componentShift,
1334                             _BinaryFileWriter &   ensightDataFile )
1335 {
1336   const bool add3dComp = false;
1337   if ( subPartValues.myNumbers )
1338   {
1339     vector< _SelectedValueIterator<T> > compValIt =
1340       subPartValues.template getAllSelectedValues<T>( componentShift, nbComponents, add3dComp );
1341     ensightDataFile.addReal( compValIt, nbValues, MED_NO_INTERLACE );
1342   }
1343   else
1344   {
1345     vector< _ValueIterator<T> > compValIt =
1346       subPartValues.template getAllValues<T>( componentShift, nbComponents, add3dComp );
1347     ensightDataFile.addReal( compValIt, nbValues, MED_NO_INTERLACE );
1348   }
1349   // add the 3-d component to a vector in 2D space
1350   if ( nbComponents == 2 ) {
1351     vector<float> values( nbValues, 0. );
1352     ensightDataFile.addReal( & values[0], nbValues );
1353   }
1354 }
1355
1356 //=======================================================================
1357 //function : write
1358 //=======================================================================
1359
1360 void ENSIGHT_FIELD_WRONLY_DRIVER::write(void) const
1361   throw (MEDEXCEPTION)
1362 {
1363   const char * LOC = "ENSIGHT_FIELD_WRONLY_DRIVER::write(void) const " ;
1364   BEGIN_OF_MED(LOC);
1365
1366   openConst(false) ; // check if can write to case file
1367
1368   // Ensight case organization requires a main file (filename.case) which defines organization
1369
1370   _CaseFileDriver caseFile( getCaseFileName(), this);
1371   caseFile.read();
1372   caseFile.addField( this );
1373   caseFile.write(); // all fields of _CaseFileDriver_User are set by this method
1374
1375   openConst(true) ; // check if can write to data file
1376
1377   const FIELD_* field     = _ptrField;
1378   const SUPPORT * support = field->getSupport();
1379   const MESH * mesh       = support->getMesh();
1380
1381   int dt              = field->getIterationNumber();
1382   int it              = field->getOrderNumber();
1383   int nbComponents    = field->getNumberOfComponents();
1384   med_type_champ type = field->getValueType() ;
1385
1386   medEntityMesh entity = support->getEntity();
1387   int totalNbValues    = support->getNumberOfElements(MED_ALL_ELEMENTS);
1388   //const int* mainNbValsByType = support->getNumberOfElements();
1389
1390   int nbValuesByType = 0;
1391   int& componentShift = (type == MED_NO_INTERLACE_BY_TYPE) ? nbValuesByType : totalNbValues;
1392
1393   bool isOnAll = support->isOnAllElements();
1394   if (!isOnAll && !mesh)
1395     throw MED_EXCEPTION(STRING("Can't write field ") << field->getName() <<
1396                         ", which is not on all elements while mesh is not set to its support");
1397   if (!isOnAll)
1398     isOnAll = ( mesh->getNumberOfElements(entity,MED_ALL_ELEMENTS) == componentShift );
1399   if (!isOnAll && entity == MED_NODE && !isGoldFormat() ) {
1400     throw MED_EXCEPTION(compatibilityPb("Can't write field ") << field->getName() <<
1401                         " which is not on all nodes of the mesh in EnSight6 format,"
1402                         " use EnSight Gold format instead");
1403   }
1404   if ( entity == MED_EDGE && mesh && !isToWriteEntity( MED_EDGE, mesh ))
1405     throw MED_EXCEPTION(STRING(LOC) << "Can't write field " << field->getName()
1406                         << ", fields on MED_EDGE in volumic mesh are not supported");
1407
1408   // part number
1409   int partNum = getPartNumber( support );
1410   string isPartial = "";
1411   if ( !partNum ) {
1412     if ( !isGoldFormat() )
1413       throw MED_EXCEPTION(STRING("Can't write field ") << field->getName()
1414                           << " in EnSight6 format, it's support " << support->getName()
1415                           << " is not stored in geo file, use EnSight Gold format instead");
1416     isPartial = " partial";
1417     SUPPORT tmpSupport;
1418     tmpSupport.setAll(true);
1419     tmpSupport.setEntity( entity );
1420     partNum = getPartNumber( &tmpSupport );
1421   }
1422
1423   // supports to write the field for
1424   map< int, const SUPPORT* > parts;
1425   map< int, const SUPPORT* >::iterator partNumIt;
1426   parts[ partNum ] = support;
1427
1428   // In Ensight 6, nodal field is visible on all parts as nodes are global, for
1429   // the rest cases we write field for all groups in order to see values on them
1430   if ( isOnAll && ( isGoldFormat() || entity != MED_NODE ) && mesh )
1431   {
1432     // In Ensight Gold, to see nodal field on all groups we write them all
1433     int ent = entity, toEntity = entity + 1;
1434     if ( isGoldFormat() && entity == MED_NODE )
1435       ent = MED_CELL, toEntity = MED_ALL_ENTITIES;
1436
1437     for ( ; ent < toEntity; ++ent ) {
1438       int nbGroups = mesh->getNumberOfGroups(ent);
1439       for ( int i=1; i<=nbGroups; i++) {
1440         const GROUP* group = mesh->getGroup( ent, i );
1441         if (  group != support && !group->isOnAllElements() ) {
1442           partNum = getPartNumber( group );
1443           if ( partNum )
1444             parts.insert( make_pair( partNum, group ));
1445         }
1446       }
1447     }
1448   }
1449
1450   // description
1451   string description = field->getDescription();
1452   if ( description.empty() )
1453     description = STRING(field->getName()) << " dt=" << dt << " it=" << it;
1454
1455   cout << "-> Creating the Ensight data file " << getDataFileName() << endl ;
1456
1457
1458   _SubPartValues geoTypeValues;
1459   geoTypeValues.myOwnValues = false;
1460   geoTypeValues.myInterlace = field->getInterlacingType();
1461   int nbWrittenValues = 0;
1462
1463   if ( isBinaryEnSightFormatForWriting() )
1464   {
1465     // ======================================================
1466     //                          Binary
1467     // ======================================================
1468
1469     // function to write values
1470     typedef void (* TWriteValues) (const _SubPartValues& , int, int, int, _BinaryFileWriter &);
1471     TWriteValues writeValues;
1472     if ( isGoldFormat() ) {
1473       if ( type == MED_REEL64 )     writeValues = writeGoldBinaryValues<double>;
1474       else if ( type == MED_INT32 ) writeValues = writeGoldBinaryValues<int>;
1475       else                          writeValues = writeGoldBinaryValues<long>;
1476     }
1477     else {
1478       if ( type == MED_REEL64 )     writeValues = write6BinaryValues<double>;
1479       else if ( type == MED_INT32 ) writeValues = write6BinaryValues<int>;
1480       else                          writeValues = write6BinaryValues<long>;
1481     }
1482     _BinaryFileWriter ensightDataFile( getDataFileName() );
1483
1484     // description
1485     ensightDataFile.addString( description );
1486
1487     for ( partNumIt = parts.begin(); partNumIt != parts.end(); ++partNumIt)
1488     {
1489       const SUPPORT* partSup = partNumIt->second;
1490       partNum                = partNumIt->first;
1491       const bool otherEntity = ( entity != partSup->getEntity() ); // Gold, nodal field
1492       nbWrittenValues        = 0;
1493
1494       // part number
1495       if ( isGoldFormat() ) {
1496         ensightDataFile.addString("part");
1497         ensightDataFile.addInt( partNum );
1498       }
1499       else if ( entity != MED_NODE ) {
1500         ensightDataFile.addString( STRING("part ") << partNum );
1501       }
1502       // loop on types
1503       int                       nbTypes = partSup->getNumberOfTypes();
1504       const medGeometryElement* geoType = partSup->getTypes();
1505       for (int i=0; i<nbTypes; i++)
1506       {
1507         const medGeometryElement medType = geoType[i];
1508         nbValuesByType = partSup->getNumberOfElements(medType);
1509
1510         // type name
1511         if ( entity != MED_NODE ) {
1512           const TEnSightElemType& ensightType = getEnSightType(medType);
1513           ensightDataFile.addString( ensightType._name + isPartial );
1514         }
1515         else if ( isGoldFormat() ) {
1516           ensightDataFile.addString( "coordinates" + isPartial );
1517         }
1518
1519         // supporting numbers (Gold only)
1520         if ( !isPartial.empty() ) {
1521           const int *number = partSup->getNumber(medType);
1522           ensightDataFile.addInt( nbValuesByType );
1523           ensightDataFile.addInt( number, nbValuesByType );
1524         }
1525
1526         // get field values
1527         if ( otherEntity ) { // Gold, nodal field, non-nodal group (partSup)
1528           if ( isOnAll && partSup->isOnAllElements() ) {
1529             geoTypeValues.myNumbers = 0;
1530             geoTypeValues.myValues  = getValuePointer( 1, field );
1531             nbValuesByType = totalNbValues;
1532           }
1533           else {
1534             geoTypeValues.myNumbers = getNodeNumbers( partSup, nbValuesByType );
1535             geoTypeValues.myValues  = getValuePointer( geoTypeValues.myNumbers[0], field );
1536           }
1537         }
1538         else if ( partSup->getNumberOfElements(MED_ALL_ELEMENTS) != totalNbValues ) {
1539           geoTypeValues.myNumbers = partSup->getNumber(medType);
1540           geoTypeValues.myValues  = getValuePointer( geoTypeValues.myNumbers[0], field );
1541         }
1542         else {
1543           geoTypeValues.myNumbers = 0;
1544           geoTypeValues.myValues  = getValuePointer( nbWrittenValues + 1, field );
1545         }
1546         // write values
1547         writeValues(geoTypeValues, nbValuesByType, nbComponents, componentShift, ensightDataFile);
1548         nbWrittenValues += nbValuesByType;
1549
1550         if ( otherEntity ) {
1551           if ( geoTypeValues.myNumbers ) delete [] geoTypeValues.myNumbers;
1552           break; // all nodal values are written at once
1553         }
1554       }
1555     }
1556   }
1557   else
1558   {
1559     // ======================================================
1560     //                           ASCII
1561     // ======================================================
1562
1563     // function to write values
1564     typedef void (* TWriteValues) (const _SubPartValues& , int, int, int, ofstream &);
1565     TWriteValues writeValues;
1566     if ( isGoldFormat() ) {
1567       if ( type == MED_REEL64 )     writeValues = writeGoldASCIIValues<double>;
1568       else if ( type == MED_INT32 ) writeValues = writeGoldASCIIValues<int>;
1569       else                          writeValues = writeGoldASCIIValues<long>;
1570     }
1571     else {
1572       if ( type == MED_REEL64 )     writeValues = write6ASCIIValues<double>;
1573       else if ( type == MED_INT32 ) writeValues = write6ASCIIValues<int>;
1574       else                          writeValues = write6ASCIIValues<long>;
1575     }
1576     const int iw = isGoldFormat() ? 10 : 8; // width of integer
1577
1578     ofstream ensightDataFile( getDataFileName().c_str(), ios::out);
1579     ensightDataFile.setf(ios::scientific);
1580     ensightDataFile.precision(5);
1581
1582     // description
1583     ensightDataFile << description << endl;
1584
1585     for ( partNumIt = parts.begin(); partNumIt != parts.end(); ++partNumIt)
1586     {
1587       const SUPPORT* partSup = partNumIt->second;
1588       partNum                = partNumIt->first;
1589       const bool otherEntity = ( entity != partSup->getEntity() ); // Gold, nodal field
1590       nbWrittenValues        = 0;
1591
1592       // part number
1593       if ( isGoldFormat() )
1594         ensightDataFile << "part" << endl << setw(iw) << partNum << endl;
1595       else if ( entity != MED_NODE ) {
1596         ensightDataFile << "part " << partNum << endl;
1597       }
1598       // loop on types
1599       int                       nbTypes = partSup->getNumberOfTypes();
1600       const medGeometryElement* geoType = partSup->getTypes();
1601       for (int i=0; i<nbTypes; i++)
1602       {
1603         const medGeometryElement medType = geoType[i];
1604         nbValuesByType = partSup->getNumberOfElements(medType);
1605
1606         // type name
1607         if ( entity != MED_NODE ) {
1608           const TEnSightElemType& ensightType = getEnSightType(medType);
1609           ensightDataFile << ensightType._name << isPartial << endl;
1610         }
1611         else if ( isGoldFormat() ) {
1612           ensightDataFile << "coordinates" << isPartial << endl;
1613         }
1614
1615         // supporting numbers (Gold only)
1616         if ( !isPartial.empty() ) {
1617           const int *number = partSup->getNumber(medType);
1618           const int *nEnd   = number + nbValuesByType;
1619           ensightDataFile << setw(iw) << nbValuesByType << endl;
1620           while ( number < nEnd )
1621             ensightDataFile << setw(iw) << *number++ << endl;
1622         }
1623
1624         // get field values
1625         if ( otherEntity ) { // Gold, nodal field, non-nodal group (partSup)
1626           if ( isOnAll && partSup->isOnAllElements() ) {
1627             geoTypeValues.myNumbers = 0;
1628             geoTypeValues.myValues  = getValuePointer( 1, field );
1629             nbValuesByType = totalNbValues;
1630           }
1631           else {
1632             geoTypeValues.myNumbers = getNodeNumbers( partSup, nbValuesByType );
1633             geoTypeValues.myValues  = getValuePointer( geoTypeValues.myNumbers[0], field );
1634           }
1635         }
1636         else if ( partSup->getNumberOfElements(MED_ALL_ELEMENTS) != totalNbValues ) {
1637           geoTypeValues.myNumbers = partSup->getNumber(medType);
1638           geoTypeValues.myValues  = getValuePointer( geoTypeValues.myNumbers[0], field );
1639         }
1640         else {
1641           geoTypeValues.myNumbers = 0;
1642           geoTypeValues.myValues  = getValuePointer( nbWrittenValues + 1, field );
1643         }
1644
1645         // write values
1646         writeValues(geoTypeValues, nbValuesByType, nbComponents, componentShift, ensightDataFile);
1647         nbWrittenValues += nbValuesByType;
1648
1649         if ( otherEntity ) {
1650           if ( geoTypeValues.myNumbers ) delete [] geoTypeValues.myNumbers;
1651           break; // all nodal values are written at once
1652         }
1653       }
1654     } // loop on supports
1655
1656     ensightDataFile.close();
1657   }
1658 }