From a27acf81cb5c78604cd0eb3bd3d43bc7c6cf7426 Mon Sep 17 00:00:00 2001 From: eap Date: Thu, 1 Sep 2011 10:08:15 +0000 Subject: [PATCH] 0021314: [CEA 490] Driver GIBI with medloader / medcoupling --- src/MEDLoader/Makefile.am | 6 +- src/MEDLoader/SauvMedConvertor.cxx | 2437 +++++++++++++++++++++++++ src/MEDLoader/SauvMedConvertor.hxx | 383 ++++ src/MEDLoader/SauvReader.cxx | 1149 ++++++++++++ src/MEDLoader/SauvReader.hxx | 102 ++ src/MEDLoader/SauvUtilities.hxx | 117 ++ src/MEDLoader/SauvWriter.cxx | 1323 ++++++++++++++ src/MEDLoader/SauvWriter.hxx | 114 ++ src/MEDLoader/Swig/MEDLoader.i | 25 + src/MEDLoader/Swig/SauvLoaderTest.py | 162 ++ src/MEDLoader/Test/Makefile.am | 10 +- src/MEDLoader/Test/SauvLoaderTest.cxx | 237 +++ src/MEDLoader/Test/SauvLoaderTest.hxx | 45 + src/MEDLoader/Test/TestSauvLoader.cxx | 25 + 14 files changed, 6129 insertions(+), 6 deletions(-) create mode 100644 src/MEDLoader/SauvMedConvertor.cxx create mode 100644 src/MEDLoader/SauvMedConvertor.hxx create mode 100644 src/MEDLoader/SauvReader.cxx create mode 100644 src/MEDLoader/SauvReader.hxx create mode 100644 src/MEDLoader/SauvUtilities.hxx create mode 100644 src/MEDLoader/SauvWriter.cxx create mode 100644 src/MEDLoader/SauvWriter.hxx create mode 100644 src/MEDLoader/Swig/SauvLoaderTest.py create mode 100644 src/MEDLoader/Test/SauvLoaderTest.cxx create mode 100644 src/MEDLoader/Test/SauvLoaderTest.hxx create mode 100644 src/MEDLoader/Test/TestSauvLoader.cxx diff --git a/src/MEDLoader/Makefile.am b/src/MEDLoader/Makefile.am index 3c1a8e523..99be5b856 100755 --- a/src/MEDLoader/Makefile.am +++ b/src/MEDLoader/Makefile.am @@ -37,12 +37,14 @@ lib_LTLIBRARIES = libmedloader.la 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 \ diff --git a/src/MEDLoader/SauvMedConvertor.cxx b/src/MEDLoader/SauvMedConvertor.cxx new file mode 100644 index 000000000..d19ec9d45 --- /dev/null +++ b/src/MEDLoader/SauvMedConvertor.cxx @@ -0,0 +1,2437 @@ +// 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 +#include +#include +#include +#include + +#include +#include +#include + +#ifdef WNT +#include +#endif + +#ifndef WNT +#define HAS_XDR +#endif + +#ifdef HAS_XDR +#include +#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 > & 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 > & 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 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 treatedGroups; + + list::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::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 & fields = isNodal ? _nodeFields : _cellFields; + for ( size_t ifi = 0; ifi < fields.size() && !name_found; ifi++) + { + if (medName.find( fields[ifi]->_name + "." ) == 0 ) + { + vector& 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 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::const_iterator elemIt, elemEnd; + vector< pair > swapVec; + + // ------------------------------------ + // fix connectivity of quadratic edges + // ------------------------------------ + set& 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 * 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::min() && + modPL > std::numeric_limits::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::const_iterator elemIt, elemEnd; + vector< pair > swapVec; + + for ( int dim = 1; dim <= 3; ++dim ) + { + CellsByDimIterator cellsIt( *this, dim ); + while ( const set * 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 faces; + map fgm; + map > linkFacesMap; + map >::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 & fList = lfIt->second; + // list::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 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 + 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 & fList = lfIt->second; + list::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 & 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 & fList = lfIt->second; + list::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::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::const_iterator elemIt, elemEnd; + vector< pair > swapVec; + + CellsByDimIterator cellsIt( *this, 3 ); + while ( const set * 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::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::const_iterator elemIt, elemEnd; + + // numbering _cells of type NORM_POINT1 by node number + { + const set& 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 * typeCells = cellsIt.nextType() ) + { + TID minNumber = std::numeric_limits::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 * 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 * 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::const_iterator elemIt, elemEnd; + for ( int dim = 3; dim > 0; --dim ) + { + CellsByDimIterator dimCells( *this, dim ); + + int nbOfCells = 0; + while ( const std::set * 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 * 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 medGroups; + vector > 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 * 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(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" <::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 ======================================================== + diff --git a/src/MEDLoader/SauvMedConvertor.hxx b/src/MEDLoader/SauvMedConvertor.hxx new file mode 100644 index 000000000..a614778a9 --- /dev/null +++ b/src/MEDLoader/SauvMedConvertor.hxx @@ -0,0 +1,383 @@ +// 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 +#include +#include +#include +#include + +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 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 _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 _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 _comp_names; // component names + std::vector _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 * 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 _coords; + std::vector _groups; + std::vector _nodeFields; + std::vector _cellFields; + + // IMP 0020434: mapping GIBI names to MED names + std::list _listGIBItoMED_mail; // to read from table "MED_MAIL" of PILE_TABLES + std::list _listGIBItoMED_cham; // to read from table "MED_CHAM" of PILE_TABLES + std::list _listGIBItoMED_comp; // to read from table "MED_COMP" of PILE_TABLES + std::map _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 diff --git a/src/MEDLoader/SauvReader.cxx b/src/MEDLoader/SauvReader.cxx new file mode 100644 index 000000000..907f6ba34 --- /dev/null +++ b/src/MEDLoader/SauvReader.cxx @@ -0,0 +1,1149 @@ +// 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 +#include +#include + +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( _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_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 objectNames(nbNamedObjects); + for ( initNameReading( nbNamedObjects ); more(); next() ) + objectNames[ index() ] = getName(); + + vector 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& objectNames, + std::vector& nameIndices) +{ + _iMed->_groups.reserve(nbObjects*2); // fields may add some groups + + char* line; + map 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::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::vector&) +{ + 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&, std::vector&) +{ + if ( isXRD() ) + { + initIntReading(1); + int nb_vals = getIntNext(); + initIntReading(nb_vals); + for(int i=0; i&, std::vector&) +{ + if ( isXRD() ) + { + initIntReading(1); + int nb_vals = getIntNext(); + initDoubleReading(nb_vals); + for(int i=0; i&, std::vector&) +{ + if ( isXRD() ) + { + initIntReading(1); + int nb_vals = getIntNext(); + initIntReading(nb_vals); + for(int i=0; i&, std::vector&) +{ + 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::vector&) +{ + 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::vector&) +{ + 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::vector&) +{ + 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& supports) +{ + SauvUtilities::Group* group = NULL; + set 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 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& fields, + const vector& objets_nommes, + const vector& 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& objectNames, + std::vector& 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 supports( nb_sub ); + vector nb_values ( nb_sub ); + vector 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& 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& objectNames, + std::vector& 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) 2 6 + 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; + } +} diff --git a/src/MEDLoader/SauvReader.hxx b/src/MEDLoader/SauvReader.hxx new file mode 100644 index 000000000..359c135b3 --- /dev/null +++ b/src/MEDLoader/SauvReader.hxx @@ -0,0 +1,102 @@ +// 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 diff --git a/src/MEDLoader/SauvUtilities.hxx b/src/MEDLoader/SauvUtilities.hxx new file mode 100644 index 000000000..a9f476c66 --- /dev/null +++ b/src/MEDLoader/SauvUtilities.hxx @@ -0,0 +1,117 @@ +// 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 diff --git a/src/MEDLoader/SauvWriter.cxx b/src/MEDLoader/SauvWriter.cxx new file mode 100644 index 000000000..c842270cb --- /dev/null +++ b/src/MEDLoader/SauvWriter.cxx @@ -0,0 +1,1323 @@ +// 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; +} diff --git a/src/MEDLoader/SauvWriter.hxx b/src/MEDLoader/SauvWriter.hxx new file mode 100644 index 000000000..0bea6ac68 --- /dev/null +++ b/src/MEDLoader/SauvWriter.hxx @@ -0,0 +1,114 @@ +// 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 diff --git a/src/MEDLoader/Swig/MEDLoader.i b/src/MEDLoader/Swig/MEDLoader.i index aa16d7c95..3bdfbd25f 100644 --- a/src/MEDLoader/Swig/MEDLoader.i +++ b/src/MEDLoader/Swig/MEDLoader.i @@ -30,6 +30,8 @@ #include "MEDFileField.hxx" #include "MEDFileData.hxx" #include "MEDLoaderTypemaps.i" +#include "SauvReader.hxx" +#include "SauvWriter.hxx" using namespace ParaMEDMEM; %} @@ -97,6 +99,10 @@ 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();" @@ -107,6 +113,8 @@ using namespace ParaMEDMEM; %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 { @@ -694,6 +702,7 @@ namespace ParaMEDMEM 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 @@ -856,4 +865,20 @@ namespace ParaMEDMEM } } }; + + 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); + }; + } diff --git a/src/MEDLoader/Swig/SauvLoaderTest.py b/src/MEDLoader/Swig/SauvLoaderTest.py new file mode 100644 index 000000000..1b0a07a7e --- /dev/null +++ b/src/MEDLoader/Swig/SauvLoaderTest.py @@ -0,0 +1,162 @@ +# -*- 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() diff --git a/src/MEDLoader/Test/Makefile.am b/src/MEDLoader/Test/Makefile.am index 3e57ba030..0d008b8f5 100755 --- a/src/MEDLoader/Test/Makefile.am +++ b/src/MEDLoader/Test/Makefile.am @@ -19,21 +19,23 @@ 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 diff --git a/src/MEDLoader/Test/SauvLoaderTest.cxx b/src/MEDLoader/Test/SauvLoaderTest.cxx new file mode 100644 index 000000000..5b1d848bf --- /dev/null +++ b/src/MEDLoader/Test/SauvLoaderTest.cxx @@ -0,0 +1,237 @@ +// 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; +} diff --git a/src/MEDLoader/Test/SauvLoaderTest.hxx b/src/MEDLoader/Test/SauvLoaderTest.hxx new file mode 100644 index 000000000..eea9b4638 --- /dev/null +++ b/src/MEDLoader/Test/SauvLoaderTest.hxx @@ -0,0 +1,45 @@ +// 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 diff --git a/src/MEDLoader/Test/TestSauvLoader.cxx b/src/MEDLoader/Test/TestSauvLoader.cxx new file mode 100644 index 000000000..b9fe7f34f --- /dev/null +++ b/src/MEDLoader/Test/TestSauvLoader.cxx @@ -0,0 +1,25 @@ +// 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" -- 2.39.2