src/MEDMEM/MEDMEM_EnsightFieldDriver.cxx

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