salomeinclude_HEADERS= \
MEDLoaderDefines.hxx \
MEDLoader.hxx MEDLoaderBase.hxx MEDFileUtilities.hxx MEDFileMesh.hxx MEDFileMeshLL.hxx \
-MEDFileMeshElt.hxx MEDFileBasis.hxx MEDFileField.hxx MEDFileData.hxx
+MEDFileMeshElt.hxx MEDFileBasis.hxx MEDFileField.hxx MEDFileData.hxx \
+SauvMedConvertor.hxx SauvReader.hxx SauvWriter.hxx SauvUtilities.hxx
dist_libmedloader_la_SOURCES= \
MEDLoader.cxx MEDLoaderBase.cxx MEDFileUtilities.cxx \
MEDFileMesh.cxx MEDFileMeshElt.cxx MEDFileBasis.cxx \
-MEDFileMeshLL.cxx MEDFileField.cxx MEDFileData.cxx
+MEDFileMeshLL.cxx MEDFileField.cxx MEDFileData.cxx \
+SauvMedConvertor.cxx SauvReader.cxx SauvWriter.cxx
libmedloader_la_CPPFLAGS= $(MED2_INCLUDES) $(HDF5_INCLUDES) @CXXTMPDPTHFLAGS@ \
-I$(srcdir)/../INTERP_KERNEL \
--- /dev/null
+// Copyright (C) 2007-2011 CEA/DEN, EDF R&D
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+// File : SauvMedConvertor.cxx
+// Created : Tue Aug 16 14:43:20 2011
+// Author : Edward AGAPOV (eap)
+//
+
+#include "SauvMedConvertor.hxx"
+
+#include "CellModel.hxx"
+#include "MEDFileMesh.hxx"
+#include "MEDFileField.hxx"
+#include "MEDFileData.hxx"
+#include "MEDCouplingFieldDouble.hxx"
+
+#include <iostream>
+#include <cassert>
+#include <cmath>
+#include <queue>
+#include <limits>
+
+#include <cstdlib>
+#include <cstring>
+#include <fcntl.h>
+
+#ifdef WNT
+#include <io.h>
+#endif
+
+#ifndef WNT
+#define HAS_XDR
+#endif
+
+#ifdef HAS_XDR
+#include <rpc/xdr.h>
+#endif
+
+using namespace SauvUtilities;
+using namespace ParaMEDMEM;
+using namespace std;
+
+namespace
+{
+ // for ASCII file reader
+ const int GIBI_MaxOutputLen = 150; // max length of a line in the sauve file
+ const int GIBI_BufferSize = 16184; // for buffered reading
+
+ using namespace INTERP_KERNEL;
+
+ const size_t MaxMedCellType = NORM_ERROR;
+ const size_t NbGibiCellTypes = 47;
+ const TCellType GibiTypeToMed[NbGibiCellTypes] =
+ {
+ /*1 */ NORM_POINT1 ,/*2 */ NORM_SEG2 ,/*3 */ NORM_SEG3 ,/*4 */ NORM_TRI3 ,/*5 */ NORM_ERROR ,
+ /*6 */ NORM_TRI6 ,/*7 */ NORM_ERROR ,/*8 */ NORM_QUAD4 ,/*9 */ NORM_ERROR ,/*10*/ NORM_QUAD8 ,
+ /*11*/ NORM_ERROR ,/*12*/ NORM_ERROR ,/*13*/ NORM_ERROR ,/*14*/ NORM_HEXA8 ,/*15*/ NORM_HEXA20 ,
+ /*16*/ NORM_PENTA6 ,/*17*/ NORM_PENTA15,/*18*/ NORM_ERROR ,/*19*/ NORM_ERROR ,/*20*/ NORM_ERROR ,
+ /*21*/ NORM_ERROR ,/*22*/ NORM_ERROR ,/*23*/ NORM_TETRA4 ,/*24*/ NORM_TETRA10,/*25*/ NORM_PYRA5 ,
+ /*26*/ NORM_PYRA13 ,/*27*/ NORM_ERROR ,/*28*/ NORM_ERROR ,/*29*/ NORM_ERROR ,/*30*/ NORM_ERROR ,
+ /*31*/ NORM_ERROR ,/*32*/ NORM_ERROR ,/*33*/ NORM_ERROR ,/*34*/ NORM_ERROR ,/*35*/ NORM_ERROR ,
+ /*36*/ NORM_ERROR ,/*37*/ NORM_ERROR ,/*38*/ NORM_ERROR ,/*39*/ NORM_ERROR ,/*40*/ NORM_ERROR ,
+ /*41*/ NORM_ERROR ,/*42*/ NORM_ERROR ,/*43*/ NORM_ERROR ,/*44*/ NORM_ERROR ,/*45*/ NORM_ERROR ,
+ /*46*/ NORM_ERROR ,/*47*/ NORM_ERROR
+ };
+
+ //================================================================================
+ /*!
+ * \brief Return dimension of a group
+ */
+ //================================================================================
+
+ unsigned getDim( const Group* grp )
+ {
+ return SauvUtilities::getDimension( grp->_groups.empty() ? grp->_cellType : grp->_groups[0]->_cellType );
+ }
+
+ //================================================================================
+ /*!
+ * \brief Converts connectivity of quadratic elements
+ */
+ //================================================================================
+
+ inline void ConvertQuadratic( const INTERP_KERNEL::NormalizedCellType type,
+ const Cell & aCell )
+ {
+ if ( const int * conn = getGibi2MedQuadraticInterlace( type ))
+ {
+ Cell* ma = (Cell*) & aCell;
+ //cout << "###### BEFORE ConvertQuadratic() " << *ma << endl;
+ vector< Node* > new_nodes( ma->_nodes.size() );
+ for ( size_t i = 0; i < new_nodes.size(); ++i )
+ new_nodes[ i ] = ma->_nodes[ conn[ i ]];
+ ma->_nodes.swap( new_nodes );
+ //cout << "###### AFTER ConvertQuadratic() " << *ma << endl;
+ }
+ }
+
+ //================================================================================
+ /*!
+ * \brief Returns a vector of pairs of node indices to inverse a med volume element
+ */
+ //================================================================================
+
+ void getReverseVector (const INTERP_KERNEL::NormalizedCellType type,
+ vector<pair<int,int> > & swapVec )
+ {
+ swapVec.clear();
+
+ switch ( type )
+ {
+ case NORM_TETRA4:
+ swapVec.resize(1);
+ swapVec[0] = make_pair( 1, 2 );
+ break;
+ case NORM_PYRA5:
+ swapVec.resize(1);
+ swapVec[0] = make_pair( 1, 3 );
+ break;
+ case NORM_PENTA6:
+ swapVec.resize(2);
+ swapVec[0] = make_pair( 1, 2 );
+ swapVec[1] = make_pair( 4, 5 );
+ break;
+ case NORM_HEXA8:
+ swapVec.resize(2);
+ swapVec[0] = make_pair( 1, 3 );
+ swapVec[1] = make_pair( 5, 7 );
+ break;
+ case NORM_TETRA10:
+ swapVec.resize(3);
+ swapVec[0] = make_pair( 1, 2 );
+ swapVec[1] = make_pair( 4, 6 );
+ swapVec[2] = make_pair( 8, 9 );
+ break;
+ case NORM_PYRA13:
+ swapVec.resize(4);
+ swapVec[0] = make_pair( 1, 3 );
+ swapVec[1] = make_pair( 5, 8 );
+ swapVec[2] = make_pair( 6, 7 );
+ swapVec[3] = make_pair( 10, 12 );
+ break;
+ case NORM_PENTA15:
+ swapVec.resize(4);
+ swapVec[0] = make_pair( 1, 2 );
+ swapVec[1] = make_pair( 4, 5 );
+ swapVec[2] = make_pair( 6, 8 );
+ swapVec[3] = make_pair( 9, 11 );
+ break;
+ case NORM_HEXA20:
+ swapVec.resize(7);
+ swapVec[0] = make_pair( 1, 3 );
+ swapVec[1] = make_pair( 5, 7 );
+ swapVec[2] = make_pair( 8, 11 );
+ swapVec[3] = make_pair( 9, 10 );
+ swapVec[4] = make_pair( 12, 15 );
+ swapVec[5] = make_pair( 13, 14 );
+ swapVec[6] = make_pair( 17, 19 );
+ break;
+ // case NORM_SEG3: no need to reverse edges
+ // swapVec.resize(1);
+ // swapVec[0] = make_pair( 1, 2 );
+ // break;
+ case NORM_TRI6:
+ swapVec.resize(2);
+ swapVec[0] = make_pair( 1, 2 );
+ swapVec[1] = make_pair( 3, 5 );
+ break;
+ case NORM_QUAD8:
+ swapVec.resize(3);
+ swapVec[0] = make_pair( 1, 3 );
+ swapVec[1] = make_pair( 4, 7 );
+ swapVec[2] = make_pair( 5, 6 );
+ break;
+ default:;
+ }
+ }
+
+ //================================================================================
+ /*!
+ * \brief Inverses element orientation using vector of indices to swap
+ */
+ //================================================================================
+
+ inline void reverse(const Cell & aCell, const vector<pair<int,int> > & swapVec )
+ {
+ Cell* ma = (Cell*) & aCell;
+ for ( unsigned i = 0; i < swapVec.size(); ++i )
+ std::swap( ma->_nodes[ swapVec[i].first ],
+ ma->_nodes[ swapVec[i].second ]);
+ if ( swapVec.empty() )
+ ma->_reverse = true;
+ else
+ ma->_reverse = false;
+ }
+ //================================================================================
+ /*!
+ * \brief Comparator of cells by number used for ordering cells thinin a med group
+ */
+ struct TCellByIDCompare
+ {
+ bool operator () (const Cell* i1, const Cell* i2)
+ {
+ return i1->_number < i2->_number;
+ }
+ };
+ typedef map< const Cell*, unsigned, TCellByIDCompare > TCellToOrderMap;
+
+ //================================================================================
+ /*!
+ * \brief Fill Group::_relocTable if necessary
+ */
+ //================================================================================
+
+ void setRelocationTable( Group* grp, TCellToOrderMap& cell2order )
+ {
+ if ( !grp->_isProfile ) return;
+
+ // check if relocation table is necessary
+ bool isRelocated = false;
+ unsigned newOrder = 0;
+ TCellToOrderMap::iterator c2oIt = cell2order.begin(), c2oEnd = cell2order.end();
+ for ( ; !isRelocated && c2oIt != c2oEnd; ++c2oIt, ++newOrder )
+ isRelocated = ( c2oIt->second != newOrder );
+
+ if ( isRelocated )
+ {
+ grp->_relocTable.resize( cell2order.size() );
+ for ( newOrder = 0, c2oIt = cell2order.begin(); c2oIt != c2oEnd; ++c2oIt, ++newOrder )
+ grp->_relocTable[ c2oIt->second ] = newOrder;
+ }
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Return dimension for the given cell type
+ */
+//================================================================================
+
+unsigned SauvUtilities::getDimension( INTERP_KERNEL::NormalizedCellType type )
+{
+ return type == NORM_ERROR ? -1 : INTERP_KERNEL::CellModel::GetCellModel( type ).getDimension();
+}
+
+//================================================================================
+/*!
+ * \brief Returns interlace array to transform a quadratic GIBI element to a MED one
+ */
+//================================================================================
+
+const int * SauvUtilities::getGibi2MedQuadraticInterlace( INTERP_KERNEL::NormalizedCellType type )
+{
+ static vector<const int*> conn;
+ static const int hexa20 [] = {0,6,4,2, 12,18,16,14, 7,5,3,1, 19,17,15,13, 8,11,10,9};
+ static const int penta15[] = {0,2,4, 9,11,13, 1,3,5, 10,12,14, 6,7,3};
+ static const int pyra13 [] = {0,2,4,6, 12, 1,3,5,7, 8,9,10,11};
+ static const int tetra10[] = {0,2,4, 9, 1,3,5, 6,7,8};
+ static const int quad8 [] = {0,2,4,6, 1,3,5,7};
+ static const int tria6 [] = {0,2,4, 1,3,5};
+ static const int seg3 [] = {0,2,1};
+ if ( conn.empty() )
+ {
+ conn.resize( MaxMedCellType + 1, 0 );
+ conn[ NORM_HEXA20 ] = hexa20;
+ conn[ NORM_PENTA15] = penta15;
+ conn[ NORM_PYRA13 ] = pyra13;
+ conn[ NORM_TETRA10] = tetra10;
+ conn[ NORM_SEG3 ] = seg3;
+ conn[ NORM_TRI6 ] = tria6;
+ conn[ NORM_QUAD8 ] = quad8;
+ }
+ return conn[ type ];
+}
+
+//================================================================================
+/*!
+ * \brief Avoid coping sortedNodeIDs
+ */
+//================================================================================
+
+Cell::Cell(const Cell& ma)
+ : _nodes(ma._nodes), _reverse(ma._reverse), _sortedNodeIDs(0), _number(ma._number)
+{
+ if ( ma._sortedNodeIDs )
+ {
+ _sortedNodeIDs = new int[ _nodes.size() ];
+ std::copy( ma._sortedNodeIDs, ma._sortedNodeIDs + _nodes.size(), _sortedNodeIDs );
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Rerturn the i-th link of face
+ */
+//================================================================================
+
+SauvUtilities::Link Cell::link(int i) const
+{
+ int i2 = ( i + 1 ) % _nodes.size();
+ if ( _reverse )
+ return make_pair( _nodes[i2]->_number, _nodes[i]->_number );
+ else
+ return make_pair( _nodes[i]->_number, _nodes[i2]->_number );
+}
+
+//================================================================================
+/*!
+ * \brief creates if needed and return _sortedNodeIDs
+ */
+//================================================================================
+
+const TID* Cell::getSortedNodes() const
+{
+ if ( !_sortedNodeIDs )
+ {
+ size_t l=_nodes.size();
+ _sortedNodeIDs = new int[ l ];
+
+ for (size_t i=0; i!=l; ++i)
+ _sortedNodeIDs[i]=_nodes[i]->_number;
+ std::sort( _sortedNodeIDs, _sortedNodeIDs + l );
+ }
+ return _sortedNodeIDs;
+}
+
+//================================================================================
+/*!
+ * \brief Compare sorted ids of cell nodes
+ */
+//================================================================================
+
+bool Cell::operator< (const Cell& ma) const
+{
+ if ( _nodes.size() == 1 )
+ return _nodes[0] < ma._nodes[0];
+
+ const int* v1 = getSortedNodes();
+ const int* v2 = ma.getSortedNodes();
+ for ( const int* vEnd = v1 + _nodes.size(); v1 < vEnd; ++v1, ++v2 )
+ if(*v1 != *v2)
+ return *v1 < *v2;
+ return false;
+}
+
+//================================================================================
+/*!
+ * \brief Dump a Cell
+ */
+//================================================================================
+
+std::ostream& SauvUtilities::operator<< (std::ostream& os, const SauvUtilities::Cell& ma)
+{
+ os << "cell " << ma._number << " (" << ma._nodes.size() << " nodes) : < " << ma._nodes[0];
+ for( size_t i=1; i!=ma._nodes.size(); ++i)
+ os << ", " << ma._nodes[0];
+#ifdef _DEBUG_
+ os << " > sortedNodes: ";
+ if ( ma._sortedNodeIDs ) {
+ os << "< ";
+ for( size_t i=0; i!=ma._nodes.size(); ++i)
+ os << ( i ? ", " : "" ) << ma._sortedNodeIDs[i];
+ os << " >";
+ }
+ else {
+ os << "NULL";
+ }
+#endif
+ return os;
+}
+
+//================================================================================
+/*!
+ * \brief Return nb of elements in the group
+ */
+//================================================================================
+
+int Group::size() const
+{
+ int size = 0;
+ if ( !_relocTable.empty() )
+ size = _relocTable.size();
+ else if ( _medGroup )
+ size = _medGroup->getNumberOfTuples();
+ else if ( !_cells.empty() )
+ size = _cells.size();
+ else
+ for ( size_t i = 0; i < _groups.size(); ++i )
+ size += _groups[i]->size();
+ return size;
+}
+
+//================================================================================
+/*!
+ * \brief Conver gibi element type to med one
+ */
+//================================================================================
+
+INTERP_KERNEL::NormalizedCellType SauvUtilities::gibi2medGeom( size_t gibiType )
+{
+ if ( gibiType < 1 || gibiType > NbGibiCellTypes )
+ return NORM_ERROR;
+
+ return GibiTypeToMed[ gibiType - 1 ];
+}
+
+//================================================================================
+/*!
+ * \brief Conver med element type to gibi one
+ */
+//================================================================================
+
+int SauvUtilities::med2gibiGeom( INTERP_KERNEL::NormalizedCellType medGeomType )
+{
+ for ( int i = 0; i < NbGibiCellTypes; i++ )
+ if ( GibiTypeToMed[ i ] == medGeomType )
+ return i + 1;
+
+ return -1;
+}
+
+//================================================================================
+/*!
+ * \brief Remember the file name
+ */
+//================================================================================
+
+FileReader::FileReader(const char* fileName):_fileName(fileName),_iRead(0),_nbToRead(0)
+{
+}
+
+//================================================================================
+/*!
+ * \brief Constructor of ASCII sauve file reader
+ */
+//================================================================================
+
+ASCIIReader::ASCIIReader(const char* fileName)
+ :FileReader(fileName),
+ _file(-1)
+{
+}
+
+//================================================================================
+/*!
+ * \brief Return true
+ */
+//================================================================================
+
+bool ASCIIReader::isASCII() const
+{
+ return true;
+}
+
+//================================================================================
+/*!
+ * \brief Try to open an ASCII file
+ */
+//================================================================================
+
+bool ASCIIReader::open()
+{
+#ifdef WNT
+ _file = ::_open (_fileName.c_str(), _O_RDONLY|_O_BINARY);
+#else
+ _file = ::open (_fileName.c_str(), O_RDONLY);
+#endif
+ if (_file >= 0)
+ {
+ _start = new char [GIBI_BufferSize]; // working buffer beginning
+ //_tmpBuf = new char [GIBI_MaxOutputLen];
+ _ptr = _start;
+ _eptr = _start;
+ _lineNb = 0;
+ }
+ else
+ {
+ //THROW_IK_EXCEPTION("Can't open file "<<_fileName << " fd: " << _file);
+ }
+ return (_file >= 0);
+}
+
+//================================================================================
+/*!
+ * \brief Close the file
+ */
+//================================================================================
+
+ASCIIReader::~ASCIIReader()
+{
+ if (_file >= 0)
+ {
+ ::close (_file);
+ if (_start != 0L)
+ {
+ delete [] _start;
+ //delete [] _tmpBuf;
+ _start = 0;
+ }
+ _file = -1;
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Return a next line of the file
+ */
+//================================================================================
+
+bool ASCIIReader::getNextLine (char* & line, bool raiseOEF /*= true*/ )
+{
+ if ( getLine( line )) return true;
+ if ( raiseOEF )
+ THROW_IK_EXCEPTION("Unexpected EOF on ln "<<_lineNb);
+ return false;
+}
+
+//================================================================================
+/*!
+ * \brief Read a next line of the file if necessary
+ */
+//================================================================================
+
+bool ASCIIReader::getLine(char* & line)
+{
+ bool aResult = true;
+ // Check the state of the buffer;
+ // if there is too little left, read the next portion of data
+ int nBytesRest = _eptr - _ptr;
+ if (nBytesRest < GIBI_MaxOutputLen)
+ {
+ if (nBytesRest > 0)
+ {
+ // move the remaining portion to the buffer beginning
+ for ( int i = 0; i < nBytesRest; ++i )
+ _start[i] = _ptr[i];
+ //memcpy (_tmpBuf, _ptr, nBytesRest);
+ //memcpy (_start, _tmpBuf, nBytesRest);
+ }
+ else
+ {
+ nBytesRest = 0;
+ }
+ _ptr = _start;
+ const int nBytesRead = ::read (_file,
+ &_start [nBytesRest],
+ GIBI_BufferSize - nBytesRest);
+ nBytesRest += nBytesRead;
+ _eptr = &_start [nBytesRest];
+ }
+ // Check the buffer for the end-of-line
+ char * ptr = _ptr;
+ while (true)
+ {
+ // Check for end-of-the-buffer, the ultimate criterion for termination
+ if (ptr >= _eptr)
+ {
+ if (nBytesRest <= 0)
+ aResult = false;
+ else
+ _eptr[-1] = '\0';
+ break;
+ }
+ // seek the line-feed character
+ if (ptr[0] == '\n')
+ {
+ if (ptr[-1] == '\r')
+ ptr[-1] = '\0';
+ ptr[0] = '\0';
+ ++ptr;
+ break;
+ }
+ ++ptr;
+ }
+ // Output the result
+ line = _ptr;
+ _ptr = ptr;
+ _lineNb++;
+
+ return aResult;
+}
+
+//================================================================================
+/*!
+ * \brief Prepare for iterating over given nb of values
+ * \param nbToRead - nb of fields to read
+ * \param nbPosInLine - nb of fields in one line
+ * \param width - field width
+ * \param shift - shift from line beginning to the field start
+ */
+//================================================================================
+
+void ASCIIReader::init( int nbToRead, int nbPosInLine, int width, int shift /*= 0*/ )
+{
+ _nbToRead = nbToRead;
+ _nbPosInLine = nbPosInLine;
+ _width = width;
+ _shift = shift;
+ _iPos = _iRead = 0;
+ if ( _nbToRead )
+ {
+ getNextLine( _curPos );
+ _curPos = _curPos + _shift;
+ }
+ else
+ {
+ _curPos = 0;
+ }
+ _curLocale.clear();
+}
+
+//================================================================================
+/*!
+ * \brief Prepare for iterating over given nb of string values
+ */
+//================================================================================
+
+void ASCIIReader::initNameReading(int nbValues, int width /*= 8*/)
+{
+ init( nbValues, 72 / ( width + 1 ), width, 1 );
+}
+
+//================================================================================
+/*!
+ * \brief Prepare for iterating over given nb of integer values
+ */
+//================================================================================
+
+void ASCIIReader::initIntReading(int nbValues)
+{
+ init( nbValues, 10, 8 );
+}
+
+//================================================================================
+/*!
+ * \brief Prepare for iterating over given nb of real values
+ */
+//================================================================================
+
+void ASCIIReader::initDoubleReading(int nbValues)
+{
+ init( nbValues, 3, 22 );
+
+ // Correction 2 of getDouble(): set "C" numeric locale to read numbers
+ // with dot decimal point separator, as it is in SAUVE files
+ _curLocale = setlocale(LC_NUMERIC, "C");
+}
+
+//================================================================================
+/*!
+ * \brief Return true if not all values have been read
+ */
+//================================================================================
+
+bool ASCIIReader::more() const
+{
+ bool result = false;
+ if ( _iRead < _nbToRead)
+ {
+ if ( _curPos ) result = true;
+ }
+ return result;
+}
+
+//================================================================================
+/*!
+ * \brief Go to the nex value
+ */
+//================================================================================
+
+void ASCIIReader::next()
+{
+ if ( !more() )
+ THROW_IK_EXCEPTION("SauvUtilities::ASCIIReader::next(): no more() values to read");
+ ++_iRead;
+ ++_iPos;
+ if ( _iRead < _nbToRead )
+ {
+ if ( _iPos >= _nbPosInLine )
+ {
+ getNextLine( _curPos );
+ _curPos = _curPos + _shift;
+ _iPos = 0;
+ }
+ else
+ {
+ _curPos = _curPos + _width + _shift;
+ }
+ }
+ else
+ {
+ _curPos = 0;
+ if ( !_curLocale.empty() )
+ {
+ setlocale(LC_NUMERIC, _curLocale.c_str());
+ _curLocale.clear();
+ }
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Return the current integer value
+ */
+//================================================================================
+
+int ASCIIReader::getInt() const
+{
+ // fix for two glued ints (issue 0021009):
+ // Line nb | File contents
+ // ------------------------------------------------------------------------------------
+ // 53619905 | 1 2 6 8
+ // 53619906 | SCALAIRE
+ // 53619907 | -63312600499 1 0 0 0 -2 0 2
+ // where -63312600499 is actualy -633 and 12600499
+ char hold=_curPos[_width];
+ _curPos[_width] = '\0';
+ int result = atoi( _curPos );
+ _curPos[_width] = hold;
+ return result;
+ //return atoi(str());
+}
+
+//================================================================================
+/*!
+ * \brief Return the current float value
+ */
+//================================================================================
+
+float ASCIIReader::getFloat() const
+{
+ return getDouble();
+}
+
+//================================================================================
+/*!
+ * \brief Return the current double value
+ */
+//================================================================================
+
+double ASCIIReader::getDouble() const
+{
+ //std::string aStr (_curPos);
+
+ // Correction: add missing 'E' specifier
+ // int aPosStart = aStr.find_first_not_of(" \t");
+ // if (aPosStart < (int)aStr.length()) {
+ // int aPosSign = aStr.find_first_of("+-", aPosStart + 1); // pass the leading symbol, as it can be a sign
+ // if (aPosSign < (int)aStr.length()) {
+ // if (aStr[aPosSign - 1] != 'e' && aStr[aPosSign - 1] != 'E')
+ // aStr.insert(aPosSign, "E", 1);
+ // }
+ // }
+
+ // Different Correction (more optimal)
+ // Sample:
+ // 0.00000000000000E+00 -2.37822406690632E+01 6.03062748797469E+01
+ // 7.70000000000000-100 7.70000000000000+100 7.70000000000000+100
+ //0123456789012345678901234567890123456789012345678901234567890123456789
+ const size_t posE = 18;
+ if ( _curPos[posE] != 'E' && _curPos[posE] != 'e' )
+ {
+ std::string aStr (_curPos);
+ if ( aStr.size() < posE )
+ THROW_IK_EXCEPTION("No more doubles (line #" << lineNb() << ")");
+ aStr.insert( posE, "E", 1 );
+ return atof(aStr.c_str());
+ }
+ return atof( _curPos );
+}
+
+//================================================================================
+/*!
+ * \brief Return the current string value
+ */
+//================================================================================
+
+string ASCIIReader::getName() const
+{
+ int len = _width;
+ while (( _curPos[len-1] == ' ' || _curPos[len-1] == 0) && len > 0 )
+ len--;
+ return string( _curPos, len );
+}
+
+//================================================================================
+/*!
+ * \brief Constructor of a binary sauve file reader
+ */
+//================================================================================
+
+XDRReader::XDRReader(const char* fileName) :FileReader(fileName), _xdrs_file(NULL)
+{
+}
+
+//================================================================================
+/*!
+ * \brief Close the XDR sauve file
+ */
+//================================================================================
+
+XDRReader::~XDRReader()
+{
+#ifdef HAS_XDR
+ if ( _xdrs_file )
+ {
+ xdr_destroy((XDR*)_xdrs);
+ free((XDR*)_xdrs);
+ ::fclose(_xdrs_file);
+ _xdrs_file = NULL;
+ }
+#endif
+}
+
+//================================================================================
+/*!
+ * \brief Return false
+ */
+//================================================================================
+
+bool XDRReader::isASCII() const
+{
+ return false;
+}
+
+//================================================================================
+/*!
+ * \brief Try to open an XRD file
+ */
+//================================================================================
+
+bool XDRReader::open()
+{
+ bool xdr_ok = false;
+#ifdef HAS_XDR
+ if ((_xdrs_file = ::fopen(_fileName.c_str(), "r")))
+ {
+ _xdrs = (XDR *)malloc(sizeof(XDR));
+ xdrstdio_create((XDR*)_xdrs, _xdrs_file, XDR_DECODE);
+
+ const int maxsize = 10;
+ char icha[maxsize+1];
+ char* icha2 = icha;
+ if (( xdr_ok = xdr_string((XDR*)_xdrs, &icha2, maxsize)))
+ {
+ icha[maxsize] = '\0';
+ xdr_ok = (strcmp(icha, "CASTEM XDR") == 0);
+ }
+ if ( !xdr_ok )
+ {
+ xdr_destroy((XDR*)_xdrs);
+ free((XDR*)_xdrs);
+ fclose(_xdrs_file);
+ _xdrs_file = NULL;
+ }
+ }
+#endif
+ return xdr_ok;
+}
+
+//================================================================================
+/*!
+ * \brief A stub
+ */
+//================================================================================
+
+bool XDRReader::getNextLine (char* &, bool )
+{
+ return true;
+}
+
+//================================================================================
+/*!
+ * \brief Prepare for iterating over given nb of values
+ * \param nbToRead - nb of fields to read
+ * \param width - field width
+ */
+//================================================================================
+
+void XDRReader::init( int nbToRead, int width/*=0*/ )
+{
+ if(_iRead < _nbToRead)
+ {
+ cout << "_iRead, _nbToRead : " << _iRead << " " << _nbToRead << endl;
+ cout << "Unfinished iteration before new one !" << endl;
+ THROW_IK_EXCEPTION("SauvUtilities::XDRReader::init(): Unfinished iteration before new one !");
+ }
+ _iRead = 0;
+ _nbToRead = nbToRead;
+ _width = width;
+}
+
+//================================================================================
+/*!
+ * \brief Prepare for iterating over given nb of string values
+ */
+//================================================================================
+
+void XDRReader::initNameReading(int nbValues, int width)
+{
+ init( nbValues, width );
+ _xdr_kind = _xdr_kind_char;
+ if(nbValues*width)
+ {
+ unsigned int nels = nbValues*width;
+ _xdr_cvals = (char*)malloc((nels+1)*sizeof(char));
+#ifdef HAS_XDR
+ xdr_string((XDR*)_xdrs, &_xdr_cvals, nels);
+#endif
+ _xdr_cvals[nels] = '\0';
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Prepare for iterating over given nb of integer values
+ */
+//================================================================================
+
+void XDRReader::initIntReading(int nbValues)
+{
+ init( nbValues );
+ _xdr_kind = _xdr_kind_int;
+ if(nbValues)
+ {
+#ifdef HAS_XDR
+ unsigned int nels = nbValues;
+ unsigned int actual_nels;
+ _xdr_ivals = (int*)malloc(nels*sizeof(int));
+ xdr_array((XDR*)_xdrs, (char **)&_xdr_ivals, &actual_nels, nels, sizeof(int), (xdrproc_t)xdr_int);
+#endif
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Prepare for iterating over given nb of real values
+ */
+//================================================================================
+
+void XDRReader::initDoubleReading(int nbValues)
+{
+ init( nbValues );
+ _xdr_kind = _xdr_kind_double;
+ if(nbValues)
+ {
+#ifdef HAS_XDR
+ unsigned int nels = nbValues;
+ unsigned int actual_nels;
+ _xdr_dvals = (double*)malloc(nels*sizeof(double));
+ xdr_array((XDR*)_xdrs, (char **)&_xdr_dvals, &actual_nels, nels, sizeof(double), (xdrproc_t)xdr_double);
+#endif
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Return true if not all values have been read
+ */
+//================================================================================
+
+bool XDRReader::more() const
+{
+ return _iRead < _nbToRead;
+}
+
+//================================================================================
+/*!
+ * \brief Go to the nex value
+ */
+//================================================================================
+
+void XDRReader::next()
+{
+ if ( !more() )
+ THROW_IK_EXCEPTION("SauvUtilities::XDRReader::next(): no more() values to read");
+
+ ++_iRead;
+ if ( _iRead < _nbToRead )
+ {
+ }
+ else
+ {
+ if(_xdr_kind == _xdr_kind_char) free(_xdr_cvals);
+ if(_xdr_kind == _xdr_kind_int) free(_xdr_ivals);
+ if(_xdr_kind == _xdr_kind_double) free(_xdr_dvals);
+ _xdr_kind = _xdr_kind_null;
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Return the current integer value
+ */
+//================================================================================
+
+int XDRReader::getInt() const
+{
+ if(_iRead < _nbToRead)
+ {
+ return _xdr_ivals[_iRead];
+ }
+ else
+ {
+ int result = 0;
+#ifdef HAS_XDR
+ xdr_int((XDR*)_xdrs, &result);
+#endif
+ return result;
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Return the current float value
+ */
+//================================================================================
+
+float XDRReader::getFloat() const
+{
+ float result = 0;
+#ifdef HAS_XDR
+ xdr_float((XDR*)_xdrs, &result);
+#endif
+ return result;
+}
+
+//================================================================================
+/*!
+ * \brief Return the current double value
+ */
+//================================================================================
+
+double XDRReader::getDouble() const
+{
+ if(_iRead < _nbToRead)
+ {
+ return _xdr_dvals[_iRead];
+ }
+ else
+ {
+ double result = 0;
+#ifdef HAS_XDR
+ xdr_double((XDR*)_xdrs, &result);
+#endif
+ return result;
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Return the current string value
+ */
+//================================================================================
+
+std::string XDRReader::getName() const
+{
+ int len = _width;
+ char* s = _xdr_cvals + _iRead*_width;
+ while (( s[len-1] == ' ' || s[len-1] == 0) && len > 0 )
+ len--;
+ return string( s, len );
+}
+
+//================================================================================
+/*!
+ * \brief Throw an exception if not all needed data is present
+ */
+//================================================================================
+
+void IntermediateMED::checkDataAvailability() const throw(INTERP_KERNEL::Exception)
+{
+ if ( _spaceDim == 0 )
+ THROW_IK_EXCEPTION("Wrong file format"); // it is the first record in the sauve file
+
+ if ( _groups.empty() )
+ THROW_IK_EXCEPTION("No elements have been read");
+
+ if ( _points.empty() || _nbNodes == 0 )
+ THROW_IK_EXCEPTION("Nodes of elements are not filled");
+
+ if ( _coords.empty() )
+ THROW_IK_EXCEPTION("Node coordinates are missing");
+
+ if ( _coords.size() < _nbNodes * _spaceDim )
+ THROW_IK_EXCEPTION("Nodes and coordinates mismatch");
+}
+
+//================================================================================
+/*!
+ * \brief Makes ParaMEDMEM::MEDFileData from self
+ */
+//================================================================================
+
+ParaMEDMEM::MEDFileData* IntermediateMED::convertInMEDFileDS()
+{
+ MEDCouplingAutoRefCountObjectPtr< MEDFileUMesh > mesh = makeMEDFileMesh();
+ MEDCouplingAutoRefCountObjectPtr< MEDFileFields > fields = makeMEDFileFields(mesh);
+
+ MEDCouplingAutoRefCountObjectPtr< MEDFileMeshes > meshes = MEDFileMeshes::New();
+ MEDCouplingAutoRefCountObjectPtr< MEDFileData > medData = MEDFileData::New();
+ meshes->pushMesh( mesh );
+ medData->setMeshes( meshes );
+ if ( fields ) medData->setFields( fields );
+
+ medData->incrRef();
+ return medData;
+}
+
+//================================================================================
+/*!
+ * \brief Creates ParaMEDMEM::MEDFileUMesh from its data
+ */
+//================================================================================
+
+ParaMEDMEM::MEDFileUMesh* IntermediateMED::makeMEDFileMesh()
+{
+ // check if all needed piles are present
+ checkDataAvailability();
+
+ // set long names
+ setGroupLongNames();
+
+ // fix element orientation
+ if ( _spaceDim == 2 )
+ orientElements2D();
+ else if ( _spaceDim == 3 )
+ orientElements3D();
+
+ // process groups
+ decreaseHierarchicalDepthOfSubgroups();
+ eraseUselessGroups();
+ detectMixDimGroups();
+
+ // assign IDs
+ _points.numberNodes();
+ numberElements();
+
+ // make the med mesh
+
+ MEDFileUMesh* mesh = MEDFileUMesh::New();
+
+ DataArrayDouble *coords = getCoords();
+ setConnectivity( mesh, coords );
+ setGroups( mesh );
+
+ coords->decrRef();
+
+ if ( !mesh->getName() || strlen( mesh->getName() ) == 0 )
+ mesh->setName( "MESH" );
+
+ return mesh;
+}
+
+//================================================================================
+/*!
+ * \brief Set long names to groups
+ */
+//================================================================================
+
+void IntermediateMED::setGroupLongNames()
+{
+ // IMP 0020434: mapping GIBI names to MED names
+ // set med names to objects (mesh, fields, support, group or other)
+
+ set<int> treatedGroups;
+
+ list<nameGIBItoMED>::iterator itGIBItoMED = _listGIBItoMED_mail.begin();
+ for (; itGIBItoMED != _listGIBItoMED_mail.end(); itGIBItoMED++)
+ {
+ if ( (int)_groups.size() < itGIBItoMED->gibi_id ) continue;
+
+ SauvUtilities::Group & grp = _groups[itGIBItoMED->gibi_id - 1];
+
+ // if there are several names for grp then the 1st name is the name
+ // of grp and the rest ones are names of groups referring grp (issue 0021311)
+ const bool isRefName = !treatedGroups.insert( itGIBItoMED->gibi_id ).second;
+ if ( !isRefName )
+ grp._name = _mapStrings[ itGIBItoMED->med_id ];
+ else
+ for ( unsigned i = 0; i < grp._refNames.size(); ++i )
+ if ( grp._refNames[i].empty() )
+ grp._refNames[i] = _mapStrings[ (*itGIBItoMED).med_id ];
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Set long names to fields
+ */
+//================================================================================
+
+void IntermediateMED::setFieldLongNames(set< string >& usedNames)
+{
+ list<nameGIBItoMED>::iterator itGIBItoMED = _listGIBItoMED_cham.begin();
+ for (; itGIBItoMED != _listGIBItoMED_cham.end(); itGIBItoMED++)
+ {
+ if (itGIBItoMED->gibi_pile == PILE_FIELD)
+ {
+ _cellFields[itGIBItoMED->gibi_id - 1]->_name = _mapStrings[itGIBItoMED->med_id];
+ }
+ else if (itGIBItoMED->gibi_pile == PILE_NODES_FIELD)
+ {
+ _nodeFields[itGIBItoMED->gibi_id - 1]->_name = _mapStrings[itGIBItoMED->med_id];
+ }
+ } // iterate on _listGIBItoMED_cham
+
+ for (itGIBItoMED =_listGIBItoMED_comp.begin(); itGIBItoMED != _listGIBItoMED_comp.end(); itGIBItoMED++)
+ {
+ string medName = _mapStrings[itGIBItoMED->med_id];
+ string gibiName = _mapStrings[itGIBItoMED->gibi_id];
+
+ bool name_found = false;
+ for ( int isNodal = 0; isNodal < 2 && !name_found; ++isNodal )
+ {
+ vector<DoubleField* > & fields = isNodal ? _nodeFields : _cellFields;
+ for ( size_t ifi = 0; ifi < fields.size() && !name_found; ifi++)
+ {
+ if (medName.find( fields[ifi]->_name + "." ) == 0 )
+ {
+ vector<DoubleField::_Sub_data>& aSubDs = fields[ifi]->_sub;
+ int nbSub = aSubDs.size();
+ for (int isu = 0; isu < nbSub; isu++)
+ for (int ico = 0; ico < aSubDs[isu].nbComponents(); ico++)
+ {
+ if (aSubDs[isu].compName(ico) == gibiName)
+ {
+ string medNameCompo = medName.substr( fields[ifi]->_name.size() + 1 );
+ fields[ifi]->_sub[isu].compName(ico) = medNameCompo;
+ }
+ }
+ }
+ }
+ }
+ } // iterate on _listGIBItoMED_comp
+
+ for ( size_t i = 0; i < _nodeFields.size() ; i++)
+ usedNames.insert( _nodeFields[i]->_name );
+ for ( size_t i = 0; i < _cellFields.size() ; i++)
+ usedNames.insert( _cellFields[i]->_name );
+}
+
+//================================================================================
+/*!
+ * \brief Decrease hierarchical depth of subgroups
+ */
+//================================================================================
+
+void IntermediateMED::decreaseHierarchicalDepthOfSubgroups()
+{
+ for (size_t i=0; i!=_groups.size(); ++i)
+ {
+ Group& grp = _groups[i];
+ for (size_t j = 0; j < grp._groups.size(); ++j )
+ {
+ Group & sub_grp = *grp._groups[j];
+ if ( !sub_grp._groups.empty() )
+ {
+ // replace j with its 1st subgroup
+ grp._groups[j] = sub_grp._groups[0];
+ // push back the rest subs
+ grp._groups.insert( grp._groups.end(), ++sub_grp._groups.begin(), sub_grp._groups.end() );
+ }
+ }
+ // remove empty sub-_groups
+ vector< Group* > newSubGroups;
+ newSubGroups.reserve( grp._groups.size() );
+ for (size_t j = 0; j < grp._groups.size(); ++j )
+ if ( !grp._groups[j]->empty() )
+ newSubGroups.push_back( grp._groups[j] );
+ if ( newSubGroups.size() < grp._groups.size() )
+ grp._groups.swap( newSubGroups );
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Erase _groups that won't be converted
+ */
+//================================================================================
+
+void IntermediateMED::eraseUselessGroups()
+{
+ // propagate _isProfile=true to sub-groups of composite groups
+ // for (size_t int i=0; i!=_groups.size(); ++i)
+ // {
+ // Group* grp = _groups[i];
+ // if ( grp->_isProfile && !grp->_groups.empty() )
+ // for (size_t j = 0; j < grp->_groups.size(); ++j )
+ // grp->_groups[j]->_isProfile=true;
+ // }
+ std::set<Group*> groups2convert;
+ // keep not named sub-groups of field supports
+ for (size_t i=0; i!=_groups.size(); ++i)
+ {
+ Group& grp = _groups[i];
+ if ( grp._isProfile && !grp._groups.empty() )
+ groups2convert.insert( grp._groups.begin(), grp._groups.end() );
+ }
+
+ // keep named groups and their subgroups
+ for (size_t i=0; i!=_groups.size(); ++i)
+ {
+ Group& grp = _groups[i];
+ if ( !grp._name.empty() && !grp.empty() )
+ {
+ groups2convert.insert( &grp );
+ groups2convert.insert( grp._groups.begin(), grp._groups.end() );
+ }
+ }
+ // erase groups that are not in groups2convert and not _isProfile
+ for (size_t i=0; i!=_groups.size(); ++i)
+ {
+ Group* grp = &_groups[i];
+ if ( !grp->_isProfile && !groups2convert.count( grp ) )
+ {
+ grp->_cells.clear();
+ grp->_groups.clear();
+ }
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Detect _groups of mixed dimension
+ */
+//================================================================================
+
+void IntermediateMED::detectMixDimGroups()
+{
+ //hasMixedCells = false;
+ for ( size_t i=0; i < _groups.size(); ++i )
+ {
+ Group& grp = _groups[i];
+ if ( grp._groups.size() < 2 )
+ continue;
+
+ // check if sub-groups have different dimension
+ unsigned dim1 = getDim( &grp );
+ for ( size_t j = 1; j < grp._groups.size(); ++j )
+ {
+ unsigned dim2 = getDim( &grp );
+ if ( dim1 != dim2 )
+ {
+ grp._cells.clear();
+ grp._groups.clear();
+ if ( !grp._name.empty() )
+ cout << "Erase a group with elements of different dim |" << grp._name << "|"<< endl;
+ break;
+ }
+ }
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Fix connectivity of elements in 2D space
+ */
+//================================================================================
+
+void IntermediateMED::orientElements2D()
+{
+ set<Cell>::const_iterator elemIt, elemEnd;
+ vector< pair<int,int> > swapVec;
+
+ // ------------------------------------
+ // fix connectivity of quadratic edges
+ // ------------------------------------
+ set<Cell>& quadEdges = _cellsByType[ INTERP_KERNEL::NORM_SEG3 ];
+ if ( !quadEdges.empty() )
+ {
+ elemIt = quadEdges.begin(), elemEnd = quadEdges.end();
+ for ( ; elemIt != elemEnd; ++elemIt )
+ ConvertQuadratic( INTERP_KERNEL::NORM_SEG3, *elemIt );
+ }
+
+ CellsByDimIterator faceIt( *this, 2 );
+ while ( const set<Cell > * faces = faceIt.nextType() )
+ {
+ TCellType cellType = faceIt.type();
+ bool isQuadratic = getGibi2MedQuadraticInterlace( cellType );
+
+ getReverseVector( cellType, swapVec );
+
+ // ------------------------------------
+ // fix connectivity of quadratic faces
+ // ------------------------------------
+ if ( isQuadratic )
+ for ( elemIt = faces->begin(), elemEnd = faces->end(); elemIt != elemEnd; elemIt++ )
+ ConvertQuadratic( cellType, *elemIt );
+
+ // --------------------------
+ // orient faces clockwise
+ // --------------------------
+ int iQuad = isQuadratic ? 2 : 1;
+ for ( elemIt = faces->begin(), elemEnd = faces->end(); elemIt != elemEnd; elemIt++ )
+ {
+ // look for index of the most left node
+ int iLeft = 0, iNode, nbNodes = elemIt->_nodes.size() / iQuad;
+ double x, minX = nodeCoords( elemIt->_nodes[0] )[0];
+ for ( iNode = 1; iNode < nbNodes; ++iNode )
+ if (( x = nodeCoords( elemIt->_nodes[ iNode ])[ 0 ]) < minX )
+ minX = x, iLeft = iNode;
+
+ // indeces of the nodes neighboring the most left one
+ int iPrev = ( iLeft - 1 < 0 ) ? nbNodes - 1 : iLeft - 1;
+ int iNext = ( iLeft + 1 == nbNodes ) ? 0 : iLeft + 1;
+ // find components of prev-left and left-next vectors
+ double xP = nodeCoords( elemIt->_nodes[ iPrev ])[ 0 ];
+ double yP = nodeCoords( elemIt->_nodes[ iPrev ])[ 1 ];
+ double xN = nodeCoords( elemIt->_nodes[ iNext ])[ 0 ];
+ double yN = nodeCoords( elemIt->_nodes[ iNext ])[ 1 ];
+ double xL = nodeCoords( elemIt->_nodes[ iLeft ])[ 0 ];
+ double yL = nodeCoords( elemIt->_nodes[ iLeft ])[ 1 ];
+ double xPL = xL - xP, yPL = yL - yP; // components of prev-left vector
+ double xLN = xN - xL, yLN = yN - yL; // components of left-next vector
+ // normalise y of the vectors
+ double modPL = sqrt ( xPL * xPL + yPL * yPL );
+ double modLN = sqrt ( xLN * xLN + yLN * yLN );
+ if ( modLN > std::numeric_limits<double>::min() &&
+ modPL > std::numeric_limits<double>::min() )
+ {
+ yPL /= modPL;
+ yLN /= modLN;
+ // summary direction of neighboring links must be positive
+ bool clockwise = ( yPL + yLN > 0 );
+ if ( !clockwise )
+ reverse( *elemIt, swapVec );
+ }
+ }
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Fix connectivity of elements in 3D space
+ */
+//================================================================================
+
+void IntermediateMED::orientElements3D()
+{
+ // set _reverse flags of faces
+ orientFaces3D();
+
+ // -----------------
+ // fix connectivity
+ // -----------------
+
+ set<Cell>::const_iterator elemIt, elemEnd;
+ vector< pair<int,int> > swapVec;
+
+ for ( int dim = 1; dim <= 3; ++dim )
+ {
+ CellsByDimIterator cellsIt( *this, dim );
+ while ( const set<Cell > * elems = cellsIt.nextType() )
+ {
+ TCellType cellType = cellsIt.type();
+ bool isQuadratic = getGibi2MedQuadraticInterlace( cellType );
+ getReverseVector( cellType, swapVec );
+
+ elemIt = elems->begin(), elemEnd = elems->end();
+ for ( ; elemIt != elemEnd; elemIt++ )
+ {
+ // GIBI connectivity -> MED one
+ if( isQuadratic )
+ ConvertQuadratic( cellType, *elemIt );
+
+ // reverse faces
+ if ( elemIt->_reverse )
+ reverse ( *elemIt, swapVec );
+ }
+ }
+ }
+
+ orientVolumes();
+}
+
+//================================================================================
+/*!
+ * \brief Orient equally (by setting _reverse flag) all connected faces in 3D space
+ */
+//================================================================================
+
+void IntermediateMED::orientFaces3D()
+{
+ // fill map of links and their faces
+ set<const Cell*> faces;
+ map<const Cell*, Group*> fgm;
+ map<Link, list<const Cell*> > linkFacesMap;
+ map<Link, list<const Cell*> >::iterator lfIt, lfIt2;
+
+ for (size_t i=0; i!=_groups.size(); ++i)
+ {
+ Group& grp = _groups[i];
+ if ( !grp._cells.empty() && getDimension( grp._cellType ) == 2 )
+ for ( size_t j = 0; j < grp._cells.size(); ++j )
+ if ( faces.insert( grp._cells[j] ).second )
+ {
+ for ( size_t k = 0; k < grp._cells[j]->_nodes.size(); ++k )
+ linkFacesMap[ grp._cells[j]->link( k ) ].push_back( grp._cells[j] );
+ fgm.insert( make_pair( grp._cells[j], &grp ));
+ }
+ }
+ // dump linkFacesMap
+ // for ( lfIt = linkFacesMap.begin(); lfIt!=linkFacesMap.end(); lfIt++) {
+ // cout<< "LINK: " << lfIt->first.first << "-" << lfIt->first.second << endl;
+ // list<const Cell*> & fList = lfIt->second;
+ // list<const Cell*>::iterator fIt = fList.begin();
+ // for ( ; fIt != fList.end(); fIt++ )
+ // cout << "\t" << **fIt << fgm[*fIt]->nom << endl;
+ // }
+
+ // Each oriented link must appear in one face only, else a face is reversed.
+
+ queue<const Cell*> faceQueue; /* the queue contains well oriented faces
+ whose neighbors orientation is to be checked */
+ bool manifold = true;
+ while ( !linkFacesMap.empty() )
+ {
+ if ( faceQueue.empty() )
+ {
+ assert( !linkFacesMap.begin()->second.empty() );
+ faceQueue.push( linkFacesMap.begin()->second.front() );
+ }
+ while ( !faceQueue.empty() )
+ {
+ const Cell* face = faceQueue.front();
+ faceQueue.pop();
+
+ // loop on links of <face>
+ for ( int i = 0; i < (int)face->_nodes.size(); ++i )
+ {
+ Link link = face->link( i );
+ // find the neighbor faces
+ lfIt = linkFacesMap.find( link );
+ int nbFaceByLink = 0;
+ list< const Cell* > ml;
+ if ( lfIt != linkFacesMap.end() )
+ {
+ list<const Cell*> & fList = lfIt->second;
+ list<const Cell*>::iterator fIt = fList.begin();
+ assert( fIt != fList.end() );
+ for ( ; fIt != fList.end(); fIt++, nbFaceByLink++ )
+ {
+ ml.push_back( *fIt );
+ if ( *fIt != face ) // wrongly oriented neighbor face
+ {
+ const Cell* badFace = *fIt;
+ // reverse and remove badFace from linkFacesMap
+ for ( int j = 0; j < (int)badFace->_nodes.size(); ++j )
+ {
+ Link badlink = badFace->link( j );
+ if ( badlink == link ) continue;
+ lfIt2 = linkFacesMap.find( badlink );
+ if ( lfIt2 != linkFacesMap.end() )
+ {
+ list<const Cell*> & ff = lfIt2->second;
+ ff.erase( find( ff.begin(), ff.end(), badFace ));
+ if ( ff.empty() )
+ linkFacesMap.erase( lfIt2 );
+ }
+ }
+ badFace->_reverse = true; // reverse
+ //INFOS_MED( "REVERSE " << *badFace );
+ faceQueue.push( badFace );
+ }
+ }
+ linkFacesMap.erase( lfIt );
+ }
+ // add good neighbors to the queue
+ Link revLink( link.second, link.first );
+ lfIt = linkFacesMap.find( revLink );
+ if ( lfIt != linkFacesMap.end() )
+ {
+ list<const Cell*> & fList = lfIt->second;
+ list<const Cell*>::iterator fIt = fList.begin();
+ for ( ; fIt != fList.end(); fIt++, nbFaceByLink++ )
+ {
+ ml.push_back( *fIt );
+ if ( *fIt != face )
+ faceQueue.push( *fIt );
+ }
+ linkFacesMap.erase( lfIt );
+ }
+ if ( nbFaceByLink > 2 )
+ {
+ if ( manifold )
+ {
+ list<const Cell*>::iterator i = ml.begin();
+ cout << nbFaceByLink << " faces by 1 link:";
+ for( ; i!= ml.end(); i++ )
+ cout << "in sub-mesh " << fgm[ *i ]->_name << endl << **i;
+ }
+ manifold = false;
+ }
+ } // loop on links of the being checked face
+ } // loop on the face queue
+ } // while ( !linkFacesMap.empty() )
+
+ if ( !manifold )
+ cout << " -> Non manifold mesh, faces orientation may be incorrect" << endl;
+}
+
+//================================================================================
+/*!
+ * \brief Orient volumes according to MED conventions:
+ * normal of a bottom (first) face should be outside
+ */
+//================================================================================
+
+void IntermediateMED::orientVolumes()
+{
+ set<Cell>::const_iterator elemIt, elemEnd;
+ vector< pair<int,int> > swapVec;
+
+ CellsByDimIterator cellsIt( *this, 3 );
+ while ( const set<Cell > * elems = cellsIt.nextType() )
+ {
+ TCellType cellType = cellsIt.type();
+ elemIt = elems->begin(), elemEnd = elems->end();
+ int nbBottomNodes = 0;
+ switch ( cellType )
+ {
+ case NORM_TETRA4:
+ case NORM_TETRA10:
+ case NORM_PENTA6:
+ case NORM_PENTA15:
+ nbBottomNodes = 3; break;
+ case NORM_PYRA5:
+ case NORM_PYRA13:
+ case NORM_HEXA8:
+ case NORM_HEXA20:
+ nbBottomNodes = 4; break;
+ default: continue;
+ }
+ getReverseVector( cellType, swapVec );
+
+ for ( ; elemIt != elemEnd; elemIt++ )
+ {
+ // find a normal to the bottom face
+ const double* n[4];
+ n[0] = nodeCoords( elemIt->_nodes[0]); // 3 bottom nodes
+ n[1] = nodeCoords( elemIt->_nodes[1]);
+ n[2] = nodeCoords( elemIt->_nodes[2]);
+ n[3] = nodeCoords( elemIt->_nodes[nbBottomNodes]); // a top node
+ double vec01[3]; // vector n[0]-n[1]
+ vec01[0] = n[1][0] - n[0][0];
+ vec01[1] = n[1][1] - n[0][1];
+ vec01[2] = n[1][2] - n[0][2];
+ double vec02 [3]; // vector n[0]-n[2]
+ vec02[0] = n[2][0] - n[0][0];
+ vec02[1] = n[2][1] - n[0][1];
+ vec02[2] = n[2][2] - n[0][2];
+ double normal [3]; // vec01 ^ vec02
+ normal[0] = vec01[1] * vec02[2] - vec01[2] * vec02[1];
+ normal[1] = vec01[2] * vec02[0] - vec01[0] * vec02[2];
+ normal[2] = vec01[0] * vec02[1] - vec01[1] * vec02[0];
+ // check if the 102 angle is convex
+ if ( nbBottomNodes > 3 )
+ {
+ const double* n3 = nodeCoords( elemIt->_nodes[nbBottomNodes-1] );// last bottom node
+ double vec03 [3]; // vector n[0]-n3
+ vec03[0] = n3[0] - n[0][0];
+ vec03[1] = n3[1] - n[0][1];
+ vec03[2] = n3[2] - n[0][2];
+ if ( fabs( normal[0]+normal[1]+normal[2] ) <= numeric_limits<double>::max() ) // vec01 || vec02
+ {
+ normal[0] = vec01[1] * vec03[2] - vec01[2] * vec03[1]; // vec01 ^ vec03
+ normal[1] = vec01[2] * vec03[0] - vec01[0] * vec03[2];
+ normal[2] = vec01[0] * vec03[1] - vec01[1] * vec03[0];
+ }
+ else
+ {
+ double vec [3]; // normal ^ vec01
+ vec[0] = normal[1] * vec01[2] - normal[2] * vec01[1];
+ vec[1] = normal[2] * vec01[0] - normal[0] * vec01[2];
+ vec[2] = normal[0] * vec01[1] - normal[1] * vec01[0];
+ double dot2 = vec[0]*vec03[0] + vec[1]*vec03[1] + vec[2]*vec03[2]; // vec*vec03
+ if ( dot2 < 0 ) // concave -> reverse normal
+ {
+ normal[0] *= -1;
+ normal[1] *= -1;
+ normal[2] *= -1;
+ }
+ }
+ }
+ // direction from top to bottom
+ double tbDir[3];
+ tbDir[0] = n[0][0] - n[3][0];
+ tbDir[1] = n[0][1] - n[3][1];
+ tbDir[2] = n[0][2] - n[3][2];
+
+ // compare 2 directions: normal and top-bottom
+ double dot = normal[0]*tbDir[0] + normal[1]*tbDir[1] + normal[2]*tbDir[2];
+ if ( dot < 0. ) // need reverse
+ reverse( *elemIt, swapVec );
+
+ } // loop on volumes of one geometry
+ } // loop on 3D geometry types
+
+}
+
+//================================================================================
+/*!
+ * \brief Assign new IDs to nodes by skipping not used nodes and return their number
+ */
+//================================================================================
+
+int NodeContainer::numberNodes()
+{
+ int id = 1;
+ for ( size_t i = 0; i < _nodes.size(); ++i )
+ for ( size_t j = 0; j < _nodes[i].size(); ++j )
+ if ( _nodes[i][j].isUsed() )
+ _nodes[i][j]._number = id++;
+ return id-1;
+}
+
+
+//================================================================================
+/*!
+ * \brief Assign new IDs to elements
+ */
+//================================================================================
+
+void IntermediateMED::numberElements()
+{
+ set<Cell>::const_iterator elemIt, elemEnd;
+
+ // numbering _cells of type NORM_POINT1 by node number
+ {
+ const set<Cell>& points = _cellsByType[ INTERP_KERNEL::NORM_POINT1 ];
+ elemIt = points.begin(), elemEnd = points.end();
+ for ( ; elemIt != elemEnd; ++elemIt )
+ elemIt->_number = elemIt->_nodes[0]->_number;
+ }
+
+ // numbering 1D-3D _cells
+ for ( int dim = 1; dim <= 3; ++dim )
+ {
+ // check if re-numeration is needed (to try to keep elem oreder as in sauve file )
+ bool ok = true, renumEntity = false;
+ CellsByDimIterator cellsIt( *this, dim );
+ int prevNbElems = 0;
+ while ( const set<Cell> * typeCells = cellsIt.nextType() )
+ {
+ TID minNumber = std::numeric_limits<TID>::max(), maxNumber = 0;
+ for ( elemIt = typeCells->begin(), elemEnd = typeCells->end(); elemIt!=elemEnd; ++elemIt)
+ {
+ if ( elemIt->_number < minNumber ) minNumber = elemIt->_number;
+ if ( elemIt->_number > maxNumber ) maxNumber = elemIt->_number;
+ }
+ TID typeSize = typeCells->size();
+ if ( typeSize != maxNumber - minNumber + 1 )
+ ok = false;
+ if ( prevNbElems != 0 ) {
+ if ( minNumber == 1 )
+ renumEntity = true;
+ else if ( prevNbElems+1 != (int)minNumber )
+ ok = false;
+ }
+ prevNbElems += typeSize;
+ }
+
+ if ( ok && renumEntity ) // each geom type was numerated separately
+ {
+ cellsIt.init( dim );
+ prevNbElems = cellsIt.nextType()->size(); // no need to renumber the first type
+ while ( const set<Cell> * typeCells = cellsIt.nextType() )
+ {
+ for ( elemIt = typeCells->begin(), elemEnd = typeCells->end(); elemIt!=elemEnd; ++elemIt)
+ elemIt->_number += prevNbElems;
+ prevNbElems += typeCells->size();
+ }
+ }
+ if ( !ok )
+ {
+ int cellID=1;
+ cellsIt.init( dim );
+ while ( const set<Cell> * typeCells = cellsIt.nextType() )
+ for ( elemIt = typeCells->begin(), elemEnd = typeCells->end(); elemIt!=elemEnd; ++elemIt)
+ elemIt->_number = cellID++;
+ }
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Creates coord array
+ */
+//================================================================================
+
+ParaMEDMEM::DataArrayDouble * IntermediateMED::getCoords()
+{
+ DataArrayDouble* coordArray = DataArrayDouble::New();
+ coordArray->alloc( _nbNodes, _spaceDim );
+ double * coordPrt = coordArray->getPointer();
+ for ( int i = 0, nb = _points.size(); i < nb; ++i )
+ {
+ Node* n = getNode( i+1 );
+ if ( n->isUsed() )
+ {
+ const double* nCoords = nodeCoords( n );
+ std::copy( nCoords, nCoords+_spaceDim, coordPrt );
+ coordPrt += _spaceDim;
+ }
+ }
+ return coordArray;
+}
+
+//================================================================================
+/*!
+ * \brief Sets connectivity of elements to the mesh
+ * \param mesh - mesh to fill in
+ * \param coords - coordinates that must be shared by all meshes of different dim
+ */
+//================================================================================
+
+void IntermediateMED::setConnectivity( ParaMEDMEM::MEDFileUMesh* mesh,
+ ParaMEDMEM::DataArrayDouble* coords )
+{
+ int meshDim = 0;
+
+ mesh->setCoords( coords );
+
+ set<Cell>::const_iterator elemIt, elemEnd;
+ for ( int dim = 3; dim > 0; --dim )
+ {
+ CellsByDimIterator dimCells( *this, dim );
+
+ int nbOfCells = 0;
+ while ( const std::set<Cell > * cells = dimCells.nextType() )
+ nbOfCells += cells->size();
+ if ( nbOfCells == 0 )
+ continue;
+
+ if ( !meshDim ) meshDim = dim;
+
+ MEDCouplingUMesh* dimMesh = MEDCouplingUMesh::New();
+ dimMesh->setCoords( coords );
+ dimMesh->setMeshDimension( dim );
+ dimMesh->allocateCells( nbOfCells );
+
+ int prevNbCells = 0;
+ dimCells.init( dim );
+ while ( const std::set<Cell > * cells = dimCells.nextType() )
+ {
+ // fill connectivity array to take into account order of elements in the sauv file
+ const int nbCellNodes = cells->begin()->_nodes.size();
+ vector< TID > connectivity( cells->size() * nbCellNodes );
+ int * nodalConnOfCell;
+ for ( elemIt = cells->begin(), elemEnd = cells->end(); elemIt != elemEnd; ++elemIt )
+ {
+ const Cell& cell = *elemIt;
+ const int index = cell._number - 1 - prevNbCells;
+ nodalConnOfCell = &connectivity[ index * nbCellNodes ];
+ if ( cell._reverse )
+ for ( int i = nbCellNodes-1; i >= 0; --i )
+ *nodalConnOfCell++ = cell._nodes[i]->_number - 1;
+ else
+ for ( int i = 0; i < nbCellNodes; ++i )
+ *nodalConnOfCell++ = cell._nodes[i]->_number - 1;
+ }
+ prevNbCells += cells->size();
+
+ // fill dimMesh
+ TCellType cellType = dimCells.type();
+ nodalConnOfCell = &connectivity[0];
+ for ( size_t i = 0; i < cells->size(); ++i, nodalConnOfCell += nbCellNodes )
+ dimMesh->insertNextCell( cellType, nbCellNodes, nodalConnOfCell );
+ }
+ dimMesh->finishInsertingCells();
+ mesh->setMeshAtLevel( dim - meshDim, dimMesh );
+ dimMesh->decrRef();
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Fill in the mesh with groups
+ * \param mesh - mesh to fill in
+ */
+//================================================================================
+
+void IntermediateMED::setGroups( ParaMEDMEM::MEDFileUMesh* mesh )
+{
+ const int meshDim = mesh->getMeshDimension();
+ for ( int dim = 0; dim <= meshDim; ++dim )
+ {
+ const int meshDimRelToMaxExt = ( dim == 0 ? 1 : dim - meshDim );
+
+ vector<const DataArrayInt *> medGroups;
+ vector<MEDCouplingAutoRefCountObjectPtr<DataArrayInt> > refGroups;
+ for ( size_t i = 0; i < _groups.size(); ++i )
+ {
+ Group& grp = _groups[i];
+ if ( getDim( &grp ) != dim )
+ continue;
+ // convert only named groups or field supports
+ if ( grp.empty() || (grp._name.empty() && !grp._isProfile ))
+ continue;
+ //if ( grp._medGroup ) continue; // already converted
+
+ // sort cells by ID and remember their initial order in the group
+ TCellToOrderMap cell2order;
+ unsigned orderInGroup = 0;
+ vector< Group* > groupVec;
+ if ( grp._groups.empty() ) groupVec.push_back( & grp );
+ else groupVec = grp._groups;
+ for ( size_t iG = 0; iG < groupVec.size(); ++iG )
+ {
+ Group* aG = groupVec[ iG ];
+ for ( size_t iC = 0; iC < aG->_cells.size(); ++iC )
+ cell2order.insert( cell2order.end(), make_pair( aG->_cells[iC], orderInGroup++ ));
+ }
+ bool isSelfIntersect = ( orderInGroup != cell2order.size() );
+ if ( isSelfIntersect ) // self intersecting group
+ {
+ ostringstream msg;
+ msg << "Self intersecting sub-mesh: id = " << i+1
+ << ", name = |" << grp._name << "|" << endl
+ << " nb unique elements = " << cell2order.size() << endl
+ << " total nb elements = " << orderInGroup;
+ if ( grp._isProfile )
+ {
+ THROW_IK_EXCEPTION( msg.str() );
+ }
+ else
+ {
+ cout << msg << endl;
+ }
+ }
+ // create a med group
+ grp._medGroup = DataArrayInt::New();
+ grp._medGroup->setName( grp._name.c_str() );
+ grp._medGroup->alloc( orderInGroup, /*nbOfCompo=*/1 );
+ int * idsPrt = grp._medGroup->getPointer();
+ TCellToOrderMap::iterator cell2orderIt, cell2orderEnd = cell2order.end();
+ for ( cell2orderIt = cell2order.begin(); cell2orderIt != cell2orderEnd; ++cell2orderIt )
+ *idsPrt++ = (*cell2orderIt).first->_number - 1;
+
+ // try to set the mesh name
+ if ( dim == meshDim &&
+ !grp._name.empty() &&
+ grp.size() == mesh->getSizeAtLevel( meshDimRelToMaxExt ))
+ {
+ mesh->setName( grp._name.c_str() );
+ }
+ else if ( !grp._name.empty() )
+ {
+ medGroups.push_back( grp._medGroup );
+ }
+ // set relocation table
+ setRelocationTable( &grp, cell2order );
+
+ // Issue 0021311. Use case: a gibi group has references (recorded in pile 1)
+ // and several names (pile 27) refer (pile 10) to this group.
+ // We create a copy of this group per each named reference
+ for ( unsigned iRef = 0 ; iRef < grp._refNames.size(); ++iRef )
+ if ( !grp._refNames[ iRef ].empty() )
+ {
+ refGroups.push_back( grp._medGroup->deepCpy() );
+ refGroups.back()->setName( grp._refNames[ iRef ].c_str() );
+ medGroups.push_back( refGroups.back() );
+ }
+ }
+ mesh->setGroupsAtLevel( meshDimRelToMaxExt, medGroups );
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Return true if the group is on all elements and return its relative dimension
+ */
+//================================================================================
+
+bool IntermediateMED::isOnAll( const Group* grp, int & dimRel ) const
+{
+ int dim = getDim( grp );
+
+ int nbElems = 0;
+ CellsByDimIterator dimCells( *this, dim );
+ while ( const set<Cell > * cells = dimCells.nextType() )
+ nbElems += cells->size();
+
+ const bool onAll = ( nbElems == grp->size() );
+
+ if ( dim == 0 )
+ dimRel = 0;
+ else
+ {
+ int meshDim = 3;
+ for ( ; meshDim > 0; --meshDim )
+ {
+ dimCells.init( meshDim );
+ if ( dimCells.nextType() )
+ break;
+ }
+ dimRel = dim - meshDim;
+ }
+ return onAll;
+}
+
+//================================================================================
+/*!
+ * \brief Makes fields from own data
+ */
+//================================================================================
+
+ParaMEDMEM::MEDFileFields * IntermediateMED::makeMEDFileFields(ParaMEDMEM::MEDFileUMesh* mesh)
+{
+ if ( _nodeFields.empty() && _cellFields.empty() ) return 0;
+
+ // set long names
+ set< string > usedFieldNames;
+ setFieldLongNames(usedFieldNames);
+
+ MEDFileFields* fields = MEDFileFields::New();
+
+ for ( size_t i = 0; i < _nodeFields.size(); ++i )
+ setFields( _nodeFields[i], fields, mesh, i+1, usedFieldNames );
+
+ for ( size_t i = 0; i < _cellFields.size(); ++i )
+ setFields( _cellFields[i], fields, mesh, i+1, usedFieldNames );
+
+ return fields;
+}
+
+//================================================================================
+/*!
+ * \brief Make med fields from a SauvUtilities::DoubleField
+ */
+//================================================================================
+
+void IntermediateMED::setFields( SauvUtilities::DoubleField* fld,
+ ParaMEDMEM::MEDFileFields* medFields,
+ ParaMEDMEM::MEDFileUMesh* mesh,
+ const TID castemID,
+ set< string >& usedFieldNames)
+{
+ if ( !fld || !fld->isMedCompatible() ) return;
+
+ // if ( !fld->hasCommonSupport() ):
+ // each sub makes MEDFileFieldMultiTS
+ // else:
+ // unite several subs into a MEDCouplingFieldDouble
+
+ const bool uniteSubs = fld->hasCommonSupport();
+ if ( !uniteSubs )
+ cout << "Castem field #" << castemID << " " << fld->_name
+ << " is incompatible with MED format, so we split it into several fields" << endl;
+
+ for ( size_t iSub = 0; iSub < fld->_sub.size(); )
+ {
+ // set field name
+ if ( !uniteSubs || fld->_name.empty() )
+ makeFieldNewName( usedFieldNames, fld );
+
+ // allocate values
+ DataArrayDouble * values = DataArrayDouble::New();
+ values->alloc( fld->getNbTuples(iSub), fld->_sub[iSub].nbComponents() );
+
+ // set values
+ double * valPtr = values->getPointer();
+ if ( uniteSubs )
+ {
+ int nbElems = fld->_group->size();
+ for ( int elemShift = 0; elemShift < nbElems; )
+ elemShift += fld->setValues( valPtr, iSub++, elemShift );
+ setTS( fld, values, medFields, mesh );
+ }
+ else
+ {
+ fld->setValues( valPtr, iSub++ );
+ setTS( fld, values, medFields, mesh, iSub );
+ }
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Store value array of a field into med fields
+ */
+//================================================================================
+
+void IntermediateMED::setTS( SauvUtilities::DoubleField* fld,
+ ParaMEDMEM::DataArrayDouble* values,
+ ParaMEDMEM::MEDFileFields* medFields,
+ ParaMEDMEM::MEDFileUMesh* mesh,
+ const int iSub)
+{
+ // analyze a field support
+ const Group* support = fld->getSupport();
+ int dimRel;
+ const bool onAll = isOnAll( support, dimRel );
+ if ( !onAll && support->_name.empty() )
+ {
+ const_cast<Group*>(support)->_name += "PFL_" + fld->_name;
+ support->_medGroup->setName( support->_name.c_str() );
+ }
+
+ // make a time-stamp
+ MEDCouplingFieldDouble * timeStamp = MEDCouplingFieldDouble::New( fld->getMedType(),
+ fld->getMedTimeDisc() );
+ timeStamp->setName( fld->_name.c_str() );
+ timeStamp->setDescription( fld->_description.c_str() );
+ MEDCouplingAutoRefCountObjectPtr< MEDCouplingUMesh > dimMesh = mesh->getMeshAtLevel( dimRel );
+ timeStamp->setMesh( dimMesh );
+ for ( size_t i = 0; i < fld->_sub[iSub].nbComponents(); ++i )
+ values->setInfoOnComponent( i, fld->_sub[iSub]._comp_names[ i ].c_str() );
+ timeStamp->setArray( values );
+ values->decrRef();
+
+ // get a field to add the time-stamp
+ bool isNewMedField = false;
+ if ( !fld->_curMedField || fld->_name != fld->_curMedField->getName() )
+ {
+ fld->_curMedField = MEDFileFieldMultiTS::New();
+ isNewMedField = true;
+ }
+
+ // set an order
+ timeStamp->setOrder( fld->_curMedField->getNumberOfTS() );
+
+ // add the time-stamp
+ if ( onAll )
+ fld->_curMedField->appendFieldNoProfileSBT( timeStamp );
+ else
+ fld->_curMedField->appendFieldProfile( timeStamp, mesh, dimRel, support->_medGroup );
+ timeStamp->decrRef();
+
+ if ( isNewMedField ) // timeStamp must be added before this
+ {
+ medFields->pushField( fld->_curMedField );
+ fld->_curMedField->decrRef();
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Make a new unique name for a field
+ */
+//================================================================================
+
+void IntermediateMED::makeFieldNewName(std::set< std::string >& usedNames,
+ SauvUtilities::DoubleField* fld )
+{
+ string base = fld->_name;
+ if ( base.empty() )
+ {
+ base = "F_";
+ }
+ else
+ {
+ string::size_type pos = base.rfind('_');
+ if ( pos != string::npos )
+ base = base.substr( 0, pos+1 );
+ else
+ base += '_';
+ }
+
+ int i = 1;
+ do
+ {
+ fld->_name = base + SauvUtilities::toString( i++ );
+ }
+ while( !usedNames.insert( fld->_name ).second );
+}
+
+//================================================================================
+/*!
+ * \brief Return a vector ready to fill in
+ */
+//================================================================================
+
+std::vector< double >& DoubleField::addComponent( int nb_values )
+{
+ _comp_values.push_back( std::vector< double >() );
+ std::vector< double >& res = _comp_values.back();
+ res.resize( nb_values );
+ return res;
+}
+//================================================================================
+/*!
+ * \brief Returns a supporting group
+ */
+//================================================================================
+
+const Group* DoubleField::getSupport( const int iSub ) const
+{
+ return _group ? _group : _sub[iSub]._support;
+}
+
+//================================================================================
+/*!
+ * \brief Return true if each sub-component is a time stamp
+ */
+//================================================================================
+
+bool DoubleField::isMultiTimeStamps() const
+{
+ if ( _sub.size() < 2 )
+ return false;
+ bool sameSupports = true;
+ Group* grp1 = _sub[0]._support;
+ for ( size_t i = 1; i < _sub.size() && sameSupports; ++i )
+ sameSupports = ( grp1 == _sub[i]._support );
+
+ return sameSupports;
+}
+
+//================================================================================
+/*!
+ * \brief True if the field can be converted into the med field
+ */
+//================================================================================
+
+bool DoubleField::isMedCompatible() const
+{
+ for ( size_t iSub = 0; iSub < _sub.size(); ++iSub )
+ {
+ if ( !getSupport(iSub) || !getSupport(iSub)->_medGroup )
+ THROW_IK_EXCEPTION("SauvReader INTERNAL ERROR: NULL field support");
+
+ if ( !_sub[iSub].isValidNbGauss() )
+ {
+ cout << "Skip field <" << _name << "> : different nb of gauss points in components" <<endl;
+ return false;
+ }
+ }
+ // check if there are no gauss or nbGauss() == nbCellNodes,
+ // else we lack info on gauss point localization
+ // TODO?
+ return true;
+}
+
+//================================================================================
+/*!
+ * \brief return true if all sub-components has same components and same nbGauss
+ */
+//================================================================================
+
+bool DoubleField::hasSameComponentsBySupport() const
+{
+ vector< _Sub_data >::const_iterator sub_data = _sub.begin();
+ const _Sub_data& first_sub_data = *sub_data;
+ for ( ++sub_data ; sub_data != _sub.end(); ++sub_data )
+ {
+ if ( first_sub_data._comp_names != sub_data->_comp_names )
+ return false; // diff names of components
+
+ if ( first_sub_data._nb_gauss != sub_data->_nb_gauss &&
+ first_sub_data._support->_cellType == sub_data->_support->_cellType)
+ return false; // diff nb of gauss points on same cell type
+ }
+ return true;
+}
+
+//================================================================================
+/*!
+ * \brief Return type of MEDCouplingFieldDouble
+ */
+//================================================================================
+
+ParaMEDMEM::TypeOfField DoubleField::getMedType( const int iSub ) const
+{
+ using namespace INTERP_KERNEL;
+
+ const Group* grp = hasCommonSupport() ? _group : _sub[iSub]._support;
+ if ( _sub[iSub].nbGauss() > 1 )
+ {
+ const CellModel& cm = CellModel::GetCellModel( _sub[iSub]._support->_cellType );
+ return cm.getNumberOfNodes() == _sub[iSub].nbGauss() ? ON_GAUSS_NE : ON_GAUSS_PT;
+ }
+ else
+ {
+ return getDim( grp ) == 0 ? ON_NODES : ON_CELLS;
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Return TypeOfTimeDiscretization
+ */
+//================================================================================
+
+ParaMEDMEM::TypeOfTimeDiscretization DoubleField::getMedTimeDisc() const
+{
+ return ONE_TIME;
+ // NO_TIME = 4,
+ // ONE_TIME = 5,
+ // LINEAR_TIME = 6,
+ // CONST_ON_TIME_INTERVAL = 7
+}
+
+//================================================================================
+/*!
+ * \brief Return nb tuples to be used to allocate DataArrayDouble
+ */
+//================================================================================
+
+int DoubleField::getNbTuples( const int iSub ) const
+{
+ int nb = 0;
+ if ( hasCommonSupport() && !_group->_groups.empty() )
+ for ( size_t i = 0; i < _group->_groups.size(); ++i )
+ nb += _sub[i].nbGauss() * _sub[i]._support->size();
+ else
+ nb = _sub[iSub].nbGauss() * getSupport(iSub)->size();
+ return nb;
+}
+
+//================================================================================
+/*!
+ * \brief Store values of a sub-component and return nb of elements in the iSub
+ */
+//================================================================================
+
+int DoubleField::setValues( double * valPtr, const int iSub, const int elemShift ) const
+{
+ // find values for iSub
+ int iComp = 0;
+ for ( int iS = 0; iS < iSub; ++iS )
+ iComp += _sub[iS].nbComponents();
+ const vector< double > * compValues = &_comp_values[ iComp ];
+
+ const vector< unsigned >& relocTable = getSupport( iSub )->_relocTable;
+
+ // Set values
+
+ const int nbElems = _sub[iSub]._support->size();
+ const int nbGauss = _sub[iSub].nbGauss();
+ const int nbComponents = _sub[iSub].nbComponents();
+ const int nbValsByElem = nbComponents * nbGauss;
+ // check nb values
+ int nbVals = 0;
+ for ( iComp = 0; iComp < nbComponents; ++iComp )
+ nbVals += compValues[iComp].size();
+ if ( nbVals != nbElems * nbValsByElem )
+ THROW_IK_EXCEPTION("SauvMedConvertor.cxx: support size mismatches field size");
+ // compute nb values in previous subs
+ int valsShift = 0;
+ for ( int iS = iSub-1, shift = elemShift; shift > 0; )
+ {
+ int nbE = _sub[iS]._support->size();
+ shift -= nbE;
+ valsShift += nbE * _sub[iS].nbComponents() * _sub[iS].nbGauss();
+ }
+ for ( int iE = 0; iE < nbElems; ++iE )
+ {
+ int iMed = valsShift + nbValsByElem * ( relocTable.empty() ? iE : relocTable[iE+elemShift]-elemShift );
+ for ( iComp = 0; iComp < nbComponents; ++iComp )
+ for ( int iG = 0; iG < nbGauss; ++iG )
+ valPtr[ iMed + iG * nbComponents + iComp ] = compValues[iComp][ iE * nbGauss + iG ];
+ }
+ return nbElems;
+}
+
+//================================================================================
+/*!
+ * \brief Destructor of IntermediateMED
+ */
+//================================================================================
+
+IntermediateMED::~IntermediateMED()
+{
+ for ( size_t i = 0; i < _nodeFields.size(); ++i )
+ if ( _nodeFields[i] )
+ delete _nodeFields[i];
+ _nodeFields.clear();
+
+ for ( size_t i = 0; i < _cellFields.size(); ++i )
+ if ( _cellFields[i] )
+ delete _cellFields[i];
+ _cellFields.clear();
+
+ for ( size_t i = 0; i < _groups.size(); ++i )
+ if ( _groups[i]._medGroup )
+ _groups[i]._medGroup->decrRef();
+}
+
+//================================================================================
+/*!
+ * \brief CellsByDimIterator constructor
+ */
+CellsByDimIterator::CellsByDimIterator( const IntermediateMED & medi, int dim)
+{
+ myImed = & medi;
+ init( dim );
+}
+/*!
+ * \brief Initialize iteration on cells of given dimention
+ */
+void CellsByDimIterator::init(const int dim)
+{
+ myCurType = -1;
+ myTypeEnd = INTERP_KERNEL::NORM_HEXA20 + 1;
+ myDim = dim;
+}
+/*!
+ * \brief return next set of Cell's of required dimension
+ */
+const std::set< Cell > * CellsByDimIterator::nextType()
+{
+ while ( ++myCurType < myTypeEnd )
+ if ( !myImed->_cellsByType[myCurType].empty() && ( myDim < 0 || dim(false) == myDim ))
+ return & myImed->_cellsByType[myCurType];
+ return 0;
+}
+/*!
+ * \brief return dimension of cells returned by the last or further next()
+ */
+int CellsByDimIterator::dim(const bool last) const
+{
+ int type = myCurType;
+ if ( !last )
+ while ( type < myTypeEnd && myImed->_cellsByType[type].empty() )
+ ++type;
+ return type < myTypeEnd ? getDimension( TCellType( type )) : 4;
+}
+// END CellsByDimIterator ========================================================
+
--- /dev/null
+// Copyright (C) 2007-2011 CEA/DEN, EDF R&D
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+// File : SauvMedConvertor.hxx
+// Created : Tue Aug 16 14:14:02 2011
+// Author : Edward AGAPOV (eap)
+//
+
+#ifndef __SauvMedConvertor_HXX__
+#define __SauvMedConvertor_HXX__
+
+#include "InterpKernelException.hxx"
+#include "NormalizedUnstructuredMesh.hxx"
+#include "MEDCouplingRefCountObject.hxx"
+#include "SauvUtilities.hxx"
+
+#include <vector>
+#include <set>
+#include <map>
+#include <list>
+#include <algorithm>
+
+namespace ParaMEDMEM
+{
+ class DataArrayDouble;
+ class DataArrayInt;
+ class MEDFileData;
+ class MEDFileFields;
+ class MEDFileFieldMultiTS;
+ class MEDFileUMesh;
+}
+
+namespace SauvUtilities
+{
+ class IntermediateMED;
+
+ // ==============================================================================
+ typedef int TID; // an ID countered from 1
+ typedef std::pair<TID,TID> Link; // a pair of node numbers
+
+ typedef INTERP_KERNEL::NormalizedCellType TCellType;
+
+ // ==============================================================================
+ struct Node
+ {
+ TID _number;
+ size_t _coordID;
+
+ Node():_number(0){}
+ bool isUsed() const { return _number != 0; }
+ };
+
+ // ==============================================================================
+ struct Cell
+ {
+ std::vector< Node* > _nodes;
+ mutable bool _reverse; // to reverse orienation of a face only
+ mutable TID* _sortedNodeIDs; // for comparison
+ mutable TID _number;
+
+ Cell(size_t nnNodes=0) : _nodes(nnNodes),_reverse(false),_sortedNodeIDs(0),_number(0) {}
+ Cell(const Cell& ma);
+ void init() const { if ( _sortedNodeIDs ) delete [] _sortedNodeIDs; _sortedNodeIDs = 0; }
+ ~Cell() { init(); }
+
+ const TID* getSortedNodes() const; // creates if needed and return _sortedNodeIDs
+ bool operator < (const Cell& ma) const;
+ Link link(int i) const;
+
+ private:
+ Cell& operator=(const Cell& ma);
+ };
+ std::ostream& operator << (std::ostream& os, const Cell& ma);
+
+ // ==============================================================================
+ struct Group
+ {
+ TCellType _cellType;
+ std::string _name;
+ std::vector<const Cell*> _cells;
+ std::vector< Group* > _groups; // des sous-groupes composant le Group
+ //bool _isShared; // true if any Cell was added to the mesh from other Group
+ bool _isProfile; // is a field support or not
+ std::vector<std::string> _refNames; /* names of groups referring this one;
+ _refNames is resized according to nb of references
+ while reading a group (pile 1) and it is filled with
+ names while reading long names (pile 27); each named
+ reference is converted into a copy of the medGroup
+ (issue 0021311)
+ */
+ ParaMEDMEM::DataArrayInt* _medGroup; // result of conversion
+ std::vector< unsigned > _relocTable; // for _cells[i] gives its index in _medGroup
+
+ bool empty() const { return _cells.empty() && _groups.empty(); }
+ int size() const;
+ Group():_cellType(INTERP_KERNEL::NORM_ERROR), _isProfile(false), _medGroup(NULL) {}
+ };
+
+ // ==============================================================================
+ struct DoubleField
+ {
+ // a field contains several subcomponents each referring to its own support and
+ // having several named components
+ // ----------------------------------------------------------------------------
+ struct _Sub_data // a subcomponent
+ // --------------------------------
+ {
+ Group* _support; // support
+ std::vector<std::string> _comp_names; // component names
+ std::vector<int> _nb_gauss; // nb values per element in a component
+
+ void setData( int nb_comp, Group* supp )
+ { _support = supp; _comp_names.resize(nb_comp); _nb_gauss.resize(nb_comp,1); }
+ int nbComponents() const { return _comp_names.size(); }
+ std::string & compName( int i_comp ) { return _comp_names[ i_comp ]; }
+ bool isValidNbGauss() const { return *std::max_element( _nb_gauss.begin(), _nb_gauss.end() ) ==
+ *std::min_element( _nb_gauss.begin(), _nb_gauss.end() ); }
+ int nbGauss() const { return _nb_gauss[0] ? _nb_gauss[0] : 1; }
+ bool hasGauss() const { return nbGauss() > 1; }
+ };
+ // ----------------------------------------------------------------------------
+ TID _idInFile;
+ std::string _name;
+ std::string _description; // field description
+ std::vector< _Sub_data > _sub;
+ Group* _group; /* if _group == NULL then each subcomponent makes a
+ separate med field, else all subcomponents
+ are converted into timestamps of one med field.
+ The latter is possible only if nb of components in all subs
+ is the same and supports of subcomponents do not overlap
+ */
+ std::vector< std::vector< double > > _comp_values;
+ ParaMEDMEM::MEDFileFieldMultiTS* _curMedField;
+
+ DoubleField( int nb_sub, int total_nb_comp )
+ : _sub(nb_sub), _group(NULL), _curMedField(NULL) { _comp_values.reserve( total_nb_comp ); }
+ std::vector< double >& addComponent( int nb_values ); // return a vector ready to fill in
+ bool hasCommonSupport() const { return _group; } // true if there is one support for all subs
+ bool hasSameComponentsBySupport() const;
+
+ bool isMultiTimeStamps() const;
+ bool isMedCompatible() const;
+ ParaMEDMEM::TypeOfField getMedType( const int iSub=0 ) const;
+ ParaMEDMEM::TypeOfTimeDiscretization getMedTimeDisc() const;
+ int getNbTuples( const int iSub=0 ) const;
+ int getNbValuesPerElement( const int iSub=0 ) const;
+ int getNbGauss( const int iSub=0 ) const;
+ const Group* getSupport( const int iSub=0 ) const;
+ int setValues( double * valPtr, const int iSub, const int elemShift=0 ) const;
+
+ //virtual void dump(std::ostream&) const;
+ //virtual ~DoubleField() {}
+ };
+ // ==============================================================================
+ /*!
+ * \if developper
+ * Iterator on set of Cell's of given dimension
+ * \endif
+ */
+ class CellsByDimIterator
+ {
+ public:
+ CellsByDimIterator( const IntermediateMED & medi, int dim=-1); // dim=-1 - for all dimensions
+ void init(const int dim=-1);
+
+ //!< return next set of Cell's of required dimension
+ const std::set<Cell > * nextType();
+ //!< return dimension of Cell's returned by the last or further next()
+ int dim(const bool last=true) const;
+ //!< return type of Cell's returned by the last next()
+ TCellType type() const { return TCellType( myCurType ); }
+
+ private:
+ const IntermediateMED* myImed;
+ int myCurType, myTypeEnd;
+ int myDim;
+ };
+
+ // ==============================================================================
+ /*!
+ * \if developper
+ * Container of Node's. Prevents re-allocation at addition of Node's
+ * \endif
+ */
+ class NodeContainer
+ {
+ std::vector< std::vector< Node > > _nodes;
+ public:
+ Node* getNode( const TID nID )
+ {
+ const size_t chunkSize = 1000;
+ const size_t chunkID = (nID-1) / chunkSize;
+ const size_t pos = (nID-1) % chunkSize;
+ if ( _nodes.size() < chunkID+1 )
+ {
+ _nodes.resize( chunkID+1 );
+ for ( size_t i = 0; i < _nodes.size(); ++i )
+ _nodes[i].resize( chunkSize );
+ }
+ return & _nodes[chunkID][pos];
+ }
+ bool empty() const { return _nodes.empty(); }
+ size_t size() const { return empty() ? 0 : _nodes.size() * _nodes[0].size(); }
+ int numberNodes();
+ };
+
+ // ==============================================================================
+ /*!
+ * \if developper
+ * Intermediate structure used to store data read from the Sauve format file.
+ * The structure provides functions that transform the stored data to the MED format
+ *
+ * The elements inserted in maillage are ordered in order to avoid duplicated elements.
+ * \endif
+ */
+ struct IntermediateMED
+ {
+ unsigned _spaceDim;
+ unsigned _nbNodes;
+ NodeContainer _points;
+ std::vector<double> _coords;
+ std::vector<Group> _groups;
+ std::vector<DoubleField* > _nodeFields;
+ std::vector<DoubleField* > _cellFields;
+
+ // IMP 0020434: mapping GIBI names to MED names
+ std::list<nameGIBItoMED> _listGIBItoMED_mail; // to read from table "MED_MAIL" of PILE_TABLES
+ std::list<nameGIBItoMED> _listGIBItoMED_cham; // to read from table "MED_CHAM" of PILE_TABLES
+ std::list<nameGIBItoMED> _listGIBItoMED_comp; // to read from table "MED_COMP" of PILE_TABLES
+ std::map<int,std::string> _mapStrings; // to read from PILE_STRINGS
+
+ IntermediateMED(): _spaceDim(0), _nbNodes(0) {}
+ ~IntermediateMED();
+
+ Node* getNode( TID nID ) { return _points.getNode( nID ); }
+ int getNbCellsOfType( TCellType type ) const { return _cellsByType[type].size(); }
+ const Cell* insert(TCellType type, const Cell& ma) { return &( *_cellsByType[type].insert( ma ).first ); }
+ ParaMEDMEM::MEDFileData* convertInMEDFileDS();
+
+ private:
+
+ ParaMEDMEM::MEDFileUMesh* makeMEDFileMesh();
+ ParaMEDMEM::DataArrayDouble * getCoords();
+ void setConnectivity( ParaMEDMEM::MEDFileUMesh* mesh, ParaMEDMEM::DataArrayDouble* coords );
+ void setGroups( ParaMEDMEM::MEDFileUMesh* mesh );
+ ParaMEDMEM::MEDFileFields * makeMEDFileFields(ParaMEDMEM::MEDFileUMesh* mesh);
+ void setFields( SauvUtilities::DoubleField* fld,
+ ParaMEDMEM::MEDFileFields* medFields,
+ ParaMEDMEM::MEDFileUMesh* mesh,
+ const TID castemID,
+ std::set< std::string >& usedNames);
+ void setTS( SauvUtilities::DoubleField* fld,
+ ParaMEDMEM::DataArrayDouble* values,
+ ParaMEDMEM::MEDFileFields* medFields,
+ ParaMEDMEM::MEDFileUMesh* mesh,
+ const int iSub=0);
+ void checkDataAvailability() const throw(INTERP_KERNEL::Exception);
+ void setGroupLongNames();
+ void setFieldLongNames(std::set< std::string >& usedNames);
+ void makeFieldNewName(std::set< std::string >& usedNames,
+ SauvUtilities::DoubleField* fld );
+ void decreaseHierarchicalDepthOfSubgroups();
+ void eraseUselessGroups();
+ void detectMixDimGroups();
+ void orientElements2D();
+ void orientElements3D();
+ void orientFaces3D();
+ void orientVolumes();
+ void numberElements();
+ bool isOnAll( const Group* grp, int & dimRel ) const;
+ const double* nodeCoords( const Node* n ) { return &_coords[ (n->_coordID-1) * _spaceDim ]; }
+
+ // IntermediateMED()
+ // { myNodesNumerated = myMaillesNumerated = myGroupsTreated = false; currentTypeMailles = 0; }
+ // ~IntermediateMED();
+
+ //bool myNodesNumerated, myMaillesNumerated;
+
+ // mailles groupped by geom type; use insert() for filling in and
+ // _CellsByDimIterator for exploring it
+ //std::set<_Cell> maillage;
+ std::set< Cell > _cellsByType[ INTERP_KERNEL::NORM_HEXA20 + 1 ];
+ friend class CellsByDimIterator;
+ };
+
+// ==============================================================================
+ /*!
+ * \brief ASCII sauve file reader
+ */
+ class ASCIIReader : public FileReader
+ {
+ public:
+ ASCIIReader(const char* fileName);
+ virtual ~ASCIIReader();
+ virtual bool isASCII() const;
+ virtual bool open();
+ virtual bool getNextLine (char* & line, bool raiseOEF = true );
+ virtual void initNameReading(int nbValues, int width = 8);
+ virtual void initIntReading(int nbValues);
+ virtual void initDoubleReading(int nbValues);
+ virtual bool more() const;
+ virtual void next();
+ virtual int getInt() const;
+ virtual float getFloat() const;
+ virtual double getDouble() const;
+ virtual std::string getName() const;
+ int lineNb() const { return _lineNb; }
+
+ private:
+
+ bool getLine(char* & line);
+ void init( int nbToRead, int nbPosInLine, int width, int shift = 0 );
+
+ // getting a line from the file
+ int _file;
+ char* _start; // working buffer beginning
+ char* _ptr;
+ char* _eptr;
+ int _lineNb;
+
+ // line parsing
+ int _iPos, _nbPosInLine, _width, _shift;
+ char* _curPos;
+ };
+// ==============================================================================
+ /*!
+ * \brief XDR (binary) sauve file reader
+ */
+ class XDRReader : public FileReader
+ {
+ public:
+ XDRReader(const char* fileName);
+ virtual ~XDRReader();
+ virtual bool isASCII() const;
+ virtual bool open();
+ virtual bool getNextLine (char* & line, bool raiseOEF = true );
+ virtual void initNameReading(int nbValues, int width = 8);
+ virtual void initIntReading(int nbValues);
+ virtual void initDoubleReading(int nbValues);
+ virtual bool more() const;
+ virtual void next();
+ virtual int getInt() const;
+ virtual float getFloat() const;
+ virtual double getDouble() const;
+ virtual std::string getName() const;
+
+ private:
+
+ void init( int nbToRead, int width = 0 );
+
+ FILE* _xdrs_file;
+ void* _xdrs;
+ int* _xdr_ivals;
+ double* _xdr_dvals;
+ char* _xdr_cvals;
+ int _width;
+ int _xdr_kind;
+ enum
+ {
+ _xdr_kind_null,
+ _xdr_kind_char,
+ _xdr_kind_int,
+ _xdr_kind_double
+ };
+ };
+}
+
+#endif
--- /dev/null
+// Copyright (C) 2007-2011 CEA/DEN, EDF R&D
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+// File : SauvReader.cxx
+// Created : Tue Aug 16 13:57:42 2011
+// Author : Edward AGAPOV (eap)
+//
+
+#include "SauvReader.hxx"
+
+#include "SauvMedConvertor.hxx"
+#include "MEDCouplingAutoRefCountObjectPtr.hxx"
+#include "NormalizedUnstructuredMesh.hxx"
+#include "MEDCouplingRefCountObject.hxx"
+
+#include <cstring>
+#include <sstream>
+#include <iostream>
+
+using namespace ParaMEDMEM;
+using namespace SauvUtilities;
+using namespace std;
+
+#define GIBI_EQUAL(var_str, stat_str) (strncmp (var_str, stat_str, strlen(stat_str)) == 0)
+
+//================================================================================
+/*!
+ * \brief Creates a reader of a given sauve file
+ */
+//================================================================================
+
+SauvReader* SauvReader::New(const char *fileName) throw(INTERP_KERNEL::Exception)
+{
+ if ( !fileName || strlen(fileName) < 1 ) THROW_IK_EXCEPTION("Invalid file name");
+
+ ParaMEDMEM::MEDCouplingAutoRefCountObjectPtr< SauvUtilities::FileReader> parser;
+
+ // try to open as XRD
+ parser = new XDRReader( fileName );
+ if ( parser->open() )
+ {
+ SauvReader* reader = new SauvReader;
+ reader->_fileReader = parser;
+ parser->incrRef();
+ return reader;
+ }
+
+ // try to open as ASCII
+ parser = new ASCIIReader( fileName );
+ if ( parser->open() )
+ {
+ SauvReader* reader = new SauvReader;
+ reader->_fileReader = parser;
+ parser->incrRef();
+ return reader;
+ }
+
+ THROW_IK_EXCEPTION("Unable to open file |"<< fileName << "|");
+}
+//================================================================================
+/*!
+ * \brief Destructor
+ */
+//================================================================================
+
+SauvReader::~SauvReader()
+{
+ _fileReader->decrRef();
+}
+
+//================================================================================
+/*!
+ * \brief Return current line of ASCII file to report an error
+ */
+//================================================================================
+
+std::string SauvReader::lineNb() const
+{
+ if ( isASCII() )
+ return string(" (line #") + SauvUtilities::toString
+ ( static_cast<SauvUtilities::ASCIIReader*>( _fileReader )->lineNb() ) + ")";
+
+ return "";
+}
+
+//================================================================================
+/*!
+ * \brief Reads contents of the sauve file and convert it to MEDFileData
+ */
+//================================================================================
+
+ParaMEDMEM::MEDFileData * SauvReader::loadInMEDFileDS() throw(INTERP_KERNEL::Exception)
+{
+ SauvUtilities::IntermediateMED iMed; // intermadiate DS
+ _iMed = &iMed;
+
+ char* line; // a current line
+ const char* enregistrement_type=" ENREGISTREMENT DE TYPE";
+
+ while ( getNextLine(line, /*raiseOEF=*/false)) // external loop looking for "ENREGISTREMENT DE TYPE"
+ {
+ if ( isASCII() && !GIBI_EQUAL( line, enregistrement_type ))
+ continue; // "ENREGISTREMENT DE TYPE" not found -> read the next line
+
+ // read the number of a record
+ int recordNumber;
+ if ( isASCII() )
+ recordNumber = atoi( line + strlen(enregistrement_type) + 1 );
+ else
+ recordNumber = getInt();
+
+ // read the record
+ switch ( recordNumber )
+ {
+ case 2:
+ readRecord2();
+ break;
+ case 4:
+ readRecord4();
+ break;
+ case 7:
+ readRecord7();
+ break;
+ case 5:
+ break;
+ default:
+ if ( !isASCII() )
+ THROW_IK_EXCEPTION("XDR : ENREGISTREMENT DE TYPE " << recordNumber << " not implemented!!!");
+ }
+ }
+
+ ParaMEDMEM::MEDFileData* medFileData = iMed.convertInMEDFileDS();
+
+ return medFileData;
+}
+
+//================================================================================
+/*!
+ * \brief Reads "ENREGISTREMENT DE TYPE 4"
+ */
+//================================================================================
+
+void SauvReader::readRecord4()
+{
+ if ( !isASCII() )
+ {
+ getInt(); // skip NIVEAU
+ getInt(); // skip ERREUR
+ _iMed->_spaceDim = getInt();
+ getFloat(); // skip DENSITE
+ }
+ else
+ {
+ char* line;
+ getNextLine(line);
+ const char* s = " NIVEAU 15 NIVEAU ERREUR 0 DIMENSION";
+ _iMed->_spaceDim = atoi( line + strlen( s ) + 1 );
+ if ( !GIBI_EQUAL( line, " NIVEAU" ))
+ THROW_IK_EXCEPTION( "Could not read space dimension" << lineNb() );
+ }
+ if ( _iMed->_spaceDim < 1 )
+ THROW_IK_EXCEPTION( "Invalid space dimension:" << _iMed->_spaceDim );
+}
+
+//================================================================================
+/*!
+ * \brief Reads "ENREGISTREMENT DE TYPE 7"
+ */
+//================================================================================
+
+void SauvReader::readRecord7()
+{
+ if ( !isASCII() )
+ {
+ getInt(); // skip NOMBRE INFO CASTEM2000
+ getInt(); // skip IFOUR
+ getInt(); // skip NIFOUR
+ getInt(); // skip IFOMOD
+ getInt(); // skip IECHO
+ getInt(); // skip IIMPI
+ getInt(); // skip IOSPI
+ getInt(); // skip ISOTYP
+ getInt(); // skip NSDPGE
+ }
+ else
+ {
+ // skip 3 lines:
+ // NOMBRE INFO CASTEM2000 8
+ // IFOUR 2 NIFOUR 0 IFOMOD 2 IECHO 1 IIMPI 0 IOSPI 0 ISOTYP 1
+ // NSDPGE 0
+ char* line;
+ getNextLine(line);
+ getNextLine(line);
+ getNextLine(line);
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Reads the pile number, nb of objects and nb named of objects
+ */
+//================================================================================
+
+int SauvReader::readPileNumber(int& nbNamedObjects, int& nbObjects)
+{
+ // FORMAT(' PILE NUMERO',I4,'NBRE ObjectS NOMMES',I8,'NBRE ObjectS',I8)
+ int pileNumber;
+ if ( !isASCII() )
+ {
+ initIntReading(3);
+ pileNumber = getInt(); next();
+ nbNamedObjects = getInt(); next();
+ nbObjects = getInt(); next();
+ }
+ else
+ {
+ char* line;
+ getNextLine(line);
+ const char *s1 = " PILE NUMERO", *s2 = "NBRE ObjectS NOMMES", *s3 = "NBRE ObjectS";
+ if ( ! GIBI_EQUAL( line, s1 ) )
+ THROW_IK_EXCEPTION("Could not read the pile number " << lineNb() );
+ line = line + strlen(s1);
+ pileNumber = atoi( line );
+ line = line + 4 + strlen(s2);
+ nbNamedObjects = atoi( line );
+ line = line + 8 + strlen(s3);
+ nbObjects = atoi( line );
+ }
+ if ( nbNamedObjects<0 )
+ THROW_IK_EXCEPTION("Invalid nb of named objects: " << nbNamedObjects << lineNb() );
+ if ( nbObjects<0)
+ THROW_IK_EXCEPTION("Invalid nb of objects: " << nbObjects << lineNb() );
+ if ( nbObjects<nbNamedObjects)
+ THROW_IK_EXCEPTION("In PILE " << pileNumber <<
+ " nb of objects is less than nb of named objects" << lineNb() );
+ return pileNumber;
+}
+
+//================================================================================
+/*!
+ * \brief Reads "ENREGISTREMENT DE TYPE 2"
+ */
+//================================================================================
+
+void SauvReader::readRecord2()
+{
+ if ( _iMed->_spaceDim == 0 )
+ THROW_IK_EXCEPTION("Missing ENREGISTREMENT DE TYPE 4");
+
+ // read a pile number
+ int pileNumber, nbNamedObjects, nbObjects;
+ pileNumber = readPileNumber(nbNamedObjects, nbObjects);
+
+ if ( !_encounteredPiles.insert( pileNumber ).second && // piles may repeat
+ isASCII())
+ return;
+
+ // read object names and their indices
+ vector<string> objectNames(nbNamedObjects);
+ for ( initNameReading( nbNamedObjects ); more(); next() )
+ objectNames[ index() ] = getName();
+
+ vector<int> nameIndices(nbNamedObjects);
+ for ( initIntReading( nbNamedObjects ); more(); next() )
+ nameIndices[ index() ] = getInt();
+
+ switch ( pileNumber )
+ {
+ case PILE_SOUS_MAILLAGE:
+ read_PILE_SOUS_MAILLAGE(nbObjects, objectNames, nameIndices);
+ break;
+ case PILE_NODES_FIELD:
+ read_PILE_NODES_FIELD(nbObjects, objectNames, nameIndices);
+ break;
+ case PILE_TABLES:
+ read_PILE_TABLES(nbObjects, objectNames, nameIndices);
+ break;
+ case PILE_LREEL:
+ read_PILE_LREEL(nbObjects, objectNames, nameIndices);
+ break;
+ case PILE_LOGIQUES:
+ read_PILE_LOGIQUES(nbObjects, objectNames, nameIndices);
+ break;
+ case PILE_FLOATS:
+ read_PILE_FLOATS(nbObjects, objectNames, nameIndices);
+ break;
+ case PILE_INTEGERS:
+ read_PILE_INTEGERS(nbObjects, objectNames, nameIndices);
+ break;
+ case PILE_STRINGS:
+ read_PILE_STRINGS(nbObjects, objectNames, nameIndices);
+ break;
+ case PILE_LMOTS:
+ read_PILE_LMOTS(nbObjects, objectNames, nameIndices);
+ break;
+ case PILE_NOEUDS:
+ read_PILE_NOEUDS(nbObjects, objectNames, nameIndices);
+ break;
+ case PILE_COORDONNEES:
+ read_PILE_COORDONNEES(nbObjects, objectNames, nameIndices);
+ break;
+ case PILE_MODL:
+ read_PILE_MODL(nbObjects, objectNames, nameIndices);
+ break;
+ case PILE_FIELD:
+ read_PILE_FIELD(nbObjects, objectNames, nameIndices);
+ break;
+ default:
+ if ( !isASCII() )
+ THROW_IK_EXCEPTION("XDR : reading PILE " << pileNumber << " not implemented !!!");
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Reads "PILE NUMERO 1": gibi sub-meshes that are converted into med groups
+ */
+//================================================================================
+
+void SauvReader::read_PILE_SOUS_MAILLAGE(const int nbObjects,
+ std::vector<std::string>& objectNames,
+ std::vector<int>& nameIndices)
+{
+ _iMed->_groups.reserve(nbObjects*2); // fields may add some groups
+
+ char* line;
+ map<int,int> strangeGroupType;
+ int i;
+
+ for (int object=0; object!=nbObjects; ++object) // loop on sub-groups
+ {
+ initIntReading( 5 );
+ int castemCellType = getIntNext();
+ int nbSubGroups = getIntNext();
+ int nbReferences = getIntNext();
+ int nbNodesPerElem = getIntNext();
+ int nbElements = getIntNext();
+
+ _iMed->_groups.push_back(Group());
+ SauvUtilities::Group & group = _iMed->_groups.back();
+
+ // Issue 0021311. Allocate places for names of referring groups
+ // that will be possibly filled after reading long names from
+ // PILE_TABLES and PILE_STRINGS
+ group._refNames.resize( nbReferences );
+
+ // castemCellType=0 corresponds to a sub-mesh composed of other sub-meshes
+ if (castemCellType==0 && nbSubGroups>0)
+ {
+ group._groups.resize( nbSubGroups );
+ for ( initIntReading( nbSubGroups ); more(); next() )
+ group._groups[ index() ] = & _iMed->_groups[ getInt() - 1 ];
+ //std::sort( group._groups.begin(), group._groups.end() ); // for _groups comparison in getFieldSupport()
+ }
+ // skip references
+ if ( isASCII() )
+ for ( i = 0; i < nbReferences; i += 10 ) // FORMAT(10I8)
+ getNextLine(line);
+ else
+ for (initIntReading(nbReferences); more(); next());
+
+ // skip colors
+ if ( isASCII() )
+ for ( i = 0; i < nbElements; i += 10 )
+ getNextLine(line);
+ else
+ for (initIntReading(nbElements); more(); next());
+
+ // not a composit group
+ if (castemCellType>0 && nbSubGroups==0)
+ {
+ group._cellType = SauvUtilities::gibi2medGeom(castemCellType);
+
+ initIntReading( nbElements * nbNodesPerElem );
+ if ( group._cellType == INTERP_KERNEL::NORM_ERROR ) // look for group end
+ {
+ for ( ; more(); next());
+ strangeGroupType.insert( make_pair( object, castemCellType ));
+ }
+ else
+ {
+ // if ( group._cellType == MED_POINT1 ) group._cellType = NORM_ERROR; // issue 21199
+
+ // read connectivity of elements of a group
+ SauvUtilities::Cell ma( nbNodesPerElem );
+ SauvUtilities::Node* pNode;
+ group._cells.resize( nbElements );
+ for ( i = 0; i < nbElements; ++i )
+ {
+ ma.init();
+ for ( int n = 0; n < nbNodesPerElem; ++n )
+ {
+ int nodeID = getIntNext();
+ pNode = _iMed->getNode( nodeID );
+ ma._nodes[n] = pNode;
+ _iMed->_nbNodes += ( !pNode->isUsed() );
+ pNode->_number = nodeID;
+ }
+ ma._number = _iMed->getNbCellsOfType( group._cellType ) + 1;
+ group._cells[i] = _iMed->insert( group._cellType, ma );
+ }
+ }
+ }
+ } // loop on groups
+
+ // set group names
+ for (i=0; i!=(int)objectNames.size(); ++i)
+ {
+ int grpID = nameIndices[i];
+ SauvUtilities::Group & grp = _iMed->_groups[ grpID-1 ];
+ if ( !grp._name.empty() ) // a group has several names
+ { // create a group with subgroup grp and named grp.name
+ _iMed->_groups.push_back(Group());
+ _iMed->_groups.back()._groups.push_back( &_iMed->_groups[ grpID-1 ]);
+ _iMed->_groups.back()._name = grp._name;
+ }
+ grp._name=objectNames[i];
+#ifdef _DEBUG
+ map<int,int>::iterator it = strangeGroupType.find( grpID - 1 );
+ if ( it != strangeGroupType.end() )
+ cout << "Skip " << grp._name << " of not supported CASTEM type: " << it->second << endl;
+#endif
+ }
+} // read_PILE_SOUS_MAILLAGE()
+
+//================================================================================
+/*!
+ * \brief Skip "PILE NUMERO 18" of XDR file
+ */
+//================================================================================
+
+void SauvReader::read_PILE_LREEL (const int nbObjects, std::vector<std::string>&, std::vector<int>&)
+{
+ if ( isXRD() )
+ {
+ for (int object=0; object!=nbObjects; ++object) // pour chaque Group
+ {
+ initIntReading(1);
+ int nb_vals = getIntNext();
+ initDoubleReading(nb_vals);
+ for(int i=0; i<nb_vals; i++) next();
+ }
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Skip "PILE NUMERO 24" of XDR file
+ */
+//================================================================================
+
+void SauvReader::read_PILE_LOGIQUES (const int, std::vector<std::string>&, std::vector<int>&)
+{
+ if ( isXRD() )
+ {
+ initIntReading(1);
+ int nb_vals = getIntNext();
+ initIntReading(nb_vals);
+ for(int i=0; i<nb_vals; i++) next();
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Skip "PILE NUMERO 25" of XDR file
+ */
+//================================================================================
+
+void SauvReader::read_PILE_FLOATS (const int, std::vector<std::string>&, std::vector<int>&)
+{
+ if ( isXRD() )
+ {
+ initIntReading(1);
+ int nb_vals = getIntNext();
+ initDoubleReading(nb_vals);
+ for(int i=0; i<nb_vals; i++) next();
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Skip "PILE NUMERO 26" of XDR file
+ */
+//================================================================================
+
+void SauvReader::read_PILE_INTEGERS (const int, std::vector<std::string>&, std::vector<int>&)
+{
+ if ( isXRD() )
+ {
+ initIntReading(1);
+ int nb_vals = getIntNext();
+ initIntReading(nb_vals);
+ for(int i=0; i<nb_vals; i++) next();
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Skip "PILE NUMERO 29" of XDR file
+ */
+//================================================================================
+
+void SauvReader::read_PILE_LMOTS (const int nbObjects, std::vector<std::string>&, std::vector<int>&)
+{
+ if ( isXRD() )
+ {
+ for (int object=0; object!=nbObjects; ++object) // pour chaque Group
+ {
+ initIntReading(2);
+ int len = getIntNext();
+ int nb_vals = getIntNext();
+ int nb_char = len*nb_vals;
+ int nb_char_tmp = 0;
+ int fixed_length = 71;
+ while (nb_char_tmp < nb_char)
+ {
+ int remain_len = nb_char - nb_char_tmp;
+ int width;
+ if ( remain_len > fixed_length )
+ {
+ width = fixed_length;
+ }
+ else
+ {
+ width = remain_len;
+ }
+ initNameReading(1, width);
+ next();
+ nb_char_tmp += width;
+ }
+ }
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Skip "PILE NUMERO 38" of XDR file
+ */
+//================================================================================
+
+void SauvReader::read_PILE_MODL (const int nbObjects, std::vector<std::string>&, std::vector<int>&)
+{
+ if ( isXRD() )
+ {
+ for (int object=0; object!=nbObjects; ++object) // pour chaque Group
+ {
+ // see wrmodl.eso
+ initIntReading(10);
+ int n1 = getIntNext();
+ int nm2 = getIntNext();
+ int nm3 = getIntNext();
+ int nm4 = getIntNext();
+ int nm5 = getIntNext();
+ int n45 = getIntNext();
+ /*int nm6 =*/ getIntNext();
+ /*int nm7 =*/ getIntNext();
+ next();
+ next();
+ int nm1 = n1 * n45;
+ int nm9 = n1 * 16;
+ for (initIntReading(nm1); more(); next());
+ for (initIntReading(nm9); more(); next());
+ for (initNameReading(nm5, 8); more(); next());
+ for (initNameReading(nm2, 8); more(); next());
+ for (initNameReading(nm3, 8); more(); next());
+ for (initIntReading(nm4); more(); next());
+ }
+ }
+} // Fin case pile 38
+
+//================================================================================
+/*!
+ * \brief Read "PILE NUMERO 32": links to node coordinates
+ */
+//================================================================================
+
+void SauvReader::read_PILE_NOEUDS (const int nbObjects, std::vector<std::string>&, std::vector<int>&)
+{
+ initIntReading(1);
+ int nb_indices = getIntNext();
+
+ if (nb_indices != nbObjects)
+ THROW_IK_EXCEPTION("Error of reading PILE NUMERO " << PILE_NOEUDS << lineNb() );
+
+ for ( initIntReading( nbObjects ); more(); next() )
+ {
+ int coordID = getInt();
+ _iMed->getNode( index()+1 )->_coordID = coordID;
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Read "PILE NUMERO 33": node coordinates
+ */
+//================================================================================
+
+void SauvReader::read_PILE_COORDONNEES (const int nbObjects, std::vector<std::string>&, std::vector<int>&)
+{
+ initIntReading(1);
+ int nbReals = getIntNext();
+
+ if ( nbReals < _iMed->_nbNodes*(_iMed->_spaceDim+1) )
+ THROW_IK_EXCEPTION("Erroor of reading PILE NUMERO " << PILE_COORDONNEES << lineNb() );
+
+ // there are coordinates + density for each node
+ _iMed->_coords.resize( nbReals - nbReals/(_iMed->_spaceDim+1));
+ double* coordPtr = &_iMed->_coords[0];
+
+ initDoubleReading( nbReals );
+ while ( more() )
+ {
+ for (unsigned j = 0; j < _iMed->_spaceDim; ++j, next())
+ *coordPtr++ = getDouble();
+ // skip density
+ getDouble();
+ next();
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Finds or create a Group equal to a given field support
+ */
+//================================================================================
+
+SauvUtilities::Group* SauvReader::getFieldSupport(const vector<SauvUtilities::Group*>& supports)
+{
+ SauvUtilities::Group* group = NULL;
+ set<SauvUtilities::Group*> sup_set( supports.begin(), supports.end() );
+ if (sup_set.size() == 1 ) // one or equal supports
+ {
+ group = supports[0];
+ }
+ else
+ {
+ // try to find an existing composite group with the same sub-groups
+ for ( size_t i = 0; i < _iMed->_groups.size() && !group; ++i )
+ {
+ Group & grp = _iMed->_groups[i];
+ if (sup_set.size() == grp._groups.size())
+ {
+ bool sameOrder = true;
+ for ( size_t j = 0; j < supports.size() && sameOrder; ++j )
+ sameOrder = ( supports[j] == grp._groups[ j % grp._groups.size() ]);
+ if ( sameOrder )
+ group = & _iMed->_groups[i];
+ }
+ }
+ if ( !group ) // no such a group, add a new one
+ {
+ vector<SauvUtilities::Group*> newGroups( supports.begin(),
+ supports.begin() + sup_set.size() );
+ // check if supports includes newGroups in the same order
+ bool sameOrder = true;
+ for ( size_t j = newGroups.size(); j < supports.size() && sameOrder; ++j )
+ sameOrder = ( supports[j] == newGroups[ j % newGroups.size() ]);
+ if ( sameOrder )
+ {
+ _iMed->_groups.push_back( SauvUtilities::Group() );
+ group = & _iMed->_groups.back();
+ group->_groups.swap( newGroups );
+ }
+ }
+ }
+ if ( group )
+ group->_isProfile = true;
+ return group;
+}
+
+//================================================================================
+/*!
+ * \brief set field names
+ */
+//================================================================================
+
+void SauvReader::setFieldNames(const vector<SauvUtilities::DoubleField* >& fields,
+ const vector<string>& objets_nommes,
+ const vector<int>& indices_objets_nommes)
+{
+ unsigned i;
+ for ( i = 0; i < indices_objets_nommes.size(); ++i )
+ {
+ int fieldIndex = indices_objets_nommes[ i ];
+ if ( fields[ fieldIndex - 1 ] )
+ fields[ fieldIndex - 1 ]->_name = objets_nommes[ i ];
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Read "PILE NUMERO 2": NODE FIELDS
+ */
+//================================================================================
+
+void SauvReader::read_PILE_NODES_FIELD (const int nbObjects,
+ std::vector<std::string>& objectNames,
+ std::vector<int>& nameIndices)
+{
+ _iMed->_nodeFields.resize( nbObjects, (SauvUtilities::DoubleField*) 0 );
+ for (int object=0; object!=nbObjects; ++object) // loop on fields
+ {
+ // EXAMPLE ( with no values )
+
+ // (1) 4 7 2 1
+ // (2) -88 0 3 -89 0 1 -90 0 2 -91
+ // (2) 0 1
+ // (3) FX FY FZ FZ FX FY FLX
+ // (4) 0 0 0 0 0 0 0
+ // (5) cree par muc pri
+ // (6)
+ // (7) 2
+
+ // (1): nb subcomponents, nb components(total), IFOUR, nb attributes
+ int nb_sub, total_nb_comp, nb_attr;
+ int i_sub, i_comp;
+ initIntReading( 4 );
+ nb_sub = getIntNext();
+ total_nb_comp = getIntNext();
+ next(); // ignore IFOUR
+ nb_attr = getIntNext();
+ if ( nb_sub < 0 || total_nb_comp < 0 || nb_attr < 0 )
+ THROW_IK_EXCEPTION("Error of field reading " << lineNb());
+
+ // (2) loop on subcomponents of a field, for each read
+ // (a) support, (b) number of values and (c) number of components
+ vector<Group*> supports( nb_sub );
+ vector<int> nb_values ( nb_sub );
+ vector<int> nb_comps ( nb_sub );
+ int total_nb_values = 0;
+ initIntReading( nb_sub * 3 );
+ for ( i_sub = 0; i_sub < nb_sub; ++i_sub )
+ {
+ int supId = -getIntNext(); // (a) reference to support
+ if ( supId < 1 || supId > (int)_iMed->_groups.size() )
+ THROW_IK_EXCEPTION("Wrong mesh reference: "<< supId << lineNb() );
+ supports[ i_sub ] = &_iMed->_groups[ supId-1 ]; // (a) reference to support
+
+ nb_values[ i_sub ] = getIntNext(); // (b) nb points
+ total_nb_values += nb_values[ i_sub ];
+ if ( nb_values[ i_sub ] < 0 )
+ THROW_IK_EXCEPTION(" Wrong nb of points: " << nb_values[ i_sub ] << lineNb() );
+ nb_comps[ i_sub ] = getInt(); next(); // (c) nb of components in i_sub
+ }
+
+ // create a field if there are values
+ SauvUtilities::DoubleField* fdouble = 0;
+ if ( total_nb_values > 0 )
+ fdouble = new DoubleField( nb_sub, total_nb_comp );
+ _iMed->_nodeFields[ object ] = fdouble;
+
+ // (3) component names
+ initNameReading( total_nb_comp, 4 );
+ for ( i_sub = 0; i_sub < nb_sub; ++i_sub )
+ {
+ // store support id and nb components of a sub
+ if ( fdouble )
+ fdouble->_sub[ i_sub ].setData( nb_comps[ i_sub ], supports[ i_sub ] );
+ for ( i_comp = 0; i_comp < nb_comps[ i_sub ]; ++i_comp, next() )
+ {
+ // store component name
+ string compName = getName();
+ if ( fdouble )
+ fdouble->_sub[ i_sub ].compName( i_comp ) = compName;
+ }
+ }
+ // (4) nb harmonics ( ignored )
+ for ( initIntReading( total_nb_comp ); more(); next() );
+ // (5) TYPE ( ignored )
+ for (initNameReading(1, /*length=*/71); more(); next());
+ // (6) TITRE ( ignored )
+ for (initNameReading(1, /*length=*/71); more(); next());
+ // (7) attributes ( ignored )
+ for ( initIntReading( nb_attr ); more(); next() );
+
+ for ( i_sub = 0; i_sub < nb_sub; ++i_sub )
+ {
+ // loop on components: read values
+ initDoubleReading( nb_values[ i_sub ] * nb_comps[ i_sub ] );
+ for ( i_comp = 0; i_comp < nb_comps[ i_sub ]; ++i_comp )
+ {
+ if ( fdouble )
+ {
+ vector<double>& vals = fdouble->addComponent( nb_values[ i_sub ] );
+ for ( int i = 0; more() && i < nb_values[ i_sub ]; next(), ++i )
+ vals[ i ] = getDouble();
+ }
+ else
+ {
+ for ( int i = 0; i < nb_values[ i_sub ]; next(), ++i );
+ }
+ }
+ } // loop on subcomponents of a field
+
+ // set a supporting group including all subs supports but only
+ // if all subs have the same components
+ if ( fdouble && fdouble->hasSameComponentsBySupport() )
+ fdouble->_group = getFieldSupport( supports );
+ else
+ for ( i_sub = 0; i_sub < nb_sub; ++i_sub )
+ fdouble->_sub[ i_sub ]._support->_isProfile;
+
+ } // end loop on field objects
+
+ // set field names
+ setFieldNames( _iMed->_nodeFields, objectNames, nameIndices );
+
+} // read_PILE_NODES_FIELD()
+
+//================================================================================
+/*!
+ * \brief Read "PILE NUMERO 39": FIELDS
+ */
+//================================================================================
+
+void SauvReader::read_PILE_FIELD (const int nbObjects,
+ std::vector<std::string>& objectNames,
+ std::vector<int>& nameIndices)
+{
+ // REAL EXAMPLE
+
+ // (1) 1 2 6 16
+ // (2) CARACTERISTIQUES
+ // (3) -15 317773 4 0 0 0 -2 0 3
+ // (4) 317581
+ // (5) 0
+ // (6) 317767 317761 317755 317815
+ // (7) YOUN NU H SIGY
+ // (8) REAL*8 REAL*8 REAL*8 REAL*8
+ // (9) 1 1 0 0
+ // (10) 2.00000000000000E+05
+ // (11) 1 1 0 0
+ // (12) 3.30000000000000E-01
+ // (13) 1 1 0 0
+ // (14) 1.00000000000000E+04
+ // (15) 6 706 0 0
+ // (16) 1.00000000000000E+02 1.00000000000000E+02 1.00000000000000E+02
+ // (17) 1.00000000000000E+02 1.00000000000000E+02 1.00000000000000E+02
+ // (18) ...
+
+ _iMed->_cellFields.resize( nbObjects, (SauvUtilities::DoubleField*) 0 );
+ for (int object=0; object!=nbObjects; ++object) // pour chaque field
+ {
+ initIntReading( 4 );
+ int i_sub, nb_sub = getIntNext(); // (1) <nb_sub> 2 6 <title length>
+ next(); // skip "2"
+ next(); // skip "6"
+ int title_length = getIntNext(); // <title length>
+ if ( nb_sub < 1 )
+ THROW_IK_EXCEPTION("Error of field reading: wrong nb of subcomponents " << nb_sub << lineNb() );
+
+ string description;
+ if ( title_length )
+ {
+ if ( isXRD() )
+ {
+ initNameReading(1, title_length);
+ description = getName();
+ next();
+ }
+ else
+ {
+ char* line;
+ getNextLine( line ); // (2) title
+ const int len = 72; // line length
+ description = string(line + len - title_length, title_length); // title is right justified
+ }
+ }
+ // look for a line starting with '-' : <reference to support>
+ if ( isXRD() )
+ {
+ initIntReading( nb_sub * 9 );
+ }
+ else
+ {
+ do {
+ initIntReading( nb_sub * 9 );
+ } while ( getInt() >= 0 );
+ }
+ int total_nb_comp = 0;
+ vector<Group*> supports( nb_sub );
+ vector<int> nb_comp( nb_sub );
+ for ( i_sub = 0; i_sub < nb_sub; ++i_sub )
+ { // (3)
+ int supportId = -getIntNext(); // <reference to support>
+ next(); // ignore <address>
+ nb_comp [ i_sub ] = getIntNext(); // <nb of components in the sub>
+ for ( int i = 0; i < 6; ++i ) next(); // ignore 6 ints, in example "0 0 0 -2 0 3"
+
+ if ( supportId < 1 || supportId > (int)_iMed->_groups.size() )
+ THROW_IK_EXCEPTION("Error of field reading: wrong mesh reference "<< supportId << lineNb() );
+ if ( nb_comp[ i_sub ] < 0 )
+ THROW_IK_EXCEPTION("Error of field reading: wrong nb of components " <<nb_comp[ i_sub ] << lineNb() );
+
+ supports[ i_sub ] = &_iMed->_groups[ supportId-1 ];
+ total_nb_comp += nb_comp[ i_sub ];
+ }
+ // (4) dummy strings
+ for ( initNameReading( nb_sub, 17 ); more(); next() );
+ // (5) dummy strings
+ for ( initNameReading( nb_sub ); more(); next() );
+
+ // loop on subcomponents of a field, each of which refers to
+ // a certain support and has its own number of components;
+ // read component values
+ SauvUtilities::DoubleField* fdouble = 0;
+ for ( i_sub = 0; i_sub < nb_sub; ++ i_sub )
+ {
+ vector<string> comp_names( nb_comp[ i_sub ]), comp_type( nb_comp[ i_sub ]);
+ // (6) nb_comp addresses of MELVAL structure
+ for ( initIntReading( nb_comp[ i_sub ] ); more(); next() );
+ // (7) component names
+ for ( initNameReading( nb_comp[ i_sub ] ); more(); next() )
+ comp_names[ index() ] = getName();
+ // (8) component type
+ for ( initNameReading( nb_comp[ i_sub ], 17 ); more(); next() ) // 17 is name width
+ {
+ comp_type[ index() ] = getName();
+ // component types must be the same
+ if ( index() > 0 && comp_type[ index() ] != comp_type[ index() - 1] )
+ THROW_IK_EXCEPTION( "Error of field reading: diff component types <"
+ << comp_type[ index() ] << "> != <" << comp_type[ index() - 1 ]
+ << ">" << lineNb() );
+ }
+ // now type is known, create a field, one for all subs
+ bool isReal = (nb_comp[i_sub] > 0) ? (comp_type[0] == "REAL*8") : true;
+ if ( !fdouble && total_nb_comp )
+ {
+ if ( !isReal )
+ cout << "Warning: read NOT REAL field, type <" << comp_type[0] << ">" << lineNb() << endl;
+ _iMed->_cellFields[ object ] = fdouble = new SauvUtilities::DoubleField( nb_sub, total_nb_comp );
+ fdouble->_description = description;
+ }
+ // store support id and nb components of a sub
+ if ( fdouble )
+ fdouble->_sub[ i_sub ].setData( nb_comp[ i_sub ], supports[ i_sub ]);
+ // loop on components: read values
+ for ( int i_comp = 0; i_comp < nb_comp[ i_sub ]; ++i_comp )
+ {
+ // (9) nb of values
+ initIntReading( 4 );
+ int nb_val_by_elem = getIntNext();
+ int nb_values = getIntNext();
+ next();
+ next();
+ fdouble->_sub[ i_sub ]._nb_gauss[ i_comp ] = nb_val_by_elem;
+
+ // (10) values
+ nb_values *= nb_val_by_elem;
+ if ( fdouble )
+ {
+ vector<double> & vals = fdouble->addComponent( nb_values );
+ for ( isReal ? initDoubleReading( nb_values ) : initIntReading( nb_values ); more(); next())
+ vals[ index() ] = getDouble();
+ // store component name
+ fdouble->_sub[ i_sub ].compName( i_comp ) = comp_names[ i_comp ];
+ }
+ else
+ {
+ for ( isReal ? initDoubleReading( nb_values ) : initIntReading( nb_values ); more(); next() ) ;
+ }
+ }
+ } // loop on subcomponents of a field
+
+ // set id of a group including all sub supports but only
+ // if all subs have the same nb of components
+ if ( fdouble && fdouble->hasSameComponentsBySupport() )
+ fdouble->_group = getFieldSupport( supports );
+ else
+ for ( i_sub = 0; i_sub < nb_sub; ++i_sub )
+ fdouble->_sub[ i_sub ]._support->_isProfile = true;
+
+ } // end loop on field objects
+
+ // set field names
+ setFieldNames( _iMed->_cellFields, objectNames, nameIndices );
+
+} // read_PILE_FIELD()
+
+//================================================================================
+/*!
+ * \brief Read "PILE NUMERO 10": TABLES
+ */
+//================================================================================
+
+void SauvReader::read_PILE_TABLES (const int nbObjects,
+ std::vector<std::string>& objectNames,
+ std::vector<int>& nameIndices)
+{
+ // IMP 0020434: mapping GIBI names to MED names
+
+ string table_med_mail = "MED_MAIL";
+ string table_med_cham = "MED_CHAM";
+ string table_med_comp = "MED_COMP";
+ int table_med_mail_id = -1;
+ int table_med_cham_id = -1;
+ int table_med_comp_id = -1;
+ for (size_t iname = 0; iname < objectNames.size(); iname++)
+ if (objectNames[iname] == table_med_mail) table_med_mail_id = nameIndices[iname];
+ else if (objectNames[iname] == table_med_cham) table_med_cham_id = nameIndices[iname];
+ else if (objectNames[iname] == table_med_comp) table_med_comp_id = nameIndices[iname];
+
+ if ( isASCII() )
+ if (table_med_mail_id < 0 && table_med_cham_id < 0 && table_med_comp_id < 0)
+ return;
+
+ for (int itable = 1; itable <= nbObjects; itable++)
+ {
+ // read tables "MED_MAIL", "MED_CHAM" and "MED_COMP", that keeps correspondence
+ // between GIBI names (8 symbols if any) and MED names (possibly longer)
+ initIntReading(1);
+ int nb_table_vals = getIntNext();
+ if (nb_table_vals < 0)
+ THROW_IK_EXCEPTION("Error of reading PILE NUMERO 10" << lineNb() );
+
+ int name_i_med_pile;
+ initIntReading(nb_table_vals);
+ for (int i = 0; i < nb_table_vals/4; i++)
+ {
+ if (itable == table_med_mail_id ||
+ itable == table_med_cham_id ||
+ itable == table_med_comp_id)
+ {
+ nameGIBItoMED name_i;
+ name_i_med_pile = getIntNext();
+ name_i.med_id = getIntNext();
+ name_i.gibi_pile = getIntNext();
+ name_i.gibi_id = getIntNext();
+
+ if (name_i_med_pile != PILE_STRINGS)
+ {
+ // error: med(long) names are always kept in PILE_STRINGS
+ }
+ if (itable == table_med_mail_id)
+ {
+ if (name_i.gibi_pile != PILE_SOUS_MAILLAGE) {
+ // error: must be PILE_SOUS_MAILLAGE
+ }
+ _iMed->_listGIBItoMED_mail.push_back(name_i);
+ }
+ else if (itable == table_med_cham_id)
+ {
+ if (name_i.gibi_pile != PILE_FIELD &&
+ name_i.gibi_pile != PILE_NODES_FIELD)
+ {
+ // error: must be PILE_FIELD or PILE_NODES_FIELD
+ }
+ _iMed->_listGIBItoMED_cham.push_back(name_i);
+ }
+ else if (itable == table_med_comp_id)
+ {
+ if (name_i.gibi_pile != PILE_STRINGS)
+ {
+ // error: gibi(short) names of components are kept in PILE_STRINGS
+ }
+ _iMed->_listGIBItoMED_comp.push_back(name_i);
+ }
+ }
+ else
+ {
+ // pass table
+ for ( int i = 0; i < 4; ++i ) next();
+ }
+ }
+ } // for (int itable = 0; itable < nbObjects; itable++)
+}
+
+//================================================================================
+/*!
+ * \brief Read "PILE NUMERO 27"
+ */
+//================================================================================
+
+void SauvReader::read_PILE_STRINGS (const int nbObjects,
+ std::vector<std::string>& objectNames,
+ std::vector<int>& nameIndices)
+{
+ // IMP 0020434: mapping GIBI names to MED names
+ initIntReading(2);
+ int stringLen = getIntNext();
+ int nbSubStrings = getIntNext();
+ if (nbSubStrings != nbObjects)
+ THROW_IK_EXCEPTION("Error of reading PILE NUMERO 27" << lineNb() );
+
+ string aWholeString;
+ if ( isXRD() )
+ {
+ const int fixedLength = 71;
+ while ((int)aWholeString.length() < stringLen)
+ {
+ int remainLen = stringLen - aWholeString.length();
+ int len;
+ if ( remainLen > fixedLength )
+ {
+ len = fixedLength;
+ }
+ else
+ {
+ len = remainLen;
+ }
+ initNameReading(1, len);
+ aWholeString += getName();
+ next();
+ }
+ }
+ else
+ {
+ char* line;
+ const int fixedLength = 71;
+ while ((int)aWholeString.length() < stringLen)
+ {
+ getNextLine( line );
+ int remainLen = stringLen - aWholeString.length();
+ if ( remainLen > fixedLength )
+ {
+ aWholeString += line + 1;
+ }
+ else
+ {
+ aWholeString += line + ( 72 - remainLen );
+ }
+ }
+ }
+ int prevOffset = 0;
+ int currOffset = 0;
+ initIntReading(nbSubStrings);
+ for (int istr = 1; istr <= nbSubStrings; istr++, next())
+ {
+ currOffset = getInt();
+ // fill mapStrings
+ _iMed->_mapStrings[istr] = aWholeString.substr(prevOffset, currOffset - prevOffset);
+ prevOffset = currOffset;
+ }
+}
--- /dev/null
+// Copyright (C) 2007-2011 CEA/DEN, EDF R&D
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+// File : SauvReader.hxx
+// Created : Tue Aug 16 13:04:25 2011
+// Author : Edward AGAPOV (eap)
+//
+#ifndef __SauvReader_HXX__
+#define __SauvReader_HXX__
+
+#include "MEDLoaderDefines.hxx"
+#include "InterpKernelException.hxx"
+#include "SauvUtilities.hxx"
+#include "MEDCouplingRefCountObject.hxx"
+
+#include <vector>
+#include <string>
+#include <set>
+
+namespace SauvUtilities
+{
+ class FileReader;
+ class IntermediateMED;
+ class Group;
+ class DoubleField;
+}
+namespace ParaMEDMEM
+{
+ class MEDFileData;
+
+class MEDLOADER_EXPORT SauvReader : public ParaMEDMEM::RefCountObject
+{
+ public:
+ static SauvReader* New(const char *fileName) throw(INTERP_KERNEL::Exception);
+ ParaMEDMEM::MEDFileData * loadInMEDFileDS() throw(INTERP_KERNEL::Exception);
+ ~SauvReader();
+
+ private:
+
+ void readRecord2();
+ void readRecord4();
+ void readRecord7();
+
+ int readPileNumber(int& nbNamedObjects, int& nbObjects);
+ void read_PILE_SOUS_MAILLAGE(const int nbObjects, std::vector<std::string>& objectNames, std::vector<int>& nameIndices);
+ void read_PILE_NODES_FIELD (const int nbObjects, std::vector<std::string>& objectNames, std::vector<int>& nameIndices);
+ void read_PILE_TABLES (const int nbObjects, std::vector<std::string>& objectNames, std::vector<int>& nameIndices);
+ void read_PILE_LREEL (const int nbObjects, std::vector<std::string>& objectNames, std::vector<int>& nameIndices);
+ void read_PILE_LOGIQUES (const int nbObjects, std::vector<std::string>& objectNames, std::vector<int>& nameIndices);
+ void read_PILE_FLOATS (const int nbObjects, std::vector<std::string>& objectNames, std::vector<int>& nameIndices);
+ void read_PILE_INTEGERS (const int nbObjects, std::vector<std::string>& objectNames, std::vector<int>& nameIndices);
+ void read_PILE_STRINGS (const int nbObjects, std::vector<std::string>& objectNames, std::vector<int>& nameIndices);
+ void read_PILE_LMOTS (const int nbObjects, std::vector<std::string>& objectNames, std::vector<int>& nameIndices);
+ void read_PILE_NOEUDS (const int nbObjects, std::vector<std::string>& objectNames, std::vector<int>& nameIndices);
+ void read_PILE_COORDONNEES (const int nbObjects, std::vector<std::string>& objectNames, std::vector<int>& nameIndices);
+ void read_PILE_MODL (const int nbObjects, std::vector<std::string>& objectNames, std::vector<int>& nameIndices);
+ void read_PILE_FIELD (const int nbObjects, std::vector<std::string>& objectNames, std::vector<int>& nameIndices);
+
+ SauvUtilities::Group* getFieldSupport(const std::vector<SauvUtilities::Group*>& fieldSupports);
+ void setFieldNames(const std::vector<SauvUtilities::DoubleField*>& fields,
+ const std::vector<std::string>& objectNames,
+ const std::vector<int>& nameIndices);
+
+ bool isASCII() const { return _fileReader->isASCII(); }
+ bool isXRD() const { return !isASCII(); }
+ bool getNextLine (char* & line, bool raiseOEF = true ) { return _fileReader->getNextLine( line, raiseOEF ); }
+ void initNameReading(int nbValues, int width = 8) { _fileReader->initNameReading( nbValues, width ); }
+ void initIntReading(int nbValues) { _fileReader->initIntReading( nbValues ); }
+ void initDoubleReading(int nbValues) { _fileReader->initDoubleReading( nbValues ); }
+ bool more() const { return _fileReader->more(); }
+ void next() { _fileReader->next(); }
+ int index() const { return _fileReader->index(); }
+ int getInt() const { return _fileReader->getInt(); }
+ int getIntNext() { int i = getInt(); next(); return i; }
+ float getFloat() const { return _fileReader->getFloat(); }
+ double getDouble() const { return _fileReader->getDouble(); }
+ std::string getName() const { return _fileReader->getName(); }
+ std::string lineNb() const;
+
+
+ std::set<int> _encounteredPiles;
+
+ SauvUtilities::FileReader* _fileReader;
+ SauvUtilities::IntermediateMED* _iMed;
+};
+}
+#endif
--- /dev/null
+// Copyright (C) 2007-2011 CEA/DEN, EDF R&D
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+// File : SauvUtilities.hxx
+// Created : Mon Aug 22 18:27:34 2011
+// Author : Edward AGAPOV (eap)
+
+#ifndef __SauvUtilities_HXX__
+#define __SauvUtilities_HXX__
+
+#include "MEDCouplingRefCountObject.hxx"
+#include "NormalizedUnstructuredMesh.hxx"
+
+#include <string>
+#include <sstream>
+
+#define THROW_IK_EXCEPTION(text) \
+ { \
+ std::ostringstream oss; oss << text; \
+ throw INTERP_KERNEL::Exception(oss.str().c_str()); \
+ }
+
+namespace SauvUtilities
+{
+ INTERP_KERNEL::NormalizedCellType gibi2medGeom( size_t gibiType );
+ int med2gibiGeom( INTERP_KERNEL::NormalizedCellType medGeomType );
+ const int * getGibi2MedQuadraticInterlace( INTERP_KERNEL::NormalizedCellType type );
+ unsigned getDimension( INTERP_KERNEL::NormalizedCellType type );
+
+ enum Readable_Piles
+ {
+ PILE_SOUS_MAILLAGE=1 ,
+ PILE_NODES_FIELD =2 ,
+ PILE_TABLES =10,
+ PILE_LREEL =18,
+ PILE_LOGIQUES =24,
+ PILE_FLOATS =25,
+ PILE_INTEGERS =26,
+ PILE_STRINGS =27,
+ PILE_LMOTS =29,
+ PILE_NOEUDS =32,
+ PILE_COORDONNEES =33,
+ PILE_MODL =38,
+ PILE_FIELD =39,
+ PILE_LAST_READABLE=39
+ };
+
+ //================================================================================
+ /*!
+ * \brief Converts anything to string
+ */
+ //================================================================================
+
+ template<class T> std::string toString(const T& anything)
+ {
+ std::ostringstream s; s << anything; return s.str();
+ }
+
+ // ==============================================================================
+ // IMP 0020434: mapping GIBI names to MED names
+ struct nameGIBItoMED
+ {
+ // GIBI value
+ int gibi_pile; // PILE_SOUS_MAILLAGE or PILE_FIELD/PILE_NODES_FIELD, or PILE_STRINGS(for components)
+ int gibi_id;
+ std::string gibi_name; // used only for components
+ // MED value
+ // med_pile = 27; // PILE_STRINGS
+ int med_id; // used only on reading
+ std::string med_name; // used only on writing
+ };
+
+ // ==============================================================================
+ /*!
+ * \brief Base class for ASCII and XDR file readers
+ */
+ class FileReader : public ParaMEDMEM::RefCountObject
+ {
+ public:
+ FileReader(const char* fileName);
+ virtual ~FileReader() {}
+ virtual bool isASCII() const = 0;
+
+ virtual bool open() = 0;
+ virtual bool getNextLine (char* & line, bool raiseOEF = true ) = 0;
+ virtual void initNameReading(int nbValues, int width = 8) = 0;
+ virtual void initIntReading(int nbValues) = 0;
+ virtual void initDoubleReading(int nbValues) = 0;
+ virtual bool more() const = 0;
+ virtual void next() = 0;
+ virtual int index() const { return _iRead; }
+ virtual int getInt() const = 0;
+ virtual float getFloat() const = 0;
+ virtual double getDouble() const = 0;
+ virtual std::string getName() const = 0;
+
+ protected:
+ std::string _fileName, _curLocale;
+ int _iRead, _nbToRead;
+ };
+}
+#endif
--- /dev/null
+// Copyright (C) 2007-2011 CEA/DEN, EDF R&D
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+// File : SauvWriter.cxx
+// Created : Wed Aug 24 12:55:55 2011
+// Author : Edward AGAPOV (eap)
+
+#include "SauvWriter.hxx"
+
+#include "InterpKernelException.hxx"
+#include "MEDFileMesh.hxx"
+#include "MEDFileField.hxx"
+#include "MEDFileData.hxx"
+#include "CellModel.hxx"
+
+#include <fstream>
+#include <sstream>
+#include <iostream>
+#include <cstdlib>
+#include <iomanip>
+
+using namespace ParaMEDMEM;
+using namespace SauvUtilities;
+using namespace std;
+
+#define INFOS_MED(txt) cout << txt << endl;
+
+namespace
+{
+ const char* zeroI8 = " 0"; // FORMAT(I8)
+
+ // ============================================================
+ // the class writes endl to the file as soon as <limit> fields
+ // have been written after the last endl
+ // ============================================================
+
+ class TFieldCounter
+ {
+ fstream& _file;
+ int _count, _limit;
+ public:
+ TFieldCounter(fstream& f, int limit=0): _file(f), _limit(limit) { init(); }
+ void init(int limit=0) // init, is done by stop() as well
+ { if (limit) _limit = limit; _count = 0; }
+ void operator++(int) // next
+ { if ( ++_count == _limit ) { _file << endl; init(); }}
+ void stop() // init() and write endl if there was no endl after the last written field
+ { if ( _count ) _file << endl; init(); }
+ ~TFieldCounter() { stop(); }
+ };
+
+ //================================================================================
+ /*!
+ * \brief Return a name of a field support on all elements
+ */
+ //================================================================================
+
+ string noProfileName( INTERP_KERNEL::NormalizedCellType type )
+ {
+ return "INTERP_KERNEL::NormalizedCellType_" + SauvUtilities::toString( type );
+ }
+
+ //================================================================================
+ /*!
+ * \brief Remove white spaces from the head and tail
+ */
+ //================================================================================
+
+ string cleanName( const string& theName )
+ {
+ string name = theName;
+ if ( !name.empty() )
+ {
+ // cut off leading white spaces
+ string::size_type firstChar = name.find_first_not_of(" \t");
+ if (firstChar < name.length())
+ {
+ name = name.substr(firstChar);
+ }
+ else
+ {
+ name = ""; // only whitespaces there - remove them
+ }
+ // cut off trailing white spaces
+ string::size_type lastChar = name.find_last_not_of(" \t");
+ if (lastChar < name.length())
+ name = name.substr(0, lastChar + 1);
+ }
+ return name;
+ }
+
+ //================================================================================
+ /*!
+ * \brief Converts MED long names into SAUVE short ones, returnes a healed long name
+ */
+ //================================================================================
+
+ string addName (map<string,int>& nameMap,
+ map<string,int>& namePrefixesMap,
+ const string& theName,
+ const int index)
+ {
+ // Converts names like:
+ // MED: GIBI:
+ // TEMPERATURE_FLUIDE -> TEMPE001
+ // TEMPERATURE_SOLIDE -> TEMPE002
+ // PRESSION -> PRESSION
+ // NU -> NU
+ // VOLUM001 -> VOLUM001
+ // VOLUMOFOBJECT -> VOLUM003
+ // VOLUM002 -> VOLUM002
+ string healedName = cleanName(theName);
+ int ind = index;
+
+ if (!healedName.empty())
+ {
+ string name = healedName;
+ int len = name.length();
+ for (int i = 0; i < len; ++i)
+ name[i] = toupper(name[i]);
+
+ bool doResave = false; // only for tracing
+
+ // I. Save a short name as it is
+ if (len <= 8)
+ {
+ INFOS_MED("Save <" << theName << "> as <" << name << ">");
+
+ map<string,int>::iterator it = nameMap.find(name);
+ if (it != nameMap.end())
+ {
+ // There is already such name in the map.
+
+ // a. Replace in the map the old pair by the current one
+ int old_ind = nameMap[name];
+ nameMap[name] = ind;
+ // b. Rebuild the old pair (which was in the map,
+ // it seems to be built automatically by step II)
+ ind = old_ind;
+ // continue with step II
+ doResave = true; // only for tracing
+ }
+ else
+ {
+ // Save in the map
+ nameMap.insert(make_pair(name, ind));
+
+ // Update loc_index for this name (if last free characters represents a number)
+ // to avoid conflicts with long names, same in first 5 characters
+ if (len == 8)
+ {
+ int new_loc_index = atoi(name.c_str() + 5);
+ if (new_loc_index > 0)
+ {
+ // prefix
+ string str = name.substr(0,5);
+ if (namePrefixesMap.find(str) != namePrefixesMap.end())
+ {
+ int old_loc_index = namePrefixesMap[str];
+ if (new_loc_index < old_loc_index) new_loc_index = old_loc_index;
+ }
+ namePrefixesMap[str] = new_loc_index;
+ }
+ }
+ return healedName;
+ }
+ } // if (len <= 8)
+
+ // II. Cut long name and add a numeric suffix
+
+ // first 5 or less characters of the name
+ if (len > 5) name = name.substr(0,5);
+
+ // numeric suffix
+ map<string,int>::iterator name2ind = namePrefixesMap.insert( make_pair( name, 0 )).first;
+ string numSuffix = SauvUtilities::toString( ++(name2ind->second) );
+
+ if ( numSuffix.size() + name.size() > 8 )
+ THROW_IK_EXCEPTION("Can't write not unique name: " << healedName);
+
+ if ( numSuffix.size() < 3 )
+ numSuffix.insert( 0, 3 - numSuffix.size(), '0' );
+
+ name += numSuffix;
+ nameMap.insert(make_pair(name, ind));
+
+ if (doResave)
+ {
+ INFOS_MED("Resave previous <" << healedName << "> as <" << name << ">");
+ }
+ else
+ {
+ INFOS_MED("Save <" << theName << "> as <" << name << ">");
+ }
+ }
+ return healedName;
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Creates SauvWriter
+ */
+//================================================================================
+
+SauvWriter* SauvWriter::New()
+{
+ return new SauvWriter;
+}
+
+//================================================================================
+/*!
+ * \brief Fills own DS by MEDFileData
+ */
+//================================================================================
+
+void SauvWriter::setMEDFileDS(const MEDFileData* medData,
+ unsigned meshIndex)
+{
+ if ( !medData) THROW_IK_EXCEPTION("NULL MEDFileData");
+
+ MEDFileMeshes * meshes = medData->getMeshes();
+ MEDFileFields * fields = medData->getFields();
+ if ( !meshes) THROW_IK_EXCEPTION("No meshes in MEDFileData");
+
+ _fileMesh = meshes->getMeshAtPos( meshIndex );
+ _fileMesh->incrRef();
+
+ if ( fields )
+ for ( int i = 0; i < fields->getNumberOfFields(); ++i )
+ {
+ MEDFileFieldMultiTS * f = fields->getFieldAtPos(i);
+ if ( f->getMeshName() == _fileMesh->getName() )
+ {
+ vector< vector<TypeOfField> > fTypes = f->getTypesOfFieldAvailable();
+ if ( fTypes[0].size() == 1 && fTypes[0][0] == ON_NODES )
+ _nodeFields.push_back( f );
+ else
+ _cellFields.push_back( f );
+ }
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Adds a submesh
+ */
+//================================================================================
+
+SauvWriter::SubMesh* SauvWriter::addSubMesh(const std::string& name, int dimRelExt)
+{
+ if ( _subs.capacity() < _subs.size() + 1 )
+ THROW_IK_EXCEPTION("SauvWriter: INTERNAL error, wrong evaluation of nb of sub-meshes");
+ _subs.resize( _subs.size() + 1 );
+ SubMesh& sm = _subs.back();
+ sm._name = name;
+ sm._dimRelExt = dimRelExt;
+ return &sm;
+}
+//================================================================================
+/*!
+ * \brief Returns nb of cell types
+ */
+//================================================================================
+
+int SauvWriter::SubMesh::nbTypes() const
+{
+ int nb = 0;
+ for (int i = 0; i < cellIDsByTypeSize(); ++i )
+ nb += int( !_cellIDsByType[i].empty() );
+ return nb;
+}
+
+//================================================================================
+/*!
+ * \brief Fill _subs
+ */
+//================================================================================
+
+void SauvWriter::fillSubMeshes( int& nbSauvObjects, map<string,int>& nameNbMap )
+{
+ // evaluate nb of _subs in order to avoid re-allocation of _subs
+ int nbSubs = 1; // for the very mesh
+ nbSubs += _fileMesh->getFamilyInfo().size() + 4; // + 4 zero families (for each dimRelExt)
+ nbSubs += _fileMesh->getGroupInfo().size();
+ nbSubs += evaluateNbProfileSubMeshes();
+ _subs.clear();
+ _subs.reserve( nbSubs );
+
+ fillFamilySubMeshes();
+ fillGroupSubMeshes();
+ fillProfileSubMeshes();
+
+ // fill names of SubMesh'es and count nb of sauv sub-meshes they will be stored into
+ nbSauvObjects = 0;
+ map<string,int> namePrefixMap;
+ for ( size_t i = 0; i < _subs.size(); ++i )
+ {
+ SubMesh& sm = _subs[i];
+
+ sm._nbSauvObjects = 0;
+ if ( sm._subs.empty() )
+ {
+ sm._nbSauvObjects = sm.nbTypes();
+ }
+ else
+ {
+ sm._nbSauvObjects = 1;
+ }
+
+ sm._id = nbSauvObjects+1;
+ nbSauvObjects += sm._nbSauvObjects;
+
+ if ( sm._nbSauvObjects )
+ sm._name = addName( nameNbMap, namePrefixMap, sm._name, sm._id );
+
+ if ( sm._nbSauvObjects && !sm._name.empty() )
+ {
+ nameGIBItoMED aMEDName;
+ aMEDName.gibi_pile = PILE_SOUS_MAILLAGE;
+ aMEDName.gibi_id = sm._id;
+ aMEDName.med_name = sm._name;
+ _longNames[ LN_MAIL ].push_back(aMEDName);
+ }
+ }
+}
+
+//================================================================================
+/*!
+ * \brief fill sub-meshes of families
+ */
+//================================================================================
+
+void SauvWriter::fillFamilySubMeshes()
+{
+ SubMesh* nilSm = (SubMesh*) 0;
+ std::vector<int> dims = _fileMesh->getNonEmptyLevelsExt();
+ for ( size_t iDim = 0; iDim < dims.size(); ++iDim )
+ {
+ int dimRelExt = dims[ iDim ];
+ MEDCouplingAutoRefCountObjectPtr< MEDCouplingMesh > mesh = _fileMesh->getGenMeshAtLevel(dimRelExt);
+ const DataArrayInt * famIds = _fileMesh->getFamilyFieldAtLevel(dimRelExt);
+
+ int curFamID = 0;
+ SubMesh* curSubMesh = addSubMesh( "", dimRelExt ); // submesh of zero family
+ _famIDs2Sub[0] = curSubMesh;
+ int sub0Index = _subs.size()-1;
+
+ const int * famID = famIds->begin(), * famIDEnd = famIds->end();
+ for ( int cellID = 0; famID < famIDEnd; ++famID, cellID++ )
+ {
+ if ( *famID != curFamID )
+ {
+ curFamID = *famID;
+ map< int, SubMesh* >::iterator f2s = _famIDs2Sub.insert( make_pair( curFamID, nilSm )).first;
+ if ( !f2s->second )
+ f2s->second = addSubMesh( "", dimRelExt ); // no names for families
+ curSubMesh = f2s->second;
+ }
+ INTERP_KERNEL::NormalizedCellType cellType =
+ dimRelExt == 1 ? INTERP_KERNEL::NORM_POINT1 : mesh->getTypeOfCell( cellID );
+ curSubMesh->_cellIDsByType[ cellType ].push_back( cellID );
+ }
+
+ if ( dimRelExt == 1 )
+ {
+ // clear submesh of nodal zero family
+ _famIDs2Sub[0]->_cellIDsByType[ INTERP_KERNEL::NORM_POINT1 ].clear();
+ }
+ else if ( dimRelExt == 0 )
+ {
+ // make a submesh including all cells
+ if ( sub0Index == _subs.size()-1 )
+ {
+ _famIDs2Sub[0]->_name = _fileMesh->getName(); // there is the zero family only
+ }
+ else
+ {
+ curSubMesh = addSubMesh( _fileMesh->getName(), dimRelExt );
+ if ( _famIDs2Sub[0]->nbTypes() == 0 )
+ sub0Index++; // skip an empty zero family
+ for ( size_t i = sub0Index; i < _subs.size()-1; ++i )
+ curSubMesh->_subs.push_back( & _subs[i] );
+ }
+ }
+ }
+}
+
+//================================================================================
+/*!
+ * \brief fill sub-mehses of groups
+ */
+//================================================================================
+
+void SauvWriter::fillGroupSubMeshes()
+{
+ const map<string, vector<string> >& grpFams = _fileMesh->getGroupInfo();
+ map<string, vector<string> >::const_iterator g2ff = grpFams.begin();
+ for ( ; g2ff != grpFams.end(); ++g2ff )
+ {
+ const string& groupName = g2ff->first;
+ const vector<string>& famNames = g2ff->second;
+ if ( famNames.empty() ) continue;
+ std::vector<SubMesh*> famSubMeshes( famNames.size() );
+ for ( size_t i = 0; i < famNames.size(); ++i )
+ {
+ int famID = _fileMesh->getFamilyId( famNames[i].c_str() );
+ map< int, SubMesh* >::iterator i2f = _famIDs2Sub.find( famID );
+ if ( i2f == _famIDs2Sub.end() )
+ THROW_IK_EXCEPTION("SauvWriter::fillGroupSubMeshes(): unknown family ID: " << famID);
+ famSubMeshes[ i ] = i2f->second;
+ }
+ SubMesh* grpSubMesh = addSubMesh( groupName, famSubMeshes[0]->_dimRelExt );
+ grpSubMesh->_subs.swap( famSubMeshes );
+ }
+}
+
+//================================================================================
+/*!
+ * \brief fill sub-mehses of profiles
+ */
+//================================================================================
+
+void SauvWriter::fillProfileSubMeshes()
+{
+ _profile2Sub.clear();
+ SubMesh* nilSm = (SubMesh*) 0;
+ for ( int isOnNodes = 0; isOnNodes < 2; ++isOnNodes )
+ {
+ vector< MEDCouplingAutoRefCountObjectPtr< MEDFileFieldMultiTS > >
+ fields = isOnNodes ? _nodeFields : _cellFields;
+ for ( size_t i = 0; i < fields.size(); ++i )
+ {
+ vector< pair<int,int> > iters = fields[i]->getIterations();
+
+ vector<INTERP_KERNEL::NormalizedCellType> types;
+ vector< vector<TypeOfField> > typesF;
+ vector< vector<string> > pfls, locs;
+ fields[i]->getFieldSplitedByType( iters[0].first, iters[0].second,
+ _fileMesh->getName(), types, typesF, pfls, locs);
+ int dimRelExt;
+ for ( size_t iType = 0; iType < types.size(); ++iType )
+ {
+ if ( types[iType] == INTERP_KERNEL::NORM_ERROR )
+ dimRelExt = 1; // on nodes
+ else
+ dimRelExt = getDimension( types[iType] ) - _fileMesh->getMeshDimension();
+ for ( size_t iPfl = 0; iPfl < pfls[iType].size(); ++iPfl )
+ {
+ bool isOnAll = pfls[iType][iPfl].empty();
+ if ( isOnAll ) pfls[iType][iPfl] = noProfileName( types[iType] );
+ map< string, SubMesh* >::iterator pfl2sm =
+ _profile2Sub.insert( make_pair( pfls[iType][iPfl], nilSm )).first;
+ if ( !pfl2sm->second )
+ {
+ SubMesh* sm = pfl2sm->second = addSubMesh( "", dimRelExt ); // no names for profiles
+ const DataArrayInt * pfl = isOnAll ? 0 : fields[i]->getProfile( pfls[iType][iPfl] );
+ makeProfileIDs( sm, types[iType], pfl );
+ }
+ }
+ }
+ }
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Return max possible nb of sub-meshes to decsribe field supports
+ */
+//================================================================================
+
+int SauvWriter::evaluateNbProfileSubMeshes() const
+{
+ int nb = 0;
+ for ( size_t i = 0; i < _nodeFields.size(); ++i )
+ nb += 1 + _nodeFields[i]->getPflsReallyUsed2().size();
+
+ for ( size_t i = 0; i < _cellFields.size(); ++i )
+ {
+ nb += _cellFields[i]->getPflsReallyUsed2().size();
+
+ vector< pair<int,int> > iters = _cellFields[i]->getIterations();
+
+ vector<INTERP_KERNEL::NormalizedCellType> types;
+ vector< vector<TypeOfField> > typesF;
+ vector< vector<string> > pfls, locs;
+ _cellFields[i]->getFieldSplitedByType( iters[0].first, iters[0].second,
+ _fileMesh->getName(), types, typesF, pfls, locs);
+ nb += 2 * types.size(); // x 2 - a type can be on nodes and on cells at the same time
+ }
+
+ return nb;
+}
+
+//================================================================================
+/*!
+ * \brief Transorm a profile into ids of mesh elements
+ */
+//================================================================================
+
+void SauvWriter::makeProfileIDs( SubMesh* sm,
+ INTERP_KERNEL::NormalizedCellType type,
+ const DataArrayInt* profile )
+{
+ MEDCouplingAutoRefCountObjectPtr< MEDCouplingMesh >
+ mesh = _fileMesh->getGenMeshAtLevel(sm->_dimRelExt);
+ const MEDCouplingUMesh* uMesh = dynamic_cast< const MEDCouplingUMesh* > ((const MEDCouplingMesh*) mesh );
+
+ if ( sm->_dimRelExt == 1 ) type = INTERP_KERNEL::NORM_POINT1;
+ vector< int >& ids = sm->_cellIDsByType[ type ];
+
+ if ( sm->_dimRelExt == 1 || !uMesh )
+ {
+ // profile on nodes or mesh is CARTESIAN
+ if ( profile )
+ {
+ ids.assign( profile->begin(), profile->end() );
+ }
+ else // on all
+ {
+ ids.resize( sm->_dimRelExt == 1 ? mesh->getNumberOfNodes() : mesh->getNumberOfCells() );
+ for ( size_t i = 0; i < ids.size(); ++i )
+ ids[i]=i;
+ }
+ }
+ else
+ {
+ // profile on cells
+ vector<int> code(3);
+ code[0] = type;
+ if ( profile ) // on all cells
+ {
+ code[1] = profile->getNumberOfTuples();
+ code[2] = 0;
+ }
+ else
+ {
+ code[1] = mesh->getNumberOfCellsWithType( type );
+ code[2] = -1;
+ }
+ vector<const DataArrayInt *> idsPerType( 1, profile );
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt>
+ resIDs = uMesh->checkTypeConsistencyAndContig( code, idsPerType );
+ ids.assign( resIDs->begin(), resIDs->end() );
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Write its data into the SAUVE file
+ */
+//================================================================================
+
+void SauvWriter::write(const char* fileName)
+{
+ std::fstream fileStream;
+ fileStream.open( fileName, ios::out);
+ if
+#ifdef WNT
+ ( !fileStream || !fileStream.is_open() )
+#else
+ ( !fileStream || !fileStream.rdbuf()->is_open() )
+#endif
+ THROW_IK_EXCEPTION("Can't open the file |"<<fileName<<"|");
+ _sauvFile = &fileStream;
+
+ _subs.clear();
+ _famIDs2Sub.clear();
+ _profile2Sub.clear();
+ _longNames[ LN_MAIL ].clear();
+ _longNames[ LN_CHAM ].clear();
+ _longNames[ LN_COMP ].clear();
+
+ map<string,int> fldNamePrefixMap;
+
+ writeFileHead();
+ writeSubMeshes();
+ writeNodes();
+ writeNodalFields(fldNamePrefixMap);
+ writeElemFields(fldNamePrefixMap);
+ writeLongNames();
+ writeLastRecord();
+
+ _sauvFile->close();
+}
+//================================================================================
+/*!
+ * \brief Writes "ENREGISTREMENT DE TYPE" 4 and 7
+ */
+//================================================================================
+
+void SauvWriter::writeFileHead()
+{
+ MEDCouplingAutoRefCountObjectPtr< MEDCouplingMesh > mesh = _fileMesh->getGenMeshAtLevel(0);
+
+ *_sauvFile
+ << " ENREGISTREMENT DE TYPE 4" << endl
+ << " NIVEAU 16 NIVEAU ERREUR 0 DIMENSION " << mesh->getSpaceDimension() <<endl
+ << " DENSITE 0.00000E+00" << endl
+ << " ENREGISTREMENT DE TYPE 7" << endl
+ << " NOMBRE INFO CASTEM2000 8" <<endl
+ << " IFOUR -1 NIFOUR 0 IFOMOD -1 IECHO 1 IIMPI 0 IOSPI 0 ISOTYP 1" << endl
+ << " NSDPGE 0" << endl;
+}
+
+//================================================================================
+/*!
+ * \brief Writes names of objects
+ */
+//================================================================================
+
+void SauvWriter::writeNames( const map<string,int>& nameNbMap )
+{
+ if ( !nameNbMap.empty() )
+ {
+ // write names of objects
+ // * 8001 FORMAT(8(1X,A8))
+ TFieldCounter fcount( *_sauvFile, 8 );
+ *_sauvFile << left;
+ map<string,int>::const_iterator nameNbIt = nameNbMap.begin();
+ for ( ; nameNbIt != nameNbMap.end(); nameNbIt++, fcount++ )
+ *_sauvFile << " " << setw(8) << nameNbIt->first;
+ fcount.stop();
+ *_sauvFile << right;
+
+ // write IDs of named objects in the pile
+ // * 8000 FORMAT(10I8)
+ nameNbIt = nameNbMap.begin();
+ for ( fcount.init(10); nameNbIt != nameNbMap.end(); nameNbIt++, fcount++ )
+ *_sauvFile << setw(8) << nameNbIt->second;
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Writes "PILE NUMERO 1"
+ */
+//================================================================================
+
+void SauvWriter::writeSubMeshes()
+{
+ int nbSauvObjects;
+ map<string,int> nameNbMap;
+ fillSubMeshes( nbSauvObjects, nameNbMap );
+
+ // * 800 FORMAT (' ENREGISTREMENT DE TYPE', I4)
+ *_sauvFile << " ENREGISTREMENT DE TYPE 2" << endl;
+ // * 801 FORMAT(' PILE NUMERO',I4,'NBRE OBJETS NOMMES',I8,'NBRE OBJETS',I8)
+ *_sauvFile << " PILE NUMERO 1NBRE OBJETS NOMMES" << setw(8) << nameNbMap.size() <<
+ "NBRE OBJETS" << setw(8) << nbSauvObjects <<endl;
+
+ writeNames( nameNbMap );
+
+ TFieldCounter fcount( *_sauvFile, 10 ); // 10 intergers per line
+
+ for ( size_t iSub = 0; iSub < _subs.size(); ++iSub )
+ {
+ SubMesh& sm = _subs[iSub];
+ if ( sm._nbSauvObjects < 1 ) continue;
+
+ // The first record of each sub-mesh writes
+ // - type of cells; zero means a compound object whose the 2nd record enumerates its components
+ // - number of components of a compound object
+ // - number of references; each reference means a "pointer" to this sub-mesh
+ // - number of nodes per cell
+ // - number of cells
+
+ if ( !sm._subs.empty() )
+ {
+ writeCompoundSubMesh(iSub);
+ }
+ else
+ {
+ // write each sub-type as a SAUV sub-mesh
+ MEDCouplingAutoRefCountObjectPtr< MEDCouplingMesh >
+ mesh = _fileMesh->getGenMeshAtLevel( sm._dimRelExt );
+ MEDCouplingAutoRefCountObjectPtr< MEDCouplingUMesh>
+ umesh = mesh->buildUnstructured();
+
+ for ( int iType=0; iType < sm.cellIDsByTypeSize(); ++iType )
+ {
+ const vector<int>& cellIDs = sm._cellIDsByType[iType];
+ if ( cellIDs.empty() ) continue;
+
+ INTERP_KERNEL::NormalizedCellType
+ cellType = INTERP_KERNEL::NormalizedCellType( iType );
+ const INTERP_KERNEL::CellModel &
+ cell = INTERP_KERNEL::CellModel::GetCellModel( cellType );
+ int castemType = SauvUtilities::med2gibiGeom( cellType );
+ unsigned nbElemNodes = cell.getNumberOfNodes();
+ unsigned nbElems = cellIDs.size();
+
+ *_sauvFile << setw(8) << castemType
+ << zeroI8
+ << zeroI8
+ << setw(8) << nbElemNodes
+ << setw(8) << nbElems << endl;
+
+ // write color of each element
+ // * 8000 FORMAT(10I8)
+ for ( size_t i = 0; i < nbElems; ++i, fcount++ ) *_sauvFile << zeroI8;
+ fcount.stop();
+
+ // write connectivity
+ // gibi IDs are in FORTRAN mode while MEDCoupling IDs are in C mode
+ if ( sm._dimRelExt == 1 ) // nodes
+ {
+ for ( size_t i = 0; i < nbElems; ++i, fcount++ )
+ *_sauvFile << setw(8) << ( cellIDs[i] + 1 );
+ }
+ else
+ {
+ // indices to transform MED connectivity to GIBI one
+ const int * toMedConn = getGibi2MedQuadraticInterlace( cellType );
+
+ vector< int > cellConn( nbElemNodes ), transformedConn( nbElemNodes );
+ for ( size_t i = 0; i < nbElems; ++i )
+ {
+ cellConn.clear();
+ umesh->getNodeIdsOfCell( cellIDs[i], cellConn );
+ if ( toMedConn )
+ {
+ for ( unsigned j = 0; j < nbElemNodes; ++j )
+ transformedConn[ j ] = cellConn[ toMedConn[ j ]];
+ cellConn.swap( transformedConn );
+ }
+ for ( unsigned j = 0; j < nbElemNodes; ++j, fcount++ )
+ *_sauvFile << setw(8) << ( cellConn[j] + 1 );
+ }
+ }
+ fcount.stop();
+
+ } // loop on cell types
+ } // not a compound object
+ } // loop on sub-meshes
+}
+
+//================================================================================
+/*!
+ * \brief Writes a sum-mesh composed of other sum-meshes
+ * This submesh corresponds to a med mesh or group composed of families
+ */
+//================================================================================
+
+void SauvWriter::writeCompoundSubMesh(int iSub)
+{
+ SubMesh& sm = _subs[iSub];
+ if ( sm._nbSauvObjects < 1 || sm._subs.empty()) return;
+
+ vector< int > subIDs;
+ for ( size_t i = 0; i < sm._subs.size(); ++i ) // loop on sub-meshes of families
+ for ( int j = 0; j < sm._subs[i]->_nbSauvObjects; ++j )
+ subIDs.push_back( sm._subs[i]->_id + j );
+
+ *_sauvFile << zeroI8
+ << setw(8) << subIDs.size()
+ << zeroI8
+ << zeroI8
+ << zeroI8 << endl;
+
+ TFieldCounter fcount( *_sauvFile, 10 ); // 10 intergers per line
+ for ( size_t i = 0; i < subIDs.size(); ++i, fcount++ )
+ *_sauvFile << setw(8) << subIDs[i];
+}
+
+//================================================================================
+/*!
+ * \brief Write piles relating to nodes
+ */
+//================================================================================
+
+void SauvWriter::writeNodes()
+{
+ MEDCouplingAutoRefCountObjectPtr< MEDCouplingMesh > mesh = _fileMesh->getGenMeshAtLevel( 1 );
+ MEDCouplingAutoRefCountObjectPtr< MEDCouplingUMesh > umesh = mesh->buildUnstructured();
+
+ // write the index connecting nodes with their coodrinates
+
+ const int nbNodes = umesh->getNumberOfNodes();
+ *_sauvFile << " ENREGISTREMENT DE TYPE 2" << endl
+ << " PILE NUMERO 32NBRE OBJETS NOMMES 0NBRE OBJETS" << setw(8) << nbNodes << endl;
+ *_sauvFile << setw(8) << nbNodes << endl;
+ //
+ TFieldCounter fcount( *_sauvFile, 10 );// * 8000 FORMAT(10I8)
+ for ( int i = 0; i < nbNodes; ++i, fcount++ )
+ *_sauvFile << setw(8) << i + 1;
+ fcount.stop();
+
+ // write coordinates and density of nodes
+
+ *_sauvFile << " ENREGISTREMENT DE TYPE 2" << endl;
+ *_sauvFile << " PILE NUMERO 33NBRE OBJETS NOMMES 0NBRE OBJETS 1" << endl;
+ //
+ const int dim = umesh->getSpaceDimension();
+ const int nbValues = nbNodes * ( dim + 1 );
+ *_sauvFile << setw(8) << nbValues << endl;
+
+ // * 8003 FORMAT(1P,3E22.14)
+ const char* density = " 0.00000000000000E+00";
+ fcount.init(3);
+ _sauvFile->precision(14);
+ _sauvFile->setf( ios_base::scientific, ios_base::floatfield );
+ _sauvFile->setf( ios_base::uppercase );
+ MEDCouplingAutoRefCountObjectPtr< DataArrayDouble> coordArray = umesh->getCoordinatesAndOwner();
+ const double precision = 1.e-99; // PAL12077
+ for ( int i = 0; i < nbNodes; ++i)
+ {
+ for ( int j = 0; j < dim; ++j, fcount++ )
+ {
+ double coo = coordArray->getIJ( i, j );
+ bool zero = ( -precision < coo && coo < precision );
+ *_sauvFile << setw(22) << ( zero ? 0.0 : coo );
+ }
+ *_sauvFile << density;
+ fcount++;
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Store correspondence between GIBI (short) and MED (long) names
+ *
+ * IMP 0020434: mapping GIBI names to MED names
+ * Store correspondence between GIBI and MED names as one PILE_STRINGS and one
+ * PILE_TABLES (in three tables: MED_MAIL, MED_CHAM and MED_COMP)
+ */
+//================================================================================
+
+void SauvWriter::writeLongNames()
+{
+ int nbTables =
+ 3 - _longNames[ LN_MAIL ].empty() - _longNames[ LN_CHAM ].empty() - _longNames[ LN_COMP ].empty();
+ if (nbTables == 0) return;
+
+ // ---------------------
+ // Write the TABLE pile
+ // ---------------------
+
+ *_sauvFile << " ENREGISTREMENT DE TYPE 2" << endl
+ << " PILE NUMERO 10NBRE OBJETS NOMMES" << setw(8) << nbTables
+ << "NBRE OBJETS" << setw(8) << nbTables << endl;
+ // table names
+ if (!_longNames[ LN_MAIL ].empty()) *_sauvFile << " MED_MAIL";
+ if (!_longNames[ LN_CHAM ].empty()) *_sauvFile << " MED_CHAM";
+ if (!_longNames[ LN_COMP ].empty()) *_sauvFile << " MED_COMP";
+ *_sauvFile << endl;
+ // table indices
+ for ( int i = 0; i < nbTables; ++i ) *_sauvFile << setw(8) << i+1;
+ *_sauvFile << endl;
+
+ string theWholeString; // concatenated long names
+ vector<int> theOffsets;
+ int iStr = 1;
+ TFieldCounter fcount (*_sauvFile, 10);
+
+ for ( int iTbl = 0; iTbl < LN_NB; ++iTbl )
+ {
+ vector<nameGIBItoMED>& longNames = _longNames[ iTbl ];
+ if ( longNames.empty() ) continue;
+ const bool isComp = ( iTbl == LN_COMP);
+
+ // to assure unique MED names
+ set<string> medUniqueNames;
+
+ *_sauvFile << setw(8) << longNames.size()*4 << endl; // Nb of table values
+
+ vector<nameGIBItoMED>::iterator itGIBItoMED = longNames.begin();
+ for (; itGIBItoMED != longNames.end(); itGIBItoMED++, iStr++)
+ {
+ // PILE of i-th key (med name)
+ *_sauvFile << setw(8) << PILE_STRINGS;
+ fcount++;
+ // ID of i-th key (med name)
+ *_sauvFile << setw(8) << iStr;
+ fcount++;
+ // PILE of i-th value (gibi name)
+ *_sauvFile << setw(8) << itGIBItoMED->gibi_pile;
+ fcount++;
+ // ID of i-th value (gibi name)
+ *_sauvFile << setw(8) << ( isComp ? ++iStr : itGIBItoMED->gibi_id );
+ fcount++;
+
+ // add a MED name to the string (while making it be unique for sub-meshes and fields)
+ string aMedName = itGIBItoMED->med_name;
+ if ( !isComp )
+ for (int ind = 1; !medUniqueNames.insert(aMedName).second; ++ind )
+ aMedName = itGIBItoMED->med_name + "_" + SauvUtilities::toString( ind );
+ theWholeString += aMedName;
+
+ // add an offset
+ theOffsets.push_back( theWholeString.size() );
+ if ( isComp )
+ {
+ theWholeString += itGIBItoMED->gibi_name;
+ theOffsets.push_back( theWholeString.size() );
+ }
+ }
+ fcount.stop();
+ }
+
+ // ----------------------
+ // Write the STRING pile
+ // ----------------------
+
+ const int nbNames = theOffsets.size();
+ *_sauvFile << " ENREGISTREMENT DE TYPE 2" << endl
+ << " PILE NUMERO 27NBRE OBJETS NOMMES" << zeroI8 << "NBRE OBJETS" << setw(8) << nbNames << endl
+ << setw(8) << theWholeString.length() << setw(8) << nbNames << endl;
+
+ // write the whole string
+ const int fixedLength = 71;
+ for ( string::size_type aPos = 0; aPos < theWholeString.length(); aPos += fixedLength)
+ *_sauvFile << setw(72) << theWholeString.substr(aPos, fixedLength) << endl;
+
+ // write the offsets
+ for ( size_t i = 0; i < theOffsets.size(); ++i, fcount++ )
+ *_sauvFile << setw(8) << theOffsets[i];
+}
+
+//================================================================================
+/*!
+ * \brief Write beginning of field record
+ */
+//================================================================================
+
+void SauvWriter::writeFieldNames( const bool isNodal,
+ std::map<std::string,int>& fldNamePrefixMap)
+{
+ vector< MEDCouplingAutoRefCountObjectPtr< MEDFileFieldMultiTS > >&
+ flds = isNodal ? _nodeFields : _cellFields;
+ map<string,int> nameNbMap;
+
+ for ( size_t iF = 0; iF < flds.size(); ++iF )
+ {
+ string name = addName( nameNbMap, fldNamePrefixMap, flds[iF]->getName(), iF+1 );
+ nameGIBItoMED aMEDName;
+ aMEDName.gibi_pile = isNodal ? PILE_NODES_FIELD : PILE_FIELD;
+ aMEDName.gibi_id = iF+1;
+ aMEDName.med_name = name;
+ _longNames[ LN_CHAM ].push_back(aMEDName);
+ }
+
+ *_sauvFile << " ENREGISTREMENT DE TYPE 2" << endl
+ << ( isNodal ? " PILE NUMERO 2" : " PILE NUMERO 39")
+ << "NBRE OBJETS NOMMES" << setw(8) << nameNbMap.size()
+ << "NBRE OBJETS" << setw(8) << flds.size() << endl;
+ writeNames( nameNbMap );
+}
+
+//================================================================================
+/*!
+ * \brief Make short names of field components
+ *
+ * IMP 0020434: mapping GIBI names to MED names
+ */
+//================================================================================
+
+void SauvWriter::makeCompNames(const string& fieldName,
+ const vector<string>& compInfo,
+ map<string, string>& mapMedToGibi)
+{
+ for ( size_t i = 0; i < compInfo.size(); ++i )
+ mapMedToGibi[compInfo[i]] = cleanName( compInfo[i] );
+
+ int compIndex = 1;
+ map<string, string>::iterator namesIt = mapMedToGibi.begin();
+ for (; namesIt != mapMedToGibi.end(); namesIt++)
+ {
+ string & compGibiName = (*namesIt).second;
+ if (compGibiName.size() > 4) {
+ // use new name in form "CXXX", where "XXX" is a number
+ do
+ {
+ compGibiName = SauvUtilities::toString( compIndex++ );
+ if ( compGibiName.size() < 3 )
+ compGibiName.insert( 0, 3 - compGibiName.size(), '0' );
+ compGibiName = "C" + compGibiName;
+ }
+ while (mapMedToGibi.count(compGibiName) > 0); // real component name could be CXXX
+ }
+
+ string compMedName = fieldName + "." + namesIt->first;
+ nameGIBItoMED aMEDName;
+ aMEDName.med_name = compMedName;
+ aMEDName.gibi_pile = PILE_STRINGS;
+ aMEDName.gibi_name = compGibiName;
+ _longNames[ LN_COMP ].push_back(aMEDName);
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Writes "PILE NUMERO 2": fields on nodes
+ */
+//================================================================================
+
+void SauvWriter::writeNodalFields(map<string,int>& fldNamePrefixMap)
+{
+ writeFieldNames( /*isNodal=*/true, fldNamePrefixMap );
+
+ TFieldCounter fcount (*_sauvFile, 10);
+
+ // EXAMPLE ( with no values )
+
+ // (1) 4 7 2 1
+ // (2) -88 0 3 -89 0 1 -90 0 2 -91
+ // (2) 0 1
+ // (3) FX FY FZ FZ FX FY FLX
+ // (4) 0 0 0 0 0 0 0
+ // (5) cree par muc pri
+ // (6)
+ // (7) 2
+ for ( size_t iF = 0; iF < _nodeFields.size(); ++iF )
+ {
+ // (1) write nb subcomponents, nb components(total)
+ vector< pair<int,int> > iters = _nodeFields[iF]->getIterations();
+ const vector<string>& compInfo = _nodeFields[iF]->getInfo();
+ const int nbSub = iters.size();
+ const int nbComp = compInfo.size();
+ const int totalNbComp = nbSub * nbComp;
+ *_sauvFile << setw(8) << nbSub
+ << setw(8) << totalNbComp
+ << setw(8) << -1 // IFOUR
+ << setw(8) << 0 << endl; // nb attributes
+
+ // (2) for each sub-component (iteration)
+ // write support, number of values and number of components
+ fcount.init(10);
+ vector< int > vals(3);
+ for ( int iIt = 0; iIt < iters.size(); ++iIt )
+ {
+ pair<int,int> it = iters[iIt];
+
+ vector<INTERP_KERNEL::NormalizedCellType> types;
+ vector< vector<TypeOfField> > typesF;
+ vector< vector<string> > pfls, locs;
+ vector< vector<const DataArrayDouble *> > valsVec;
+ valsVec = _nodeFields[iF]->getFieldSplitedByType( it.first, it.second, _fileMesh->getName(),
+ types, typesF, pfls, locs);
+ // believe that there can be only one type in a nodal field,
+ // so do not use a loop on types
+ if ( pfls[0][0].empty() ) pfls[0][0] = noProfileName( types[0] );
+ map< string, SubMesh* >::iterator pfl2Sub = _profile2Sub.find( pfls[0][0] );
+ if ( pfl2Sub == _profile2Sub.end() )
+ THROW_IK_EXCEPTION( "SauvWriter::writeNodalFields(): no sub-mesh for profile |"
+ << pfls[0][0] << "|");
+ vals[0] = -pfl2Sub->second->_id;
+ vals[1] = valsVec[0][0]->getNumberOfTuples() / nbComp;
+ vals[2] = compInfo.size();
+ for ( size_t i = 0; i < vals.size(); ++i, fcount++ )
+ *_sauvFile << setw(8) << vals[i];
+ }
+ fcount.stop();
+
+ // (3) Write names of components
+ map<string, string> mapMedToGibi;
+ makeCompNames( _nodeFields[iF]->getName(), compInfo, mapMedToGibi );
+ fcount.init(8);
+ *_sauvFile << left;
+ for ( int iIt = 0; iIt < iters.size(); ++iIt )
+ for ( size_t i = 0; i < compInfo.size(); ++i, fcount++ )
+ *_sauvFile << " " << setw(4) << mapMedToGibi[compInfo[i]];
+ *_sauvFile << right;
+ fcount.stop();
+
+ // (4) nb harmonics
+ fcount.init(10);
+ for ( size_t i = 0; i < totalNbComp; ++i, fcount++ )
+ *_sauvFile << " " << setw(8) << 0;
+ fcount.stop();
+
+ string description = _nodeFields[iF]->getName();
+ *_sauvFile << endl; // (5) TYPE
+ *_sauvFile << setw(72) << description.substr(0,71) << endl; // (6) TITRE
+ //*_sauvFile << endl; // (7) 0 attributes
+
+ // write values of each component
+ fcount.init( 3 ); // 3 values per a line
+ for ( int iIt = 0; iIt < iters.size(); ++iIt )
+ {
+ pair<int,int> it = iters[iIt];
+
+ vector<INTERP_KERNEL::NormalizedCellType> types;
+ vector< vector<TypeOfField> > typesF;
+ vector< vector<string> > pfls, locs;
+ vector< vector<const DataArrayDouble *> > valsVec;
+ valsVec = _nodeFields[iF]->getFieldSplitedByType( it.first, it.second, _fileMesh->getName(),
+ types, typesF, pfls, locs);
+ // believe that there can be only one type in a nodal field,
+ // so do not perform a loop on types
+ const DataArrayDouble* valsArray = valsVec[0][0];
+ for ( size_t j = 0; j < compInfo.size(); ++j )
+ {
+ for ( size_t i = 0; i < valsArray->getNumberOfTuples(); ++i, fcount++ )
+ *_sauvFile << setw(22) << valsArray->getIJ( i, j );
+ fcount.stop();
+ }
+ }
+ } // loop on fiels
+}
+
+//================================================================================
+/*!
+ * \brief Writes "PILE NUMERO 39": fields on cells
+ */
+//================================================================================
+
+void SauvWriter::writeElemFields(map<string,int>& fldNamePrefixMap)
+{
+ writeFieldNames( /*isNodal=*/false, fldNamePrefixMap );
+
+ TFieldCounter fcount (*_sauvFile, 10);
+
+ // REAL EXAMPLE
+
+ // (1) 1 2 6 16
+ // (2) CARACTERISTIQUES
+ // (3) -15 317773 4 0 0 0 -2 0 3
+ // (4) 317581
+ // (5) 0
+ // (6) 317767 317761 317755 317815
+ // (7) YOUN NU H SIGY
+ // (8) REAL*8 REAL*8 REAL*8 REAL*8
+ // (9) 1 1 0 0
+ // (10) 2.00000000000000E+05
+ // (9) 1 1 0 0
+ // (10) 3.30000000000000E-01
+ // (9) 1 1 0 0
+ // (10) 1.00000000000000E+04
+ // (9) 6 706 0 0
+ // (10) 1.00000000000000E+02 1.00000000000000E+02 1.00000000000000E+02
+ // (10) 1.00000000000000E+02 1.00000000000000E+02 1.00000000000000E+02
+ // (10) ...
+
+ for ( size_t iF = 0; iF < _cellFields.size(); ++iF )
+ {
+ // count nb of sub-components
+ int iSub, nbSub = 0;
+ vector< pair<int,int> > iters = _cellFields[iF]->getIterations();
+ for ( int iIt = 0; iIt < iters.size(); ++iIt )
+ {
+ pair<int,int> it = iters[iIt];
+
+ vector<INTERP_KERNEL::NormalizedCellType> types;
+ vector< vector<TypeOfField> > typesF;
+ vector< vector<string> > pfls, locs;
+ vector< vector<const DataArrayDouble *> > valsVec;
+ valsVec = _cellFields[iF]->getFieldSplitedByType( it.first, it.second, _fileMesh->getName(),
+ types, typesF, pfls, locs);
+ for ( size_t i = 0; i < valsVec.size(); ++i )
+ nbSub += valsVec[i].size();
+ }
+ // (1) write nb sub-components, title length
+ *_sauvFile << setw(8) << nbSub
+ << setw(8) << -1 // whatever
+ << setw(8) << 6 // whatever
+ << setw(8) << 72 << endl; // title length
+ // (2) title
+ string title = _cellFields[iF]->getName();
+ *_sauvFile << setw(72) << title.substr(0,71) << endl;
+ *_sauvFile << setw(72) << " " << endl;
+
+ // (3) support, nb components
+ vector<int> vals(9, 0);
+ const vector<string>& compInfo = _cellFields[iF]->getInfo();
+ vals[2] = compInfo.size();
+ fcount.init(10);
+ for ( int iIt = 0; iIt < iters.size(); ++iIt )
+ {
+ pair<int,int> it = iters[iIt];
+
+ vector<INTERP_KERNEL::NormalizedCellType> types;
+ vector< vector<TypeOfField> > typesF;
+ vector< vector<string> > pfls, locs;
+ vector< vector<const DataArrayDouble *> > valsVec;
+ valsVec = _cellFields[iF]->getFieldSplitedByType( it.first, it.second, _fileMesh->getName(),
+ types, typesF, pfls, locs);
+ for ( size_t iType = 0; iType < pfls.size(); ++iType )
+ for ( size_t iP = 0; iP < pfls[iType].size(); ++iP )
+ {
+ if ( pfls[iType][iP].empty() ) pfls[iType][iP] = noProfileName( types[iType] );
+ map< string, SubMesh* >::iterator pfl2Sub = _profile2Sub.find( pfls[iType][iP] );
+ if ( pfl2Sub == _profile2Sub.end() )
+ THROW_IK_EXCEPTION( "SauvWriter::writeElemFields(): no sub-mesh for profile |"
+ << pfls[iType][iP] << "|");
+ const int supportID = pfl2Sub->second->_id;
+ vals[0] = -supportID;
+
+ for ( size_t i = 0; i < vals.size(); ++i, fcount++ )
+ *_sauvFile << setw(8) << vals[ i ];
+ }
+ }
+ fcount.stop();
+
+ // (4) dummy strings
+ for ( fcount.init(4), iSub = 0; iSub < nbSub; ++iSub, fcount++ )
+ *_sauvFile << " ";
+ fcount.stop();
+
+ // (5) dummy strings
+ for ( fcount.init(8), iSub = 0; iSub < nbSub; ++iSub, fcount++ )
+ *_sauvFile << " ";
+ fcount.stop();
+
+ // loop on sub-components of a field, each of which refers to
+ // a certain support and has its own number of components
+ for ( int iIt = 0; iIt < iters.size(); ++iIt )
+ {
+ pair<int,int> it = iters[iIt];
+ writeElemTimeStamp( iF, it.first, it.second );
+ }
+ } // loop on cell fields
+}
+
+//================================================================================
+/*!
+ * \brief Write one elemental time stamp
+ */
+//================================================================================
+
+void SauvWriter::writeElemTimeStamp(int iF, int iter, int order)
+{
+ // (6) 317767 317761 317755 317815
+ // (7) YOUN NU H SIGY
+ // (8) REAL*8 REAL*8 REAL*8 REAL*8
+ // (9) 1 1 0 0
+ // (10) 2.00000000000000E+05
+ // (9) 1 1 0 0
+ // (10) 3.30000000000000E-01
+ // (9) 1 1 0 0
+ // (10) 1.00000000000000E+04
+ // (9) 6 706 0 0
+ // (10) 1.00000000000000E+02 1.00000000000000E+02 1.00000000000000E+02
+ // (10) 1.00000000000000E+02 1.00000000000000E+02 1.00000000000000E+02
+
+ TFieldCounter fcount (*_sauvFile, 10);
+
+ vector<INTERP_KERNEL::NormalizedCellType> types;
+ vector< vector<TypeOfField> > typesF;
+ vector< vector<string> > pfls, locs;
+ vector< vector<const DataArrayDouble *> > valsVec;
+ valsVec = _cellFields[iF]->getFieldSplitedByType( iter, order, _fileMesh->getName(),
+ types, typesF, pfls, locs);
+ for ( size_t iType = 0; iType < pfls.size(); ++iType )
+ for ( size_t iP = 0; iP < pfls[iType].size(); ++iP )
+ {
+ const vector<string>& compInfo = _cellFields[iF]->getInfo();
+
+ // (6) component addresses
+ int iComp = 0, nbComp = compInfo.size();
+ for ( fcount.init(10); iComp < nbComp; ++iComp, fcount++ )
+ *_sauvFile << setw(8) << 777; // a good number
+ fcount.stop();
+
+ // (7) component names
+ map<string, string> mapMedToGibi;
+ makeCompNames( _cellFields[iF]->getName(), compInfo, mapMedToGibi );
+ *_sauvFile << left;
+ for ( fcount.init(8), iComp = 0; iComp < nbComp; ++iComp, fcount++ )
+ *_sauvFile << " " << setw(8) << mapMedToGibi[compInfo[iComp]];
+ fcount.stop();
+
+ // (8) component types
+ for ( fcount.init(4), iComp = 0; iComp < nbComp; ++iComp, fcount++ )
+ *_sauvFile << " " << setw(17) << "REAL*8";
+ fcount.stop();
+ *_sauvFile << right;
+
+ // (9) nb values per element, nb of elements
+ int nbPntPerCell = 1;
+ if ( !locs[iType][iP].empty() )
+ {
+ int locID = _cellFields[iF]->getLocalizationId( locs[iType][iP].c_str() );
+ nbPntPerCell = _cellFields[iF]->getNbOfGaussPtPerCell( locID );
+ }
+ else if ( typesF[iType][iP] == ON_GAUSS_NE )
+ {
+ nbPntPerCell = INTERP_KERNEL::CellModel::GetCellModel(types[iType]).getNumberOfNodes();
+ }
+
+ // (10) values
+ const DataArrayDouble* valArray = valsVec[iType][iP];
+ for ( iComp = 0; iComp < nbComp; ++iComp )
+ {
+ *_sauvFile << setw(8) << nbPntPerCell
+ << setw(8) << valArray->getNbOfElems() / nbPntPerCell / nbComp
+ << setw(8) << 0
+ << setw(8) << 0
+ << endl;
+ fcount.init(3);
+ for ( size_t i = 0; i < valArray->getNumberOfTuples(); ++i, fcount++ )
+ *_sauvFile << setw(22) << valArray->getIJ( i, iComp );
+ fcount.stop();
+ }
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Write the last record of the SAUV file
+ */
+//================================================================================
+
+void SauvWriter::writeLastRecord()
+{
+ *_sauvFile << " ENREGISTREMENT DE TYPE 5" << endl;
+ *_sauvFile << "LABEL AUTOMATIQUE : 1" << endl;
+}
--- /dev/null
+// Copyright (C) 2007-2011 CEA/DEN, EDF R&D
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+// File : SauvWriter.hxx
+// Created : Wed Aug 24 11:18:49 2011
+// Author : Edward AGAPOV (eap)
+
+#ifndef __SauvWriter_HXX__
+#define __SauvWriter_HXX__
+
+#include "MEDLoaderDefines.hxx"
+#include "MEDCouplingRefCountObject.hxx"
+#include "NormalizedUnstructuredMesh.hxx"
+#include "SauvUtilities.hxx"
+#include "MEDCouplingAutoRefCountObjectPtr.hxx"
+
+#include <vector>
+#include <string>
+#include <map>
+
+namespace ParaMEDMEM
+{
+ class MEDFileData;
+ class MEDFileMesh;
+ class MEDFileFieldMultiTS;
+ class DataArrayInt;
+
+ /*!
+ * \brief Class to write a MEDFileData into a SAUVE format file
+ */
+ class MEDLOADER_EXPORT SauvWriter : public ParaMEDMEM::RefCountObject
+ {
+ public:
+ static SauvWriter * New();
+ void setMEDFileDS(const MEDFileData* medData, unsigned meshIndex = 0);
+ void write(const char* fileName);
+
+ private:
+
+ /*!
+ * \brief Class representing a GIBI sub-mesh (described in the pile 1 of the SAUVE file).
+ * It stands for a named med sub-mesh (family, etc) and contains either cell IDs or other sub-meshes.
+ */
+ struct SubMesh
+ {
+ std::vector<int> _cellIDsByType[ INTERP_KERNEL::NORM_MAXTYPE+1 ];
+ std::vector<SubMesh*> _subs;
+ std::string _name;
+ int _id;
+ int _nbSauvObjects;
+ int _dimRelExt;
+ int nbTypes() const;
+ static int cellIDsByTypeSize() { return INTERP_KERNEL::NORM_MAXTYPE+1; }
+ };
+ SubMesh* addSubMesh(const std::string& name, int dimRelExt);
+
+ void fillSubMeshes(int& nbSauvObjects, std::map<std::string,int>& nameNbMap);
+ void fillFamilySubMeshes();
+ void fillGroupSubMeshes();
+ void fillProfileSubMeshes();
+ int evaluateNbProfileSubMeshes() const;
+ void makeCompNames(const std::string& fieldName,
+ const std::vector<std::string>& compInfo,
+ std::map<std::string, std::string>& compMedToGibi );
+ void makeProfileIDs( SubMesh* sm,
+ INTERP_KERNEL::NormalizedCellType type,
+ const DataArrayInt* profile );
+ void writeFileHead();
+ void writeSubMeshes();
+ void writeCompoundSubMesh(int iSub);
+ void writeNodes();
+ void writeLongNames();
+ void writeNodalFields(std::map<std::string,int>& fldNamePrefixMap);
+ void writeElemFields(std::map<std::string,int>& fldNamePrefixMap);
+ void writeFieldNames(const bool isNodal, std::map<std::string,int>& fldNamePrefixMap);
+ void writeElemTimeStamp(int iF, int iter, int order);
+ void writeLastRecord();
+ void writeNames( const std::map<std::string,int>& nameNbMap );
+
+ private:
+
+ MEDCouplingAutoRefCountObjectPtr< MEDFileMesh > _fileMesh;
+ std::vector< MEDCouplingAutoRefCountObjectPtr< MEDFileFieldMultiTS > > _nodeFields;
+ std::vector< MEDCouplingAutoRefCountObjectPtr< MEDFileFieldMultiTS > > _cellFields;
+
+ std::vector<SubMesh> _subs;
+ std::map< int, SubMesh* > _famIDs2Sub;
+ std::map< std::string, SubMesh* > _profile2Sub;
+ enum
+ {
+ LN_MAIL=0, LN_CHAM, LN_COMP, LN_NB
+ };
+ std::vector<SauvUtilities::nameGIBItoMED> _longNames[ LN_NB ];
+
+ std::fstream* _sauvFile;
+ };
+}
+
+#endif
#include "MEDFileField.hxx"
#include "MEDFileData.hxx"
#include "MEDLoaderTypemaps.i"
+#include "SauvReader.hxx"
+#include "SauvWriter.hxx"
using namespace ParaMEDMEM;
%}
%newobject ParaMEDMEM::MEDFileData::getMeshes;
%newobject ParaMEDMEM::MEDFileData::getFields;
+%newobject ParaMEDMEM::SauvWriter::New;
+%newobject ParaMEDMEM::SauvReader::New;
+%newobject ParaMEDMEM::SauvReader::loadInMEDFileDS;
+
%feature("unref") MEDFileMesh "$this->decrRef();"
%feature("unref") MEDFileUMesh "$this->decrRef();"
%feature("unref") MEDFileCMesh "$this->decrRef();"
%feature("unref") MEDFileFieldMultiTS "$this->decrRef();"
%feature("unref") MEDFileFields "$this->decrRef();"
%feature("unref") MEDFileData "$this->decrRef();"
+%feature("unref") SauvReader "$this->decrRef();"
+%feature("unref") SauvWriter "$this->decrRef();"
class MEDLoader
{
int getNumberOfTS() const;
std::string getName() const;
std::string getMeshName() const throw(INTERP_KERNEL::Exception);
+ const std::vector<std::string>& getInfo() const;
%extend
{
PyObject *getIterations() const
}
}
};
+
+ class SauvReader : public RefCountObject
+ {
+ public:
+ static SauvReader* New(const char *fileName) throw(INTERP_KERNEL::Exception);
+ MEDFileData * loadInMEDFileDS() throw(INTERP_KERNEL::Exception);
+ };
+
+ class SauvWriter : public RefCountObject
+ {
+ public:
+ static SauvWriter * New();
+ void setMEDFileDS(const MEDFileData* medData, unsigned meshIndex = 0);
+ void write(const char* fileName);
+ };
+
}
--- /dev/null
+# -*- coding: iso-8859-1 -*-
+# Copyright (C) 2007-2011 CEA/DEN, EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+
+from MEDLoader import *
+import unittest, os
+from MEDLoaderDataForTest import MEDLoaderDataForTest
+
+class SauvLoaderTest(unittest.TestCase):
+
+ def testSauv2Med(self):
+ # get a file containing all types of readable piles
+ self.assertTrue( os.getenv("MED_ROOT_DIR") )
+ sauvFile = os.path.join( os.getenv("MED_ROOT_DIR"), "share","salome",
+ "resources","med","allPillesTest.sauv")
+ self.assertTrue( os.access( sauvFile, os.F_OK))
+
+ # read SAUV and write MED
+ medFile = "SauvLoaderTest.med"
+ sr=SauvReader.New(sauvFile);
+ d2=sr.loadInMEDFileDS();
+ d2.write(medFile,0);
+
+ # check
+ self.assertEqual(1,d2.getNumberOfMeshes())
+ self.assertEqual(8+97,d2.getNumberOfFields())
+ mm = d2.getMeshes()
+ m = mm.getMeshAtPos(0)
+ self.assertEqual(17,len(m.getGroupsNames()))
+
+ os.remove( medFile )
+ pass
+
+ def testMed2Sauv(self):
+ # read pointe.med
+ self.assertTrue( os.getenv("MED_ROOT_DIR") )
+ medFile = os.path.join( os.getenv("MED_ROOT_DIR"), "share","salome",
+ "resources","med","pointe.med")
+ self.assertTrue( os.access( medFile, os.F_OK))
+ pointeMed = MEDFileData.New( medFile )
+
+ # add 3 faces to pointeMed
+ pointeMedMesh = pointeMed.getMeshes().getMeshAtPos(0)
+ pointeM1D = MEDCouplingUMesh.New()
+ pointeM1D.setCoords( pointeMedMesh.getCoords() )
+ pointeM1D.setMeshDimension( 2 )
+ pointeM1D.allocateCells( 3 )
+ pointeM1D.insertNextCell( NORM_TRI3, 3, [0,1,2])
+ pointeM1D.insertNextCell( NORM_TRI3, 3, [0,1,3])
+ pointeM1D.insertNextCell( NORM_QUAD4, 4, [10,11,12,13])
+ pointeM1D.finishInsertingCells()
+ pointeMedMesh.setMeshAtLevel( -1, pointeM1D )
+ pointeMed.getMeshes().setMeshAtPos( 0, pointeMedMesh )
+
+ # add a field on 2 faces to pointeMed
+ ff1=MEDFileFieldMultiTS.New()
+ f1=MEDCouplingFieldDouble.New(ON_GAUSS_NE,ONE_TIME)
+ #f1.setMesh( pointeM1D )
+ f1.setName("Field on 2 faces")
+ d=DataArrayDouble.New()
+ d.alloc(3+4,2)
+ d.setInfoOnComponent(0,"sigX [MPa]")
+ d.setInfoOnComponent(1,"sigY [GPa]")
+ d.setValues([311,312,321,322,331,332,411,412,421,422,431,432,441,442],3+4,2)
+ f1.setArray(d)
+ da=DataArrayInt.New()
+ da.setValues([0,2],2,1)
+ da.setName("sup2")
+ ff1.appendFieldProfile(f1,pointeMedMesh,-1,da)
+ pointeMed.getFields().pushField( ff1 )
+
+ #remove fieldnodeint
+ pointeFields = pointeMed.getFields()
+ for i in range( pointeFields.getNumberOfFields() ):
+ if pointeFields.getFieldAtPos(i).getName() == "fieldnodeint":
+ pointeFields.destroyFieldAtPos( i )
+ break
+
+ # write pointeMed to SAUV
+ sauvFile = "SauvLoaderTest.sauv"
+ sw=SauvWriter.New();
+ sw.setMEDFileDS(pointeMed);
+ sw.write(sauvFile);
+
+ # read SAUV and check
+ sr=SauvReader.New(sauvFile);
+ d2=sr.loadInMEDFileDS();
+ self.assertEqual(1,d2.getNumberOfMeshes())
+ self.assertEqual(4,d2.getNumberOfFields())
+ m = d2.getMeshes().getMeshAtPos(0)
+ self.assertEqual("maa1",m.getName())
+ self.assertEqual(5,len(m.getGroupsNames()))
+ self.assertEqual(3,m.getMeshDimension())
+ groups = m.getGroupsNames()
+ self.assertTrue( "groupe1" in groups )
+ self.assertTrue( "groupe2" in groups )
+ self.assertTrue( "groupe3" in groups )
+ self.assertTrue( "groupe4" in groups )
+ self.assertTrue( "groupe5" in groups )
+ self.assertEqual(16,m.getSizeAtLevel(0))
+ um0 = m.getGenMeshAtLevel(0)
+ self.assertEqual(12, um0.getNumberOfCellsWithType( NORM_TETRA4 ))
+ self.assertEqual(2, um0.getNumberOfCellsWithType( NORM_PYRA5 ))
+ self.assertEqual(2, um0.getNumberOfCellsWithType( NORM_HEXA8 ))
+ um1 = m.getGenMeshAtLevel(-1)
+ self.assertEqual(2, um1.getNumberOfCellsWithType( NORM_TRI3 ))
+ pointeUM0 = pointeMedMesh.getGenMeshAtLevel(0)
+ self.assertTrue(m.getCoords().isEqualWithoutConsideringStr(pointeMedMesh.getCoords(),1e-12))
+ self.assertEqual( um0.getMeasureField(0).accumulate(0),
+ pointeUM0.getMeasureField(0).accumulate(0),1e-12)
+ # check fields
+ # fieldnodedouble
+ fieldnodedoubleTS1 = pointeMed.getFields().getFieldWithName("fieldnodedouble")
+ fieldnodedoubleTS2 = d2.getFields().getFieldWithName("fieldnodedouble")
+ self.assertEqual( fieldnodedoubleTS1.getInfo(), fieldnodedoubleTS2.getInfo())
+ self.assertEqual( fieldnodedoubleTS1.getNumberOfTS(), fieldnodedoubleTS2.getNumberOfTS())
+ io1 = fieldnodedoubleTS1.getIterations()
+ io2 = fieldnodedoubleTS2.getIterations()
+ for i in range(fieldnodedoubleTS1.getNumberOfTS() ):
+ fnd1 = fieldnodedoubleTS1.getFieldOnMeshAtLevel(ON_NODES, io1[i][0],io1[i][1],pointeUM0)
+ fnd2 = fieldnodedoubleTS2.getFieldOnMeshAtLevel(ON_NODES, io2[i][0],io2[i][1],um0)
+ self.assertTrue( fnd1.getArray().isEqual( fnd2.getArray(), 1e-12 ))
+ # fieldcelldoublevector
+ fieldnodedoubleTS1 = pointeMed.getFields().getFieldWithName("fieldcelldoublevector")
+ fieldnodedoubleTS2 = d2.getFields().getFieldWithName("fieldcelldoublevector")
+ self.assertEqual( fieldnodedoubleTS1.getInfo(), fieldnodedoubleTS2.getInfo())
+ self.assertEqual( fieldnodedoubleTS1.getNumberOfTS(), fieldnodedoubleTS2.getNumberOfTS())
+ io1 = fieldnodedoubleTS1.getIterations()
+ io2 = fieldnodedoubleTS2.getIterations()
+ for i in range(fieldnodedoubleTS1.getNumberOfTS() ):
+ fnd1 = fieldnodedoubleTS1.getFieldOnMeshAtLevel(ON_CELLS, io1[i][0],io1[i][1],pointeUM0)
+ fnd2 = fieldnodedoubleTS2.getFieldOnMeshAtLevel(ON_CELLS, io2[i][0],io2[i][1],um0)
+ self.assertAlmostEqual( fnd1.accumulate(0), fnd2.accumulate(0) )
+ self.assertAlmostEqual( fnd1.accumulate(1), fnd2.accumulate(1) )
+ self.assertAlmostEqual( fnd1.accumulate(2), fnd2.accumulate(2) )
+ # Field on 2 faces
+ fieldOnFaces = d2.getFields().getFieldWithName(f1.getName())
+ io1 = fieldOnFaces.getIterations()
+ fof = fieldOnFaces.getFieldOnMeshAtLevel(f1.getTypeOfField(),io1[i][0],io1[i][1],um1)
+ self.assertTrue( d.isEqual( fof.getArray(), 1e-12 ))
+
+ os.remove( sauvFile )
+ pass
+ pass
+
+unittest.main()
include $(top_srcdir)/adm_local/unix/make_common_starter.am
-bin_PROGRAMS= TestMEDLoader
+bin_PROGRAMS= TestMEDLoader TestSauvLoader
TestMEDLoader_CPPFLAGS=@CPPUNIT_INCLUDES@ @PTHREAD_CFLAGS@ -I$(srcdir)/.. -I$(srcdir)/../../INTERP_KERNEL/Bases -I$(srcdir)/../../INTERP_KERNELTest -I$(srcdir)/../../INTERP_KERNEL -I$(srcdir)/../../INTERP_KERNEL/Geometric2D -I$(srcdir)/../../MEDCoupling
+TestSauvLoader_CPPFLAGS=$(TestMEDLoader_CPPFLAGS) $(MED2_INCLUDES) $(HDF5_INCLUDES)
TestMEDLoader_LDFLAGS = @CPPUNIT_LIBS@ ../libmedloader.la ../../MEDCoupling/libmedcoupling.la ../../INTERP_KERNEL/libinterpkernel.la
+TestSauvLoader_LDFLAGS=$(TestMEDLoader_LDFLAGS)
dist_TestMEDLoader_SOURCES = TestMEDLoader.cxx MEDLoaderTest.cxx MEDLoaderTest.hxx
+dist_TestSauvLoader_SOURCES = TestSauvLoader.cxx SauvLoaderTest.cxx SauvLoaderTest.hxx
-UNIT_TEST_PROG = TestMEDLoader
+UNIT_TEST_PROG = TestMEDLoader TestSauvLoader
check : tests
CLEANFILES = \
UnitTestsResult
-
+
clean-local:
rm -rf *.med
-
\ No newline at end of file
--- /dev/null
+// Copyright (C) 2007-2011 CEA/DEN, EDF R&D
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+
+#include "SauvLoaderTest.hxx"
+
+#include "SauvReader.hxx"
+#include "SauvWriter.hxx"
+#include "MEDFileData.hxx"
+#include "MEDCouplingFieldDouble.hxx"
+#include "MEDCouplingMemArray.hxx"
+
+#ifdef WNT
+#include <windows.h>
+#endif
+
+#include <vector>
+#include <string>
+
+using namespace ParaMEDMEM;
+
+void SauvLoaderTest::testSauv2Med()
+{
+ // read a file containing all types of readable piles
+ std::string file = getResourceFile("allPillesTest.sauv");
+ MEDCouplingAutoRefCountObjectPtr<SauvReader> sr=SauvReader::New(file.c_str());
+ MEDCouplingAutoRefCountObjectPtr<MEDFileData> d2=sr->loadInMEDFileDS();
+ // write MED
+ d2->write("allPillesTest.med",0);
+ // check
+ CPPUNIT_ASSERT_EQUAL(1,d2->getNumberOfMeshes());
+ CPPUNIT_ASSERT_EQUAL(8+97,d2->getNumberOfFields());
+ MEDFileMesh * m = d2->getMeshes()->getMeshAtPos(0);
+ CPPUNIT_ASSERT_EQUAL(17,int(m->getGroupsNames().size()));
+}
+
+void SauvLoaderTest::testMed2Sauv()
+{
+ // read pointe.med
+ std::string file = getResourceFile("pointe.med");
+ MEDCouplingAutoRefCountObjectPtr<MEDFileData> pointeMed=MEDFileData::New(file.c_str());
+
+ // add 3 faces to pointeMed
+ MEDFileUMesh* pointeMedMesh = static_cast<MEDFileUMesh*>(pointeMed->getMeshes()->getMeshAtPos(0));
+ MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> pointeM1D = MEDCouplingUMesh::New();
+ MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> coords = pointeMedMesh->getCoords();
+ pointeM1D->setCoords( coords );
+ pointeM1D->setMeshDimension( 2 );
+ pointeM1D->allocateCells( 3 );
+ int conn[]=
+ {
+ 0,1,2, 0,1,3, 10,11,12,13
+ };
+ pointeM1D->insertNextCell( INTERP_KERNEL::NORM_TRI3, 3, conn);
+ pointeM1D->insertNextCell( INTERP_KERNEL::NORM_TRI3, 3, conn+3);
+ pointeM1D->insertNextCell( INTERP_KERNEL::NORM_QUAD4, 4, conn+6);
+ pointeM1D->finishInsertingCells();
+ pointeMedMesh->setMeshAtLevel( -1, pointeM1D );
+ pointeMed->getMeshes()->setMeshAtPos( 0, pointeMedMesh );
+
+ // add a field on 2 faces to pointeMed
+ MEDCouplingAutoRefCountObjectPtr<MEDFileFieldMultiTS> ff1=MEDFileFieldMultiTS::New();
+ MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> f1=MEDCouplingFieldDouble::New(ON_GAUSS_NE,ONE_TIME);
+ f1->setMesh( pointeM1D );
+ f1->setName("Field on 2 faces");
+ MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> d=DataArrayDouble::New();
+ d->alloc(3+4,2);
+ d->setInfoOnComponent(0,"sigX [MPa]");
+ d->setInfoOnComponent(1,"sigY [GPa]");
+ double vals[2*(3+4)] =
+ {
+ 311,312,321,322,331,332,411,412,421,422,431,432,441,442
+ };
+ std::copy(vals,vals+d->getNbOfElems(),d->getPointer());
+ f1->setArray(d);
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> da=DataArrayInt::New();
+ int ids[] =
+ {
+ 0,2
+ };
+ da->alloc(2,1);
+ std::copy(ids,ids+da->getNbOfElems(),da->getPointer());
+ da->setName("sup2");
+ ff1->appendFieldProfile(f1,pointeMedMesh,-1,da);
+ pointeMed->getFields()->pushField( ff1 );
+
+ // remove "fieldnodeint"
+ MEDFileFields* pointeFields = pointeMed->getFields();
+ for ( int i = 0; i < pointeFields->getNumberOfFields(); ++i )
+ {
+ MEDCouplingAutoRefCountObjectPtr<MEDFileFieldMultiTS> ts = pointeFields->getFieldAtPos(i);
+ if ( std::string("fieldnodeint") == ts->getName())
+ {
+ pointeFields->destroyFieldAtPos( i );
+ break;
+ }
+ }
+ // write pointeMed to SAUV
+ const char* sauvFile = "pointe.sauv";
+ MEDCouplingAutoRefCountObjectPtr<SauvWriter> sw=SauvWriter::New();
+ sw->setMEDFileDS(pointeMed);
+ sw->write(sauvFile);
+
+ // read SAUV and check
+ MEDCouplingAutoRefCountObjectPtr<SauvReader> sr=SauvReader::New(sauvFile);
+ MEDCouplingAutoRefCountObjectPtr<MEDFileData> d2=sr->loadInMEDFileDS();
+ CPPUNIT_ASSERT_EQUAL(1,d2->getNumberOfMeshes());
+ CPPUNIT_ASSERT_EQUAL(4,d2->getNumberOfFields());
+ MEDFileUMesh * m = static_cast<MEDFileUMesh*>( d2->getMeshes()->getMeshAtPos(0) );
+ CPPUNIT_ASSERT_EQUAL(std::string("maa1"),std::string(m->getName() ));
+ CPPUNIT_ASSERT_EQUAL(3,m->getMeshDimension());
+ std::vector<std::string > groups = m->getGroupsNames();
+ CPPUNIT_ASSERT_EQUAL(5,(int)groups.size());
+ CPPUNIT_ASSERT( std::find(groups.begin(),groups.end(),"groupe1") != groups.end() );
+ CPPUNIT_ASSERT( std::find(groups.begin(),groups.end(),"groupe2") != groups.end() );
+ CPPUNIT_ASSERT( std::find(groups.begin(),groups.end(),"groupe3") != groups.end() );
+ CPPUNIT_ASSERT( std::find(groups.begin(),groups.end(),"groupe4") != groups.end() );
+ CPPUNIT_ASSERT( std::find(groups.begin(),groups.end(),"groupe5") != groups.end() );
+ CPPUNIT_ASSERT_EQUAL(16,m->getSizeAtLevel(0));
+ MEDCouplingAutoRefCountObjectPtr<MEDCouplingMesh> um0 = m->getGenMeshAtLevel(0);
+ CPPUNIT_ASSERT_EQUAL(12, um0->getNumberOfCellsWithType( INTERP_KERNEL::NORM_TETRA4 ));
+ CPPUNIT_ASSERT_EQUAL(2, um0->getNumberOfCellsWithType( INTERP_KERNEL::NORM_PYRA5 ));
+ CPPUNIT_ASSERT_EQUAL(2, um0->getNumberOfCellsWithType( INTERP_KERNEL::NORM_HEXA8 ));
+ MEDCouplingAutoRefCountObjectPtr<MEDCouplingMesh> um1 = m->getGenMeshAtLevel(-1);
+ CPPUNIT_ASSERT_EQUAL(2, um1->getNumberOfCellsWithType( INTERP_KERNEL::NORM_TRI3 ));
+ MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> pointeUM0 =
+ static_cast<MEDCouplingUMesh*>( pointeMedMesh->getGenMeshAtLevel(0));
+ MEDCouplingAutoRefCountObjectPtr< DataArrayDouble > coo = m->getCoords();
+ MEDCouplingAutoRefCountObjectPtr< DataArrayDouble > pointeCoo = pointeMedMesh->getCoords();
+ CPPUNIT_ASSERT(coo->isEqualWithoutConsideringStr(*pointeCoo,1e-12));
+ MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> vol = um0->getMeasureField(0);
+ MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> pointeVol = pointeUM0->getMeasureField(0);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL( vol->accumulate(0), pointeVol->accumulate(0),1e-12);
+ // check fields
+ // fieldnodedouble
+ MEDCouplingAutoRefCountObjectPtr<MEDFileFieldMultiTS> fieldnodedoubleTS1 =
+ pointeMed->getFields()->getFieldWithName("fieldnodedouble");
+ MEDCouplingAutoRefCountObjectPtr<MEDFileFieldMultiTS> fieldnodedoubleTS2 =
+ d2->getFields()->getFieldWithName("fieldnodedouble");
+ CPPUNIT_ASSERT_EQUAL( fieldnodedoubleTS1->getInfo().size(), fieldnodedoubleTS2->getInfo().size());
+ for ( size_t i = 0; i < fieldnodedoubleTS1->getInfo().size(); ++i )
+ CPPUNIT_ASSERT_EQUAL( fieldnodedoubleTS1->getInfo()[i], fieldnodedoubleTS2->getInfo()[i]);
+ CPPUNIT_ASSERT_EQUAL( fieldnodedoubleTS1->getNumberOfTS(), fieldnodedoubleTS2->getNumberOfTS());
+ std::vector< std::pair<int,int> > io1 = fieldnodedoubleTS1->getIterations();
+ std::vector< std::pair<int,int> > io2 = fieldnodedoubleTS2->getIterations();
+ for ( int i =0; i < fieldnodedoubleTS1->getNumberOfTS(); ++i )
+ {
+ MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> fnd1 =
+ fieldnodedoubleTS1->getFieldOnMeshAtLevel(ON_NODES, io1[i].first,io1[i].second,pointeUM0);
+ MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> fnd2 =
+ fieldnodedoubleTS2->getFieldOnMeshAtLevel(ON_NODES, io2[i].first,io2[i].second,um0);
+ CPPUNIT_ASSERT( fnd1->getArray()->isEqual( *fnd2->getArray(), 1e-12 ));
+ }
+ // fieldcelldoublevector
+ MEDCouplingAutoRefCountObjectPtr<MEDFileFieldMultiTS> fieldcelldoublevectorTS1 =
+ pointeMed->getFields()->getFieldWithName("fieldcelldoublevector");
+ MEDCouplingAutoRefCountObjectPtr<MEDFileFieldMultiTS> fieldcelldoublevectorTS2 =
+ d2->getFields()->getFieldWithName("fieldcelldoublevector");
+ CPPUNIT_ASSERT_EQUAL( fieldcelldoublevectorTS1->getInfo().size(), fieldcelldoublevectorTS2->getInfo().size());
+ for ( size_t i = 0; i < fieldcelldoublevectorTS1->getInfo().size(); ++i )
+ CPPUNIT_ASSERT_EQUAL( fieldcelldoublevectorTS1->getInfo()[i], fieldcelldoublevectorTS2->getInfo()[i]);
+ CPPUNIT_ASSERT_EQUAL( fieldcelldoublevectorTS1->getNumberOfTS(), fieldcelldoublevectorTS2->getNumberOfTS());
+ io1 = fieldcelldoublevectorTS1->getIterations();
+ io2 = fieldcelldoublevectorTS2->getIterations();
+ for ( int i =0; i < fieldcelldoublevectorTS1->getNumberOfTS(); ++i )
+ {
+ MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> fnd1 =
+ fieldcelldoublevectorTS1->getFieldOnMeshAtLevel(ON_CELLS, io1[i].first,io1[i].second,pointeUM0);
+ MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> fnd2 =
+ fieldcelldoublevectorTS2->getFieldOnMeshAtLevel(ON_CELLS, io2[i].first,io2[i].second,um0);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL( fnd1->accumulate(0), fnd2->accumulate(0), 1e-12 );
+ CPPUNIT_ASSERT_DOUBLES_EQUAL( fnd1->accumulate(1), fnd2->accumulate(1), 1e-12 );
+ CPPUNIT_ASSERT_DOUBLES_EQUAL( fnd1->accumulate(2), fnd2->accumulate(2), 1e-12 );
+ }
+ // "Field on 2 faces"
+ MEDCouplingAutoRefCountObjectPtr<MEDFileFieldMultiTS> fieldOnFaces =
+ d2->getFields()->getFieldWithName(f1->getName());
+ io1 = fieldOnFaces->getIterations();
+ MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> fof =
+ fieldOnFaces->getFieldOnMeshAtLevel(f1->getTypeOfField(),io1[0].first,io1[0].second,um1);
+ CPPUNIT_ASSERT( d->isEqual( *fof->getArray(), 1e-12 ));
+}
+
+void SauvLoaderTest::tearDown()
+{
+ const int nbFilesToRemove = 2;
+ const char* fileToRemove[nbFilesToRemove] = { "allPillesTest.med", "pointe.sauv" };
+ for ( int i = 0; i < nbFilesToRemove; ++i )
+ {
+#ifdef WNT
+ if (GetFileAttributes(fileToRemove[i]) != INVALID_FILE_ATTRIBUTES)
+#else
+ if (access(fileToRemove[i], F_OK) == 0)
+#endif
+ remove(fileToRemove[i]);
+ }
+}
+
+std::string SauvLoaderTest::getResourceFile( const std::string& filename )
+{
+ std::string resourceFile = "";
+
+ if ( getenv("top_srcdir") ) {
+ // we are in 'make check' step
+ resourceFile = getenv("top_srcdir");
+ resourceFile += "/resources/";
+ }
+ else if ( getenv("MED_ROOT_DIR") ) {
+ // use MED_ROOT_DIR env.var
+ resourceFile = getenv("MED_ROOT_DIR");
+ resourceFile += "/share/salome/resources/med/";
+ }
+ resourceFile += filename;
+#ifdef WNT
+ std::string fixedpath = resourceFile;
+ for ( int i=0; i < path.size(); ++i )
+ if (path[i] == '/')
+ fixedpath[i] = '\\';
+ return fixedpath;
+#endif
+ return resourceFile;
+}
--- /dev/null
+// Copyright (C) 2007-2011 CEA/DEN, EDF R&D
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+
+#ifndef __SauvLoaderTest_HXX__
+#define __SauvLoaderTest_HXX__
+
+#include <cppunit/extensions/HelperMacros.h>
+
+#include "MEDLoaderTest.hxx"
+
+namespace ParaMEDMEM
+{
+ class SauvLoaderTest : public CppUnit::TestFixture
+ {
+ CPPUNIT_TEST_SUITE(SauvLoaderTest);
+ CPPUNIT_TEST( testSauv2Med );
+ CPPUNIT_TEST( testMed2Sauv );
+ CPPUNIT_TEST_SUITE_END();
+ public:
+ void testSauv2Med();
+ void testMed2Sauv();
+
+ public:
+ void tearDown();
+ static std::string getResourceFile( const std::string& filename );
+ };
+}
+
+#endif
--- /dev/null
+// Copyright (C) 2007-2011 CEA/DEN, EDF R&D
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+
+#include "CppUnitTest.hxx"
+#include "SauvLoaderTest.hxx"
+
+CPPUNIT_TEST_SUITE_REGISTRATION( ParaMEDMEM::SauvLoaderTest );
+
+#include "BasicMainTest.hxx"
