Back to index

salome-med  6.5.0
SauvMedConvertor.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D
00002 //
00003 // This library is free software; you can redistribute it and/or
00004 // modify it under the terms of the GNU Lesser General Public
00005 // License as published by the Free Software Foundation; either
00006 // version 2.1 of the License.
00007 //
00008 // This library is distributed in the hope that it will be useful,
00009 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00010 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00011 // Lesser General Public License for more details.
00012 //
00013 // You should have received a copy of the GNU Lesser General Public
00014 // License along with this library; if not, write to the Free Software
00015 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00016 //
00017 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00018 //
00019 // File      : SauvMedConvertor.cxx
00020 // Created   : Tue Aug 16 14:43:20 2011
00021 // Author    : Edward AGAPOV (eap)
00022 //
00023 
00024 #include "SauvMedConvertor.hxx"
00025 
00026 #include "CellModel.hxx"
00027 #include "MEDFileMesh.hxx"
00028 #include "MEDFileField.hxx"
00029 #include "MEDFileData.hxx"
00030 #include "MEDCouplingFieldDouble.hxx"
00031 
00032 #include <iostream>
00033 #include <cassert>
00034 #include <cmath>
00035 #include <queue>
00036 #include <limits>
00037 
00038 #include <cstdlib>
00039 #include <cstring>
00040 #include <fcntl.h>
00041 
00042 #ifdef WNT
00043 #include <io.h>
00044 #endif
00045 
00046 #ifndef WNT
00047 #define HAS_XDR
00048 #endif
00049 
00050 #ifdef HAS_XDR
00051 #include <rpc/xdr.h>
00052 #endif
00053 
00054 using namespace SauvUtilities;
00055 using namespace ParaMEDMEM;
00056 using namespace std;
00057 
00058 namespace
00059 {
00060   // for ASCII file reader
00061   const int GIBI_MaxOutputLen = 150; // max length of a line in the sauve file
00062   const int GIBI_BufferSize   = 16184; // for buffered reading
00063 
00064   using namespace INTERP_KERNEL;
00065 
00066   const size_t MaxMedCellType = NORM_ERROR;
00067   const size_t NbGibiCellTypes = 47;
00068   const TCellType GibiTypeToMed[NbGibiCellTypes] =
00069     {
00070       /*1 */ NORM_POINT1 ,/*2 */ NORM_SEG2   ,/*3 */ NORM_SEG3   ,/*4 */ NORM_TRI3   ,/*5 */ NORM_ERROR  ,
00071       /*6 */ NORM_TRI6   ,/*7 */ NORM_ERROR  ,/*8 */ NORM_QUAD4  ,/*9 */ NORM_ERROR  ,/*10*/ NORM_QUAD8  ,
00072       /*11*/ NORM_ERROR  ,/*12*/ NORM_ERROR  ,/*13*/ NORM_ERROR  ,/*14*/ NORM_HEXA8  ,/*15*/ NORM_HEXA20 ,
00073       /*16*/ NORM_PENTA6 ,/*17*/ NORM_PENTA15,/*18*/ NORM_ERROR  ,/*19*/ NORM_ERROR  ,/*20*/ NORM_ERROR  ,
00074       /*21*/ NORM_ERROR  ,/*22*/ NORM_ERROR  ,/*23*/ NORM_TETRA4 ,/*24*/ NORM_TETRA10,/*25*/ NORM_PYRA5  ,
00075       /*26*/ NORM_PYRA13 ,/*27*/ NORM_ERROR  ,/*28*/ NORM_ERROR  ,/*29*/ NORM_ERROR  ,/*30*/ NORM_ERROR  ,
00076       /*31*/ NORM_ERROR  ,/*32*/ NORM_ERROR  ,/*33*/ NORM_ERROR  ,/*34*/ NORM_ERROR  ,/*35*/ NORM_ERROR  ,
00077       /*36*/ NORM_ERROR  ,/*37*/ NORM_ERROR  ,/*38*/ NORM_ERROR  ,/*39*/ NORM_ERROR  ,/*40*/ NORM_ERROR  ,
00078       /*41*/ NORM_ERROR  ,/*42*/ NORM_ERROR  ,/*43*/ NORM_ERROR  ,/*44*/ NORM_ERROR  ,/*45*/ NORM_ERROR  ,
00079       /*46*/ NORM_ERROR  ,/*47*/ NORM_ERROR
00080     };
00081 
00082   //================================================================================
00086   //================================================================================
00087 
00088   unsigned getDim( const Group* grp )
00089   {
00090     return SauvUtilities::getDimension( grp->_groups.empty() ? grp->_cellType : grp->_groups[0]->_cellType );
00091   }
00092 
00093   //================================================================================
00097   //================================================================================
00098 
00099   inline void ConvertQuadratic( const INTERP_KERNEL::NormalizedCellType type,
00100                                 const Cell &                     aCell )
00101   {
00102     if ( const int * conn = getGibi2MedQuadraticInterlace( type ))
00103       {
00104         Cell* ma = const_cast<Cell*>(&aCell);
00105         //cout << "###### BEFORE ConvertQuadratic() " << *ma << endl;
00106         vector< Node* > new_nodes( ma->_nodes.size() );
00107         for ( size_t i = 0; i < new_nodes.size(); ++i )
00108           new_nodes[ i ] = ma->_nodes[ conn[ i ]];
00109         ma->_nodes.swap( new_nodes );
00110         //cout << "###### AFTER ConvertQuadratic() " << *ma << endl;
00111       }
00112   }
00113 
00114   //================================================================================
00118   //================================================================================
00119 
00120   void getReverseVector (const INTERP_KERNEL::NormalizedCellType type,
00121                          vector<pair<int,int> > &                swapVec )
00122   {
00123     swapVec.clear();
00124 
00125     switch ( type )
00126       {
00127       case NORM_TETRA4:
00128         swapVec.resize(1);
00129         swapVec[0] = make_pair( 1, 2 );
00130         break;
00131       case NORM_PYRA5:
00132         swapVec.resize(1);
00133         swapVec[0] = make_pair( 1, 3 );
00134         break;
00135       case NORM_PENTA6:
00136         swapVec.resize(2);
00137         swapVec[0] = make_pair( 1, 2 );
00138         swapVec[1] = make_pair( 4, 5 );
00139         break;
00140       case NORM_HEXA8:
00141         swapVec.resize(2);
00142         swapVec[0] = make_pair( 1, 3 );
00143         swapVec[1] = make_pair( 5, 7 );
00144         break;
00145       case NORM_TETRA10:
00146         swapVec.resize(3);
00147         swapVec[0] = make_pair( 1, 2 );
00148         swapVec[1] = make_pair( 4, 6 );
00149         swapVec[2] = make_pair( 8, 9 );
00150         break;
00151       case NORM_PYRA13:
00152         swapVec.resize(4);
00153         swapVec[0] = make_pair( 1, 3 );
00154         swapVec[1] = make_pair( 5, 8 );
00155         swapVec[2] = make_pair( 6, 7 );
00156         swapVec[3] = make_pair( 10, 12 );
00157         break;
00158       case NORM_PENTA15:
00159         swapVec.resize(4);
00160         swapVec[0] = make_pair( 1, 2 );
00161         swapVec[1] = make_pair( 4, 5 );
00162         swapVec[2] = make_pair( 6, 8 );
00163         swapVec[3] = make_pair( 9, 11 );
00164         break;
00165       case NORM_HEXA20:
00166         swapVec.resize(7);
00167         swapVec[0] = make_pair( 1, 3 );
00168         swapVec[1] = make_pair( 5, 7 );
00169         swapVec[2] = make_pair( 8, 11 );
00170         swapVec[3] = make_pair( 9, 10 );
00171         swapVec[4] = make_pair( 12, 15 );
00172         swapVec[5] = make_pair( 13, 14 );
00173         swapVec[6] = make_pair( 17, 19 );
00174         break;
00175         //   case NORM_SEG3: no need to reverse edges
00176         //     swapVec.resize(1);
00177         //     swapVec[0] = make_pair( 1, 2 );
00178         //     break;
00179       case NORM_TRI6:
00180         swapVec.resize(2);
00181         swapVec[0] = make_pair( 1, 2 );
00182         swapVec[1] = make_pair( 3, 5 );
00183         break;
00184       case NORM_QUAD8:
00185         swapVec.resize(3);
00186         swapVec[0] = make_pair( 1, 3 );
00187         swapVec[1] = make_pair( 4, 7 );
00188         swapVec[2] = make_pair( 5, 6 );
00189         break;
00190       default:;
00191       }
00192   }
00193 
00194   //================================================================================
00198   //================================================================================
00199 
00200   inline void reverse(const Cell & aCell, const vector<pair<int,int> > & swapVec )
00201   {
00202     Cell* ma = const_cast<Cell*>(&aCell);
00203     for ( unsigned i = 0; i < swapVec.size(); ++i )
00204       std::swap( ma->_nodes[ swapVec[i].first ],
00205                  ma->_nodes[ swapVec[i].second ]);
00206     if ( swapVec.empty() )
00207       ma->_reverse = true;
00208     else
00209       ma->_reverse = false;
00210   }
00211   //================================================================================
00215   struct TCellByIDCompare
00216   {
00217     bool operator () (const Cell* i1, const Cell* i2)
00218     {
00219       return i1->_number < i2->_number;
00220     }
00221   };
00222   typedef map< const Cell*, unsigned, TCellByIDCompare > TCellToOrderMap;
00223 
00224   //================================================================================
00228   //================================================================================
00229 
00230   void setRelocationTable( Group* grp, TCellToOrderMap& cell2order )
00231   {
00232     if ( !grp->_isProfile ) return;
00233 
00234     // check if relocation table is necessary
00235     bool isRelocated = false;
00236     unsigned newOrder = 0;
00237     TCellToOrderMap::iterator c2oIt = cell2order.begin(), c2oEnd = cell2order.end();
00238     for ( ; !isRelocated && c2oIt != c2oEnd; ++c2oIt, ++newOrder )
00239       isRelocated = ( c2oIt->second != newOrder );
00240 
00241     if ( isRelocated )
00242     {
00243       grp->_relocTable.resize( cell2order.size() );
00244       for ( newOrder = 0, c2oIt = cell2order.begin(); c2oIt != c2oEnd; ++c2oIt, ++newOrder )
00245         grp->_relocTable[ c2oIt->second ] = newOrder;
00246     }
00247   }
00248 }
00249 
00250 //================================================================================
00254 //================================================================================
00255 
00256 unsigned SauvUtilities::getDimension( INTERP_KERNEL::NormalizedCellType type )
00257 {
00258   return type == NORM_ERROR ? -1 : INTERP_KERNEL::CellModel::GetCellModel( type ).getDimension();
00259 }
00260 
00261 //================================================================================
00265 //================================================================================
00266 
00267 const int * SauvUtilities::getGibi2MedQuadraticInterlace( INTERP_KERNEL::NormalizedCellType type )
00268 {
00269   static vector<const int*> conn;
00270   static const int hexa20 [] = {0,6,4,2, 12,18,16,14, 7,5,3,1, 19,17,15,13, 8,11,10,9};
00271   static const int penta15[] = {0,2,4, 9,11,13, 1,3,5, 10,12,14, 6,7,3};
00272   static const int pyra13 [] = {0,2,4,6, 12, 1,3,5,7, 8,9,10,11};
00273   static const int tetra10[] = {0,2,4, 9, 1,3,5, 6,7,8};
00274   static const int quad8  [] = {0,2,4,6, 1,3,5,7};
00275   static const int tria6  [] = {0,2,4, 1,3,5};
00276   static const int seg3   [] = {0,2,1};
00277   if ( conn.empty() )
00278   {
00279     conn.resize( MaxMedCellType + 1, 0 );
00280     conn[ NORM_HEXA20 ] = hexa20;
00281     conn[ NORM_PENTA15] = penta15;
00282     conn[ NORM_PYRA13 ] = pyra13;
00283     conn[ NORM_TETRA10] = tetra10;
00284     conn[ NORM_SEG3   ] = seg3;
00285     conn[ NORM_TRI6   ] = tria6;
00286     conn[ NORM_QUAD8  ] = quad8;
00287   }
00288   return conn[ type ];
00289 }
00290 
00291 //================================================================================
00295 //================================================================================
00296 
00297 Cell::Cell(const Cell& ma)
00298   : _nodes(ma._nodes), _reverse(ma._reverse), _sortedNodeIDs(0), _number(ma._number)
00299 {
00300   if ( ma._sortedNodeIDs )
00301     {
00302       _sortedNodeIDs = new int[ _nodes.size() ];
00303       std::copy( ma._sortedNodeIDs, ma._sortedNodeIDs + _nodes.size(), _sortedNodeIDs );
00304     }
00305 }
00306 
00307 //================================================================================
00311 //================================================================================
00312 
00313 SauvUtilities::Link Cell::link(int i) const
00314 {
00315   int i2 = ( i + 1 ) % _nodes.size();
00316   if ( _reverse )
00317     return make_pair( _nodes[i2]->_number, _nodes[i]->_number );
00318   else
00319     return make_pair( _nodes[i]->_number, _nodes[i2]->_number );
00320 }
00321 
00322 //================================================================================
00326 //================================================================================
00327 
00328 const TID* Cell::getSortedNodes() const
00329 {
00330   if ( !_sortedNodeIDs )
00331   {
00332     size_t l=_nodes.size();
00333     _sortedNodeIDs = new int[ l ];
00334 
00335     for (size_t i=0; i!=l; ++i)
00336       _sortedNodeIDs[i]=_nodes[i]->_number;
00337     std::sort( _sortedNodeIDs, _sortedNodeIDs + l );
00338   }
00339   return _sortedNodeIDs;
00340 }
00341 
00342 //================================================================================
00346 //================================================================================
00347 
00348 bool Cell::operator< (const Cell& ma) const
00349 {
00350   if ( _nodes.size() == 1 )
00351     return _nodes[0] < ma._nodes[0];
00352 
00353   const int* v1 = getSortedNodes();
00354   const int* v2 = ma.getSortedNodes();
00355   for ( const int* vEnd = v1 + _nodes.size(); v1 < vEnd; ++v1, ++v2 )
00356     if(*v1 != *v2)
00357       return *v1 < *v2;
00358   return false;
00359 }
00360 
00361 //================================================================================
00365 //================================================================================
00366 
00367 std::ostream& SauvUtilities::operator<< (std::ostream& os, const SauvUtilities::Cell& ma)
00368 {
00369   os << "cell " << ma._number << " (" << ma._nodes.size() << " nodes) : < " << ma._nodes[0]->_number;
00370   for( size_t i=1; i!=ma._nodes.size(); ++i)
00371     os << ", " << ma._nodes[i]->_number;
00372 #ifdef _DEBUG_
00373   os << " > sortedNodes: ";
00374   if ( ma._sortedNodeIDs ) {
00375     os << "< ";
00376     for( size_t i=0; i!=ma._nodes.size(); ++i)
00377       os << ( i ? ", " : "" ) << ma._sortedNodeIDs[i];
00378     os << " >";
00379   }
00380   else {
00381     os << "NULL";
00382   }
00383 #endif
00384   return os;
00385 }
00386 
00387 //================================================================================
00391 //================================================================================
00392 
00393 int Group::size() const
00394 {
00395   int sizze = 0;
00396   if ( !_relocTable.empty() )
00397     sizze =  _relocTable.size();
00398   else if ( _medGroup )
00399     sizze = _medGroup->getNumberOfTuples();
00400   else if ( !_cells.empty() )
00401     sizze = _cells.size();
00402   else
00403     for ( size_t i = 0; i < _groups.size(); ++i )
00404       sizze += _groups[i]->size();
00405   return sizze;
00406 }
00407 
00408 //================================================================================
00412 //================================================================================
00413 
00414 INTERP_KERNEL::NormalizedCellType SauvUtilities::gibi2medGeom( size_t gibiType )
00415 {
00416   if ( gibiType < 1 || gibiType > NbGibiCellTypes )
00417     return NORM_ERROR;
00418 
00419   return GibiTypeToMed[ gibiType - 1 ];
00420 }
00421 
00422 //================================================================================
00426 //================================================================================
00427 
00428 int SauvUtilities::med2gibiGeom( INTERP_KERNEL::NormalizedCellType medGeomType )
00429 {
00430   for ( unsigned int i = 0; i < NbGibiCellTypes; i++ )
00431     if ( GibiTypeToMed[ i ] == medGeomType )
00432       return i + 1;
00433 
00434   return -1;
00435 }
00436 
00437 //================================================================================
00441 //================================================================================
00442 
00443 FileReader::FileReader(const char* fileName):_fileName(fileName),_iRead(0),_nbToRead(0)
00444 {
00445 }
00446 
00447 //================================================================================
00451 //================================================================================
00452 
00453 ASCIIReader::ASCIIReader(const char* fileName)
00454   :FileReader(fileName),
00455    _file(-1)
00456 {
00457 }
00458 
00459 //================================================================================
00463 //================================================================================
00464 
00465 bool ASCIIReader::isASCII() const
00466 {
00467   return true;
00468 }
00469 
00470 //================================================================================
00474 //================================================================================
00475 
00476 bool ASCIIReader::open()
00477 {
00478 #ifdef WNT
00479   _file = ::_open (_fileName.c_str(), _O_RDONLY|_O_BINARY);
00480 #else
00481   _file = ::open (_fileName.c_str(), O_RDONLY);
00482 #endif
00483   if (_file >= 0)
00484     {
00485       _start  = new char [GIBI_BufferSize]; // working buffer beginning
00486       //_tmpBuf = new char [GIBI_MaxOutputLen];
00487       _ptr    = _start;
00488       _eptr   = _start;
00489       _lineNb = 0;
00490     }
00491   else
00492     {
00493       //THROW_IK_EXCEPTION("Can't open file "<<_fileName << " fd: " << _file);
00494     }
00495   return (_file >= 0);
00496 }
00497 
00498 //================================================================================
00502 //================================================================================
00503 
00504 ASCIIReader::~ASCIIReader()
00505 {
00506   if (_file >= 0)
00507     {
00508       ::close (_file);
00509       if (_start != 0L)
00510         {
00511           delete [] _start;
00512           //delete [] _tmpBuf;
00513           _start = 0;
00514         }
00515       _file = -1;
00516     }
00517 }
00518 
00519 //================================================================================
00523 //================================================================================
00524 
00525 bool ASCIIReader::getNextLine (char* & line, bool raiseOEF /*= true*/ )
00526 {
00527   if ( getLine( line )) return true;
00528   if ( raiseOEF )
00529     THROW_IK_EXCEPTION("Unexpected EOF on ln "<<_lineNb);
00530   return false;
00531 }
00532 
00533 //================================================================================
00537 //================================================================================
00538 
00539 bool ASCIIReader::getLine(char* & line)
00540 {
00541   bool aResult = true;
00542   // Check the state of the buffer;
00543   // if there is too little left, read the next portion of data
00544   int nBytesRest = _eptr - _ptr;
00545   if (nBytesRest < GIBI_MaxOutputLen)
00546     {
00547       if (nBytesRest > 0)
00548         {
00549           // move the remaining portion to the buffer beginning
00550           for ( int i = 0; i < nBytesRest; ++i )
00551             _start[i] = _ptr[i];
00552           //memcpy (_tmpBuf, _ptr, nBytesRest);
00553           //memcpy (_start, _tmpBuf, nBytesRest);
00554         }
00555       else
00556         {
00557           nBytesRest = 0;
00558         }
00559       _ptr = _start;
00560       const int nBytesRead = ::read (_file,
00561                                      &_start [nBytesRest],
00562                                      GIBI_BufferSize - nBytesRest);
00563       nBytesRest += nBytesRead;
00564       _eptr = &_start [nBytesRest];
00565     }
00566   // Check the buffer for the end-of-line
00567   char * ptr = _ptr;
00568   while (true)
00569     {
00570       // Check for end-of-the-buffer, the ultimate criterion for termination
00571       if (ptr >= _eptr)
00572         {
00573           if (nBytesRest <= 0)
00574             aResult = false;
00575           else
00576             _eptr[-1] = '\0';
00577           break;
00578         }
00579       // seek the line-feed character
00580       if (ptr[0] == '\n')
00581         {
00582           if (ptr[-1] == '\r')
00583             ptr[-1] = '\0';
00584           ptr[0] = '\0';
00585           ++ptr;
00586           break;
00587         }
00588       ++ptr;
00589     }
00590   // Output the result
00591   line = _ptr;
00592   _ptr = ptr;
00593   _lineNb++;
00594 
00595   return aResult;
00596 }
00597 
00598 //================================================================================
00606 //================================================================================
00607 
00608 void ASCIIReader::init( int nbToRead, int nbPosInLine, int width, int shift /*= 0*/ )
00609 {
00610   _nbToRead    = nbToRead;
00611   _nbPosInLine = nbPosInLine;
00612   _width       = width;
00613   _shift       = shift;
00614   _iPos = _iRead = 0;
00615   if ( _nbToRead )
00616     {
00617       getNextLine( _curPos );
00618       _curPos = _curPos + _shift;
00619     }
00620   else
00621     {
00622       _curPos = 0;
00623     }
00624   _curLocale.clear();
00625 }
00626 
00627 //================================================================================
00631 //================================================================================
00632 
00633 void ASCIIReader::initNameReading(int nbValues, int width /*= 8*/)
00634 {
00635   init( nbValues, 72 / ( width + 1 ), width, 1 );
00636 }
00637 
00638 //================================================================================
00642 //================================================================================
00643 
00644 void ASCIIReader::initIntReading(int nbValues)
00645 {
00646   init( nbValues, 10, 8 );
00647 }
00648 
00649 //================================================================================
00653 //================================================================================
00654 
00655 void ASCIIReader::initDoubleReading(int nbValues)
00656 {
00657   init( nbValues, 3, 22 );
00658 
00659   // Correction 2 of getDouble(): set "C" numeric locale to read numbers
00660   // with dot decimal point separator, as it is in SAUVE files
00661   _curLocale = setlocale(LC_NUMERIC, "C");
00662 }
00663 
00664 //================================================================================
00668 //================================================================================
00669 
00670 bool ASCIIReader::more() const
00671 {
00672   bool result = false;
00673   if ( _iRead < _nbToRead)
00674     {
00675       if ( _curPos ) result = true;
00676     }
00677   return result;
00678 }
00679 
00680 //================================================================================
00684 //================================================================================
00685 
00686 void ASCIIReader::next()
00687 {
00688   if ( !more() )
00689     THROW_IK_EXCEPTION("SauvUtilities::ASCIIReader::next(): no more() values to read");
00690   ++_iRead;
00691   ++_iPos;
00692   if ( _iRead < _nbToRead )
00693     {
00694       if ( _iPos >= _nbPosInLine )
00695         {
00696           getNextLine( _curPos );
00697           _curPos = _curPos + _shift;
00698           _iPos = 0;
00699         }
00700       else
00701         {
00702           _curPos = _curPos + _width + _shift;
00703         }
00704     }
00705   else
00706     {
00707       _curPos = 0;
00708       if ( !_curLocale.empty() )
00709         {
00710           setlocale(LC_NUMERIC, _curLocale.c_str());
00711           _curLocale.clear();
00712         }
00713     }
00714 }
00715 
00716 //================================================================================
00720 //================================================================================
00721 
00722 int ASCIIReader::getInt() const
00723 {
00724   // fix for two glued ints (issue 0021009):
00725   // Line nb    |   File contents
00726   // ------------------------------------------------------------------------------------
00727   // 53619905   |       1       2       6       8
00728   // 53619906   |                                                                SCALAIRE
00729   // 53619907   |    -63312600499       1       0       0       0      -2       0       2
00730   //   where -63312600499 is actualy -633 and 12600499
00731   char hold=_curPos[_width];
00732   _curPos[_width] = '\0';
00733   int result = atoi( _curPos );
00734   _curPos[_width] = hold;
00735   return result;
00736   //return atoi(str());
00737 }
00738 
00739 //================================================================================
00743 //================================================================================
00744 
00745 float ASCIIReader::getFloat() const
00746 {
00747   return getDouble();
00748 }
00749 
00750 //================================================================================
00754 //================================================================================
00755 
00756 double ASCIIReader::getDouble() const
00757 {
00758   //std::string aStr (_curPos);
00759 
00760   // Correction: add missing 'E' specifier
00761   // int aPosStart = aStr.find_first_not_of(" \t");
00762   // if (aPosStart < (int)aStr.length()) {
00763   //   int aPosSign = aStr.find_first_of("+-", aPosStart + 1); // pass the leading symbol, as it can be a sign
00764   //   if (aPosSign < (int)aStr.length()) {
00765   //     if (aStr[aPosSign - 1] != 'e' && aStr[aPosSign - 1] != 'E')
00766   //       aStr.insert(aPosSign, "E", 1);
00767   //   }
00768   // }
00769 
00770   // Different Correction (more optimal)
00771   // Sample:
00772   //  0.00000000000000E+00 -2.37822406690632E+01  6.03062748797469E+01
00773   //  7.70000000000000-100  7.70000000000000+100  7.70000000000000+100
00774   //0123456789012345678901234567890123456789012345678901234567890123456789
00775   const size_t posE = 18;
00776   if ( _curPos[posE] != 'E' && _curPos[posE] != 'e' )
00777     {
00778       std::string aStr (_curPos);
00779       if ( aStr.size() < posE+1 )
00780         THROW_IK_EXCEPTION("No more doubles (line #" << lineNb() << ")");
00781       aStr.insert( posE, "E", 1 );
00782       return atof(aStr.c_str());
00783     }
00784   return atof( _curPos );
00785 }
00786 
00787 //================================================================================
00791 //================================================================================
00792 
00793 string ASCIIReader::getName() const
00794 {
00795   int len = _width;
00796   while (( _curPos[len-1] == ' ' || _curPos[len-1] == 0) && len > 0 )
00797     len--;
00798   return string( _curPos, len );
00799 }
00800 
00801 //================================================================================
00805 //================================================================================
00806 
00807 XDRReader::XDRReader(const char* fileName) :FileReader(fileName), _xdrs_file(NULL)
00808 {
00809 }
00810 
00811 //================================================================================
00815 //================================================================================
00816 
00817 XDRReader::~XDRReader()
00818 {
00819 #ifdef HAS_XDR  
00820   if ( _xdrs_file )
00821     {
00822       xdr_destroy((XDR*)_xdrs);
00823       free((XDR*)_xdrs);
00824       ::fclose(_xdrs_file);
00825       _xdrs_file = NULL;
00826     }
00827 #endif
00828 }
00829 
00830 //================================================================================
00834 //================================================================================
00835 
00836 bool XDRReader::isASCII() const
00837 {
00838   return false;
00839 }
00840 
00841 //================================================================================
00845 //================================================================================
00846 
00847 bool XDRReader::open()
00848 {
00849   bool xdr_ok = false;
00850 #ifdef HAS_XDR
00851   if ((_xdrs_file = ::fopen(_fileName.c_str(), "r")))
00852     {
00853       _xdrs = (XDR *)malloc(sizeof(XDR));
00854       xdrstdio_create((XDR*)_xdrs, _xdrs_file, XDR_DECODE);
00855 
00856       const int maxsize = 10;
00857       char icha[maxsize+1];
00858       char* icha2 = icha;
00859       if (( xdr_ok = xdr_string((XDR*)_xdrs, &icha2, maxsize)))
00860         {
00861           icha[maxsize] = '\0';
00862           xdr_ok = (strcmp(icha, "CASTEM XDR") == 0);
00863         }
00864       if ( !xdr_ok )
00865         {
00866           xdr_destroy((XDR*)_xdrs);
00867           free((XDR*)_xdrs);
00868           fclose(_xdrs_file);
00869           _xdrs_file = NULL;
00870         }
00871     }
00872 #endif
00873   return xdr_ok;
00874 }
00875 
00876 //================================================================================
00880 //================================================================================
00881 
00882 bool XDRReader::getNextLine (char* &, bool )
00883 {
00884   return true;
00885 }
00886 
00887 //================================================================================
00893 //================================================================================
00894 
00895 void XDRReader::init( int nbToRead, int width/*=0*/ )
00896 {
00897   if(_iRead < _nbToRead)
00898     {
00899       cout << "_iRead, _nbToRead : " << _iRead << " " << _nbToRead << endl;
00900       cout << "Unfinished iteration before new one !" << endl;
00901       THROW_IK_EXCEPTION("SauvUtilities::XDRReader::init(): Unfinished iteration before new one !");
00902     }
00903   _iRead    = 0;
00904   _nbToRead = nbToRead;
00905   _width    = width;
00906 }
00907 
00908 //================================================================================
00912 //================================================================================
00913 
00914 void XDRReader::initNameReading(int nbValues, int width)
00915 {
00916   init( nbValues, width );
00917   _xdr_kind = _xdr_kind_char;
00918   if(nbValues*width)
00919     {
00920       unsigned int nels = nbValues*width;
00921       _xdr_cvals = (char*)malloc((nels+1)*sizeof(char));
00922 #ifdef HAS_XDR
00923       xdr_string((XDR*)_xdrs, &_xdr_cvals, nels);
00924 #endif
00925       _xdr_cvals[nels] = '\0';
00926     }
00927 }
00928 
00929 //================================================================================
00933 //================================================================================
00934 
00935 void XDRReader::initIntReading(int nbValues)
00936 {
00937   init( nbValues );
00938   _xdr_kind = _xdr_kind_int;
00939   if(nbValues)
00940     {
00941 #ifdef HAS_XDR
00942       unsigned int nels = nbValues;
00943       unsigned int actual_nels;
00944       _xdr_ivals = (int*)malloc(nels*sizeof(int));
00945       xdr_array((XDR*)_xdrs, (char **)&_xdr_ivals, &actual_nels, nels, sizeof(int), (xdrproc_t)xdr_int);
00946 #endif
00947     }
00948 }
00949 
00950 //================================================================================
00954 //================================================================================
00955 
00956 void XDRReader::initDoubleReading(int nbValues)
00957 {
00958   init( nbValues );
00959   _xdr_kind = _xdr_kind_double;
00960   if(nbValues)
00961     {
00962 #ifdef HAS_XDR
00963       unsigned int nels = nbValues;
00964       unsigned int actual_nels;
00965       _xdr_dvals = (double*)malloc(nels*sizeof(double));
00966       xdr_array((XDR*)_xdrs, (char **)&_xdr_dvals, &actual_nels, nels, sizeof(double), (xdrproc_t)xdr_double);
00967 #endif
00968     }
00969 }
00970 
00971 //================================================================================
00975 //================================================================================
00976 
00977 bool XDRReader::more() const
00978 {
00979   return _iRead < _nbToRead;
00980 }
00981 
00982 //================================================================================
00986 //================================================================================
00987 
00988 void XDRReader::next()
00989 {
00990   if ( !more() )
00991     THROW_IK_EXCEPTION("SauvUtilities::XDRReader::next(): no more() values to read");
00992 
00993   ++_iRead;
00994   if ( _iRead < _nbToRead )
00995     {
00996     }
00997   else
00998     {
00999       if(_xdr_kind == _xdr_kind_char) free(_xdr_cvals);
01000       if(_xdr_kind == _xdr_kind_int) free(_xdr_ivals);
01001       if(_xdr_kind == _xdr_kind_double) free(_xdr_dvals);
01002       _xdr_kind = _xdr_kind_null;
01003     }
01004 }
01005 
01006 //================================================================================
01010 //================================================================================
01011 
01012 int XDRReader::getInt() const
01013 {
01014   if(_iRead < _nbToRead)
01015     {
01016       return _xdr_ivals[_iRead];
01017     }
01018   else
01019     {
01020       int result = 0;
01021 #ifdef HAS_XDR
01022       xdr_int((XDR*)_xdrs, &result);
01023 #endif
01024       return result;
01025     }
01026 }
01027 
01028 //================================================================================
01032 //================================================================================
01033 
01034 float  XDRReader::getFloat() const
01035 {
01036   float result = 0;
01037 #ifdef HAS_XDR
01038   xdr_float((XDR*)_xdrs, &result);
01039 #endif
01040   return result;
01041 }
01042 
01043 //================================================================================
01047 //================================================================================
01048 
01049 double XDRReader::getDouble() const
01050 {
01051   if(_iRead < _nbToRead)
01052     {
01053       return _xdr_dvals[_iRead];
01054     }
01055   else
01056     {
01057       double result = 0;
01058 #ifdef HAS_XDR
01059       xdr_double((XDR*)_xdrs, &result);
01060 #endif
01061       return result;
01062     }
01063 }
01064 
01065 //================================================================================
01069 //================================================================================
01070 
01071 std::string XDRReader::getName() const
01072 {
01073   int len = _width;
01074   char* s = _xdr_cvals + _iRead*_width;
01075   while (( s[len-1] == ' ' || s[len-1] == 0) && len > 0 )
01076     len--;
01077   return string( s, len );
01078 }
01079 
01080 //================================================================================
01084 //================================================================================
01085 
01086 void IntermediateMED::checkDataAvailability() const throw(INTERP_KERNEL::Exception)
01087 {
01088   if ( _spaceDim == 0 )
01089     THROW_IK_EXCEPTION("Wrong file format"); // it is the first record in the sauve file
01090 
01091   if ( _groups.empty() )
01092     THROW_IK_EXCEPTION("No elements have been read");
01093 
01094   if ( _points.empty() || _nbNodes == 0 )
01095     THROW_IK_EXCEPTION("Nodes of elements are not filled");
01096 
01097   if ( _coords.empty() )
01098     THROW_IK_EXCEPTION("Node coordinates are missing");
01099 
01100   if ( _coords.size() < _nbNodes * _spaceDim )
01101     THROW_IK_EXCEPTION("Nodes and coordinates mismatch");
01102 }
01103 
01104 //================================================================================
01108 //================================================================================
01109 
01110 ParaMEDMEM::MEDFileData* IntermediateMED::convertInMEDFileDS()
01111 {
01112   MEDCouplingAutoRefCountObjectPtr< MEDFileUMesh >  mesh   = makeMEDFileMesh();
01113   MEDCouplingAutoRefCountObjectPtr< MEDFileFields > fields = makeMEDFileFields(mesh);
01114 
01115   MEDCouplingAutoRefCountObjectPtr< MEDFileMeshes > meshes = MEDFileMeshes::New();
01116   MEDCouplingAutoRefCountObjectPtr< MEDFileData >  medData = MEDFileData::New();
01117   meshes->pushMesh( mesh );
01118   medData->setMeshes( meshes );
01119   if ( fields ) medData->setFields( fields );
01120 
01121   medData->incrRef();
01122   return medData;
01123 }
01124 
01125 //================================================================================
01129 //================================================================================
01130 
01131 ParaMEDMEM::MEDFileUMesh* IntermediateMED::makeMEDFileMesh()
01132 {
01133   // check if all needed piles are present
01134   checkDataAvailability();
01135 
01136   // set long names
01137   setGroupLongNames();
01138 
01139   // fix element orientation
01140   if ( _spaceDim == 2 )
01141     orientElements2D();
01142   else if ( _spaceDim == 3 )
01143     orientElements3D();
01144 
01145   // process groups
01146   decreaseHierarchicalDepthOfSubgroups();
01147   eraseUselessGroups();
01148   detectMixDimGroups();
01149 
01150   // assign IDs
01151   _points.numberNodes();
01152   numberElements();
01153 
01154   // make the med mesh
01155 
01156   MEDFileUMesh* mesh = MEDFileUMesh::New();
01157 
01158   DataArrayDouble *coords = getCoords();
01159   setConnectivity( mesh, coords );
01160   setGroups( mesh );
01161 
01162   coords->decrRef();
01163 
01164   if ( !mesh->getName() || strlen( mesh->getName() ) == 0 )
01165     mesh->setName( "MESH" );
01166 
01167   return mesh;
01168 }
01169 
01170 //================================================================================
01174 //================================================================================
01175 
01176 void IntermediateMED::setGroupLongNames()
01177 {
01178   // IMP 0020434: mapping GIBI names to MED names
01179   // set med names to objects (mesh, fields, support, group or other)
01180 
01181   set<int> treatedGroups;
01182 
01183   list<nameGIBItoMED>::iterator itGIBItoMED = _listGIBItoMED_mail.begin();
01184   for (; itGIBItoMED != _listGIBItoMED_mail.end(); itGIBItoMED++)
01185     {
01186       if ( (int)_groups.size() < itGIBItoMED->gibi_id ) continue;
01187 
01188       SauvUtilities::Group & grp = _groups[itGIBItoMED->gibi_id - 1];
01189 
01190       // if there are several names for grp then the 1st name is the name
01191       // of grp and the rest ones are names of groups referring grp (issue 0021311)
01192       const bool isRefName = !treatedGroups.insert( itGIBItoMED->gibi_id ).second;
01193       if ( !isRefName )
01194         grp._name = _mapStrings[ itGIBItoMED->med_id ];
01195       else
01196         for ( unsigned i = 0; i < grp._refNames.size(); ++i )
01197           if ( grp._refNames[i].empty() )
01198             grp._refNames[i] = _mapStrings[ (*itGIBItoMED).med_id ];
01199     }
01200 }
01201 
01202 //================================================================================
01206 //================================================================================
01207 
01208 void IntermediateMED::setFieldLongNames(set< string >& usedNames)
01209 {
01210   list<nameGIBItoMED>::iterator itGIBItoMED = _listGIBItoMED_cham.begin();
01211   for (; itGIBItoMED != _listGIBItoMED_cham.end(); itGIBItoMED++)
01212     {
01213       if (itGIBItoMED->gibi_pile == PILE_FIELD)
01214         {
01215           _cellFields[itGIBItoMED->gibi_id - 1]->_name = _mapStrings[itGIBItoMED->med_id];
01216         }
01217       else if (itGIBItoMED->gibi_pile == PILE_NODES_FIELD)
01218         {
01219           _nodeFields[itGIBItoMED->gibi_id - 1]->_name = _mapStrings[itGIBItoMED->med_id];
01220         }
01221     } // iterate on _listGIBItoMED_cham
01222 
01223   for (itGIBItoMED =_listGIBItoMED_comp.begin(); itGIBItoMED != _listGIBItoMED_comp.end(); itGIBItoMED++)
01224     {
01225       string medName  = _mapStrings[itGIBItoMED->med_id];
01226       string gibiName = _mapStrings[itGIBItoMED->gibi_id];
01227 
01228       bool name_found = false;
01229       for ( int isNodal = 0; isNodal < 2 && !name_found; ++isNodal )
01230         {
01231           vector<DoubleField* > & fields = isNodal ? _nodeFields : _cellFields;
01232           for ( size_t ifi = 0; ifi < fields.size() && !name_found; ifi++)
01233             {
01234               if (medName.find( fields[ifi]->_name + "." ) == 0 )
01235                 {
01236                   vector<DoubleField::_Sub_data>& aSubDs = fields[ifi]->_sub;
01237                   int nbSub = aSubDs.size();
01238                   for (int isu = 0; isu < nbSub; isu++)
01239                     for (int ico = 0; ico < aSubDs[isu].nbComponents(); ico++)
01240                       {
01241                         if (aSubDs[isu].compName(ico) == gibiName)
01242                           {
01243                             string medNameCompo = medName.substr( fields[ifi]->_name.size() + 1 );
01244                             fields[ifi]->_sub[isu].compName(ico) = medNameCompo;
01245                           }
01246                       }
01247                 }
01248             }
01249         }
01250     } // iterate on _listGIBItoMED_comp
01251 
01252   for ( size_t i = 0; i < _nodeFields.size() ; i++)
01253     usedNames.insert( _nodeFields[i]->_name );
01254   for ( size_t i = 0; i < _cellFields.size() ; i++)
01255     usedNames.insert( _cellFields[i]->_name );
01256 }
01257 
01258 //================================================================================
01262 //================================================================================
01263 
01264 void IntermediateMED::decreaseHierarchicalDepthOfSubgroups()
01265 {
01266   for (size_t i=0; i!=_groups.size(); ++i)
01267   {
01268     Group& grp = _groups[i];
01269     for (size_t j = 0; j < grp._groups.size(); ++j )
01270     {
01271       Group & sub_grp = *grp._groups[j];
01272       if ( !sub_grp._groups.empty() )
01273       {
01274         // replace j with its 1st subgroup
01275         grp._groups[j] = sub_grp._groups[0];
01276         // push back the rest subs
01277         grp._groups.insert( grp._groups.end(), ++sub_grp._groups.begin(), sub_grp._groups.end() );
01278       }
01279     }
01280     // remove empty sub-_groups
01281     vector< Group* > newSubGroups;
01282     newSubGroups.reserve( grp._groups.size() );
01283     for (size_t j = 0; j < grp._groups.size(); ++j )
01284       if ( !grp._groups[j]->empty() )
01285         newSubGroups.push_back( grp._groups[j] );
01286     if ( newSubGroups.size() < grp._groups.size() )
01287       grp._groups.swap( newSubGroups );
01288   }
01289 }
01290 
01291 //================================================================================
01295 //================================================================================
01296 
01297 void IntermediateMED::eraseUselessGroups()
01298 {
01299   // propagate _isProfile=true to sub-groups of composite groups
01300   // for (size_t int i=0; i!=_groups.size(); ++i)
01301   // {
01302   //   Group* grp = _groups[i];
01303   //   if ( grp->_isProfile && !grp->_groups.empty() )
01304   //     for (size_t j = 0; j < grp->_groups.size(); ++j )
01305   //       grp->_groups[j]->_isProfile=true;
01306   // }
01307   std::set<Group*> groups2convert;
01308   // keep not named sub-groups of field supports
01309   for (size_t i=0; i!=_groups.size(); ++i)
01310   {
01311     Group& grp = _groups[i];
01312     if ( grp._isProfile && !grp._groups.empty() )
01313       groups2convert.insert( grp._groups.begin(), grp._groups.end() );
01314   }
01315 
01316   // keep named groups and their subgroups
01317   for (size_t i=0; i!=_groups.size(); ++i)
01318   {
01319     Group& grp = _groups[i];
01320     if ( !grp._name.empty() && !grp.empty() )
01321     {
01322       groups2convert.insert( &grp );
01323       groups2convert.insert( grp._groups.begin(), grp._groups.end() );
01324     }
01325   }
01326   // erase groups that are not in groups2convert and not _isProfile
01327   for (size_t i=0; i!=_groups.size(); ++i)
01328   {
01329     Group* grp = &_groups[i];
01330     if ( !grp->_isProfile && !groups2convert.count( grp ) )
01331     {
01332       grp->_cells.clear();
01333       grp->_groups.clear();
01334     }
01335   }
01336 }
01337 
01338 //================================================================================
01342 //================================================================================
01343 
01344 void IntermediateMED::detectMixDimGroups()
01345 {
01346   //hasMixedCells = false;
01347   for ( size_t i=0; i < _groups.size(); ++i )
01348   {
01349     Group& grp = _groups[i];
01350     if ( grp._groups.size() < 2 )
01351       continue;
01352 
01353     // check if sub-groups have different dimension
01354     unsigned dim1 = getDim( &grp );
01355     for ( size_t j = 1; j  < grp._groups.size(); ++j )
01356     {
01357       unsigned dim2 = getDim( &grp );
01358       if ( dim1 != dim2 )
01359       {
01360         grp._cells.clear();
01361         grp._groups.clear();
01362         if ( !grp._name.empty() )
01363           cout << "Erase a group with elements of different dim |" << grp._name << "|"<< endl;
01364         break;
01365       }
01366     }
01367   }
01368 }
01369 
01370 //================================================================================
01374 //================================================================================
01375 
01376 void IntermediateMED::orientElements2D()
01377 {
01378   set<Cell>::const_iterator elemIt, elemEnd;
01379   vector< pair<int,int> > swapVec;
01380 
01381   // ------------------------------------
01382   // fix connectivity of quadratic edges
01383   // ------------------------------------
01384   set<Cell>& quadEdges = _cellsByType[ INTERP_KERNEL::NORM_SEG3 ];
01385   if ( !quadEdges.empty() )
01386     {
01387       elemIt = quadEdges.begin(), elemEnd = quadEdges.end();
01388       for ( ; elemIt != elemEnd; ++elemIt )
01389         ConvertQuadratic( INTERP_KERNEL::NORM_SEG3, *elemIt );
01390     }
01391 
01392   CellsByDimIterator faceIt( *this, 2 );
01393   while ( const set<Cell > * faces = faceIt.nextType() )
01394     {
01395       TCellType cellType = faceIt.type();
01396       bool isQuadratic = getGibi2MedQuadraticInterlace( cellType );
01397 
01398       getReverseVector( cellType, swapVec );
01399 
01400       // ------------------------------------
01401       // fix connectivity of quadratic faces
01402       // ------------------------------------
01403       if ( isQuadratic )
01404         for ( elemIt = faces->begin(), elemEnd = faces->end(); elemIt != elemEnd; elemIt++ )
01405           ConvertQuadratic( cellType, *elemIt );
01406 
01407       // --------------------------
01408       // orient faces clockwise
01409       // --------------------------
01410       int iQuad = isQuadratic ? 2 : 1;
01411       for ( elemIt = faces->begin(), elemEnd = faces->end(); elemIt != elemEnd; elemIt++ )
01412         {
01413           // look for index of the most left node
01414           int iLeft = 0, iNode, nbNodes = elemIt->_nodes.size() / iQuad;
01415           double x, minX = nodeCoords( elemIt->_nodes[0] )[0];
01416           for ( iNode = 1; iNode < nbNodes; ++iNode )
01417             if (( x = nodeCoords( elemIt->_nodes[ iNode ])[ 0 ]) < minX )
01418               minX = x, iLeft = iNode;
01419 
01420           // indeces of the nodes neighboring the most left one
01421           int iPrev = ( iLeft - 1 < 0 ) ? nbNodes - 1 : iLeft - 1;
01422           int iNext = ( iLeft + 1 == nbNodes ) ? 0 : iLeft + 1;
01423           // find components of prev-left and left-next vectors
01424           double xP = nodeCoords( elemIt->_nodes[ iPrev ])[ 0 ];
01425           double yP = nodeCoords( elemIt->_nodes[ iPrev ])[ 1 ];
01426           double xN = nodeCoords( elemIt->_nodes[ iNext ])[ 0 ];
01427           double yN = nodeCoords( elemIt->_nodes[ iNext ])[ 1 ];
01428           double xL = nodeCoords( elemIt->_nodes[ iLeft ])[ 0 ];
01429           double yL = nodeCoords( elemIt->_nodes[ iLeft ])[ 1 ];
01430           double xPL = xL - xP, yPL = yL - yP; // components of prev-left vector
01431           double xLN = xN - xL, yLN = yN - yL; // components of left-next vector
01432           // normalise y of the vectors
01433           double modPL = sqrt ( xPL * xPL + yPL * yPL );
01434           double modLN = sqrt ( xLN * xLN + yLN * yLN );
01435           if ( modLN > std::numeric_limits<double>::min() &&
01436                modPL > std::numeric_limits<double>::min() )
01437             {
01438               yPL /= modPL;
01439               yLN /= modLN;
01440               // summary direction of neighboring links must be positive
01441               bool clockwise = ( yPL + yLN > 0 );
01442               if ( !clockwise )
01443                 reverse( *elemIt, swapVec );
01444             }
01445         }
01446     }
01447 }
01448 
01449 //================================================================================
01453 //================================================================================
01454 
01455 void IntermediateMED::orientElements3D()
01456 {
01457   // set _reverse flags of faces
01458   orientFaces3D();
01459 
01460   // -----------------
01461   // fix connectivity
01462   // -----------------
01463 
01464   set<Cell>::const_iterator elemIt, elemEnd;
01465   vector< pair<int,int> > swapVec;
01466 
01467   for ( int dim = 1; dim <= 3; ++dim )
01468   {
01469     CellsByDimIterator cellsIt( *this, dim );
01470     while ( const set<Cell > * elems = cellsIt.nextType() )
01471     {
01472       TCellType cellType = cellsIt.type();
01473       bool isQuadratic = getGibi2MedQuadraticInterlace( cellType );
01474       getReverseVector( cellType, swapVec );
01475 
01476       elemIt = elems->begin(), elemEnd = elems->end();
01477       for ( ; elemIt != elemEnd; elemIt++ )
01478       {
01479         // GIBI connectivity -> MED one
01480         if( isQuadratic )
01481           ConvertQuadratic( cellType, *elemIt );
01482 
01483         // reverse faces
01484         if ( elemIt->_reverse )
01485           reverse ( *elemIt, swapVec );
01486       }
01487     }
01488   }
01489 
01490   orientVolumes();
01491 }
01492 
01493 //================================================================================
01497 //================================================================================
01498 
01499 void IntermediateMED::orientFaces3D()
01500 {
01501   // fill map of links and their faces
01502   set<const Cell*> faces;
01503   map<const Cell*, Group*> fgm;
01504   map<Link, list<const Cell*> > linkFacesMap;
01505   map<Link, list<const Cell*> >::iterator lfIt, lfIt2;
01506 
01507   for (size_t i=0; i!=_groups.size(); ++i)
01508     {
01509       Group& grp = _groups[i];
01510       if ( !grp._cells.empty() && getDimension( grp._cellType ) == 2 )
01511         for ( size_t j = 0; j < grp._cells.size(); ++j )
01512           if ( faces.insert( grp._cells[j] ).second )
01513             {
01514               for ( size_t k = 0; k < grp._cells[j]->_nodes.size(); ++k )
01515                 linkFacesMap[ grp._cells[j]->link( k ) ].push_back( grp._cells[j] );
01516               fgm.insert( make_pair( grp._cells[j], &grp ));
01517             }
01518     }
01519   // dump linkFacesMap
01520   //     for ( lfIt = linkFacesMap.begin(); lfIt!=linkFacesMap.end(); lfIt++) {
01521   //       cout<< "LINK: " << lfIt->first.first << "-" << lfIt->first.second << endl;
01522   //       list<const Cell*> & fList = lfIt->second;
01523   //       list<const Cell*>::iterator fIt = fList.begin();
01524   //       for ( ; fIt != fList.end(); fIt++ )
01525   //         cout << "\t" << **fIt << fgm[*fIt]->nom << endl;
01526   //     }
01527 
01528   // Each oriented link must appear in one face only, else a face is reversed.
01529 
01530   queue<const Cell*> faceQueue; /* the queue contains well oriented faces
01531                                      whose neighbors orientation is to be checked */
01532   bool manifold = true;
01533   while ( !linkFacesMap.empty() )
01534     {
01535       if ( faceQueue.empty() )
01536         {
01537           assert( !linkFacesMap.begin()->second.empty() );
01538           faceQueue.push( linkFacesMap.begin()->second.front() );
01539         }
01540       while ( !faceQueue.empty() )
01541         {
01542           const Cell* face = faceQueue.front();
01543           faceQueue.pop();
01544 
01545           // loop on links of <face>
01546           for ( int i = 0; i < (int)face->_nodes.size(); ++i )
01547             {
01548               Link link = face->link( i );
01549               // find the neighbor faces
01550               lfIt = linkFacesMap.find( link );
01551               int nbFaceByLink = 0;
01552               list< const Cell* > ml;
01553               if ( lfIt != linkFacesMap.end() )
01554                 {
01555                   list<const Cell*> & fList = lfIt->second;
01556                   list<const Cell*>::iterator fIt = fList.begin();
01557                   assert( fIt != fList.end() );
01558                   for ( ; fIt != fList.end(); fIt++, nbFaceByLink++ )
01559                     {
01560                       ml.push_back( *fIt );
01561                       if ( *fIt != face ) // wrongly oriented neighbor face
01562                         {
01563                           const Cell* badFace = *fIt;
01564                           // reverse and remove badFace from linkFacesMap
01565                           for ( int j = 0; j < (int)badFace->_nodes.size(); ++j )
01566                             {
01567                               Link badlink = badFace->link( j );
01568                               if ( badlink == link ) continue;
01569                               lfIt2 = linkFacesMap.find( badlink );
01570                               if ( lfIt2 != linkFacesMap.end() )
01571                                 {
01572                                   list<const Cell*> & ff = lfIt2->second;
01573                                   ff.erase( find( ff.begin(), ff.end(), badFace ));
01574                                   if ( ff.empty() )
01575                                     linkFacesMap.erase( lfIt2 );
01576                                 }
01577                             }
01578                           badFace->_reverse = true; // reverse
01579                           //INFOS_MED( "REVERSE " << *badFace );
01580                           faceQueue.push( badFace );
01581                         }
01582                     }
01583                   linkFacesMap.erase( lfIt );
01584                 }
01585               // add good neighbors to the queue
01586               Link revLink( link.second, link.first );
01587               lfIt = linkFacesMap.find( revLink );
01588               if ( lfIt != linkFacesMap.end() )
01589                 {
01590                   list<const Cell*> & fList = lfIt->second;
01591                   list<const Cell*>::iterator fIt = fList.begin();
01592                   for ( ; fIt != fList.end(); fIt++, nbFaceByLink++ )
01593                     {
01594                       ml.push_back( *fIt );
01595                       if ( *fIt != face )
01596                         faceQueue.push( *fIt );
01597                     }
01598                   linkFacesMap.erase( lfIt );
01599                 }
01600               if ( nbFaceByLink > 2 )
01601                 {
01602                   if ( manifold )
01603                     {
01604                       list<const Cell*>::iterator ii = ml.begin();
01605                       cout << nbFaceByLink << " faces by 1 link:";
01606                       for( ; ii!= ml.end(); ii++ )
01607                         cout << "in sub-mesh " << fgm[ *ii ]->_name << endl << **ii;
01608                     }
01609                   manifold = false;
01610                 }
01611             } // loop on links of the being checked face
01612         } // loop on the face queue
01613     } // while ( !linkFacesMap.empty() )
01614 
01615   if ( !manifold )
01616     cout << " -> Non manifold mesh, faces orientation may be incorrect" << endl;
01617 }
01618 
01619 //================================================================================
01624 //================================================================================
01625 
01626 void IntermediateMED::orientVolumes()
01627 {
01628   set<Cell>::const_iterator elemIt, elemEnd;
01629   vector< pair<int,int> > swapVec;
01630 
01631   CellsByDimIterator cellsIt( *this, 3 );
01632   while ( const set<Cell > * elems = cellsIt.nextType() )
01633     {
01634       TCellType cellType = cellsIt.type();
01635       elemIt = elems->begin(), elemEnd = elems->end();
01636       int nbBottomNodes = 0;
01637       switch ( cellType )
01638         {
01639         case NORM_TETRA4:
01640         case NORM_TETRA10:
01641         case NORM_PENTA6:
01642         case NORM_PENTA15:
01643           nbBottomNodes = 3; break;
01644         case NORM_PYRA5:
01645         case NORM_PYRA13:
01646         case NORM_HEXA8:
01647         case NORM_HEXA20:
01648           nbBottomNodes = 4; break;
01649         default: continue;
01650         }
01651       getReverseVector( cellType, swapVec );
01652 
01653       for ( ; elemIt != elemEnd; elemIt++ )
01654         {
01655           // find a normal to the bottom face
01656           const double* n[4];
01657           n[0] = nodeCoords( elemIt->_nodes[0]); // 3 bottom nodes
01658           n[1] = nodeCoords( elemIt->_nodes[1]);
01659           n[2] = nodeCoords( elemIt->_nodes[2]);
01660           n[3] = nodeCoords( elemIt->_nodes[nbBottomNodes]); // a top node
01661           double vec01[3]; // vector n[0]-n[1]
01662           vec01[0] = n[1][0] - n[0][0];
01663           vec01[1] = n[1][1] - n[0][1];
01664           vec01[2] = n[1][2] - n[0][2];
01665           double vec02 [3]; // vector n[0]-n[2]
01666           vec02[0] = n[2][0] - n[0][0];
01667           vec02[1] = n[2][1] - n[0][1];
01668           vec02[2] = n[2][2] - n[0][2];
01669           double normal [3]; // vec01 ^ vec02
01670           normal[0] = vec01[1] * vec02[2] - vec01[2] * vec02[1];
01671           normal[1] = vec01[2] * vec02[0] - vec01[0] * vec02[2];
01672           normal[2] = vec01[0] * vec02[1] - vec01[1] * vec02[0];
01673           // check if the 102 angle is convex
01674           if ( nbBottomNodes > 3 )
01675             {
01676               const double* n3 = nodeCoords( elemIt->_nodes[nbBottomNodes-1] );// last bottom node
01677               double vec03 [3];  // vector n[0]-n3
01678               vec03[0] = n3[0] - n[0][0];
01679               vec03[1] = n3[1] - n[0][1];
01680               vec03[2] = n3[2] - n[0][2];
01681               if ( fabs( normal[0]+normal[1]+normal[2] ) <= numeric_limits<double>::max() ) // vec01 || vec02
01682                 {
01683                   normal[0] = vec01[1] * vec03[2] - vec01[2] * vec03[1]; // vec01 ^ vec03
01684                   normal[1] = vec01[2] * vec03[0] - vec01[0] * vec03[2];
01685                   normal[2] = vec01[0] * vec03[1] - vec01[1] * vec03[0];
01686                 }
01687               else
01688                 {
01689                   double vec [3]; // normal ^ vec01
01690                   vec[0] = normal[1] * vec01[2] - normal[2] * vec01[1];
01691                   vec[1] = normal[2] * vec01[0] - normal[0] * vec01[2];
01692                   vec[2] = normal[0] * vec01[1] - normal[1] * vec01[0];
01693                   double dot2 = vec[0]*vec03[0] + vec[1]*vec03[1] + vec[2]*vec03[2]; // vec*vec03
01694                   if ( dot2 < 0 ) // concave -> reverse normal
01695                     {
01696                       normal[0] *= -1;
01697                       normal[1] *= -1;
01698                       normal[2] *= -1;
01699                     }
01700                 }
01701             }
01702           // direction from top to bottom
01703           double tbDir[3];
01704           tbDir[0] = n[0][0] - n[3][0];
01705           tbDir[1] = n[0][1] - n[3][1];
01706           tbDir[2] = n[0][2] - n[3][2];
01707 
01708           // compare 2 directions: normal and top-bottom
01709           double dot = normal[0]*tbDir[0] + normal[1]*tbDir[1] + normal[2]*tbDir[2];
01710           if ( dot < 0. ) // need reverse
01711             reverse( *elemIt, swapVec );
01712 
01713         } // loop on volumes of one geometry
01714     } // loop on 3D geometry types
01715 
01716 }
01717 
01718 //================================================================================
01722 //================================================================================
01723 
01724 int NodeContainer::numberNodes()
01725 {
01726   int id = 1;
01727   for ( size_t i = 0; i < _nodes.size(); ++i )
01728     for ( size_t j = 0; j < _nodes[i].size(); ++j )
01729       if ( _nodes[i][j].isUsed() )
01730         _nodes[i][j]._number = id++;
01731   return id-1;
01732 }
01733 
01734 
01735 //================================================================================
01739 //================================================================================
01740 
01741 void IntermediateMED::numberElements()
01742 {
01743   set<Cell>::const_iterator elemIt, elemEnd;
01744 
01745   // numbering _cells of type NORM_POINT1 by node number
01746   {
01747     const set<Cell>& points = _cellsByType[ INTERP_KERNEL::NORM_POINT1 ];
01748     elemIt = points.begin(), elemEnd = points.end();
01749     for ( ; elemIt != elemEnd; ++elemIt )
01750       elemIt->_number = elemIt->_nodes[0]->_number;
01751   }
01752 
01753   // numbering 1D-3D _cells
01754   for ( int dim = 1; dim <= 3; ++dim )
01755     {
01756       // check if re-numeration is needed (to try to keep elem oreder as in sauve file )
01757       bool ok = true, renumEntity = false;
01758       CellsByDimIterator cellsIt( *this, dim );
01759       int prevNbElems = 0;
01760       while ( const set<Cell> * typeCells = cellsIt.nextType() )
01761         {
01762           TID minNumber = std::numeric_limits<TID>::max(), maxNumber = 0;
01763           for ( elemIt = typeCells->begin(), elemEnd = typeCells->end(); elemIt!=elemEnd; ++elemIt)
01764             {
01765               if ( elemIt->_number < minNumber ) minNumber = elemIt->_number;
01766               if ( elemIt->_number > maxNumber ) maxNumber = elemIt->_number;
01767             }
01768           TID typeSize = typeCells->size();
01769           if ( typeSize != maxNumber - minNumber + 1 )
01770             ok = false;
01771           if ( prevNbElems != 0 ) {
01772             if ( minNumber == 1 )
01773               renumEntity = true;
01774             else if ( prevNbElems+1 != (int)minNumber )
01775               ok = false;
01776           }
01777           prevNbElems += typeSize;
01778         }
01779 
01780       if ( ok && renumEntity ) // each geom type was numerated separately
01781         {
01782           cellsIt.init( dim );
01783           prevNbElems = cellsIt.nextType()->size(); // no need to renumber the first type
01784           while ( const set<Cell> * typeCells = cellsIt.nextType() )
01785             {
01786               for ( elemIt = typeCells->begin(), elemEnd = typeCells->end(); elemIt!=elemEnd; ++elemIt)
01787                 elemIt->_number += prevNbElems;
01788               prevNbElems += typeCells->size();
01789             }
01790         }
01791       if ( !ok )
01792         {
01793           int cellID=1;
01794           cellsIt.init( dim );
01795           while ( const set<Cell> * typeCells = cellsIt.nextType() )
01796             for ( elemIt = typeCells->begin(), elemEnd = typeCells->end(); elemIt!=elemEnd; ++elemIt)
01797               elemIt->_number = cellID++;
01798         }
01799     }
01800 }
01801 
01802 //================================================================================
01806 //================================================================================
01807 
01808 ParaMEDMEM::DataArrayDouble * IntermediateMED::getCoords()
01809 {
01810   DataArrayDouble* coordArray = DataArrayDouble::New();
01811   coordArray->alloc( _nbNodes, _spaceDim );
01812   double * coordPrt = coordArray->getPointer();
01813   for ( int i = 0, nb = _points.size(); i < nb; ++i )
01814     {
01815       Node* n = getNode( i+1 );
01816       if ( n->isUsed() )
01817         {
01818           const double* nCoords = nodeCoords( n );
01819           std::copy( nCoords, nCoords+_spaceDim, coordPrt );
01820           coordPrt += _spaceDim;
01821         }
01822     }
01823   return coordArray;
01824 }
01825 
01826 //================================================================================
01832 //================================================================================
01833 
01834 void IntermediateMED::setConnectivity( ParaMEDMEM::MEDFileUMesh*    mesh,
01835                                        ParaMEDMEM::DataArrayDouble* coords )
01836 {
01837   int meshDim = 0;
01838 
01839   mesh->setCoords( coords );
01840 
01841   set<Cell>::const_iterator elemIt, elemEnd;
01842   for ( int dim = 3; dim > 0; --dim )
01843   {
01844     CellsByDimIterator dimCells( *this, dim );
01845 
01846     int nbOfCells = 0;
01847     while ( const std::set<Cell > * cells = dimCells.nextType() )
01848       nbOfCells += cells->size();
01849     if ( nbOfCells == 0 )
01850       continue;
01851 
01852     if ( !meshDim ) meshDim = dim;
01853 
01854     MEDCouplingUMesh* dimMesh = MEDCouplingUMesh::New();
01855     dimMesh->setCoords( coords );
01856     dimMesh->setMeshDimension( dim );
01857     dimMesh->allocateCells( nbOfCells );
01858 
01859     int prevNbCells = 0;
01860     dimCells.init( dim );
01861     while ( const std::set<Cell > * cells = dimCells.nextType() )
01862       {
01863         // fill connectivity array to take into account order of elements in the sauv file
01864         const int nbCellNodes = cells->begin()->_nodes.size();
01865         vector< TID > connectivity( cells->size() * nbCellNodes );
01866         int * nodalConnOfCell;
01867         for ( elemIt = cells->begin(), elemEnd = cells->end(); elemIt != elemEnd; ++elemIt )
01868           {
01869             const Cell& cell = *elemIt;
01870             const int index = cell._number - 1 - prevNbCells;
01871             nodalConnOfCell = &connectivity[ index * nbCellNodes ];
01872             if ( cell._reverse )
01873               for ( int i = nbCellNodes-1; i >= 0; --i )
01874                 *nodalConnOfCell++ = cell._nodes[i]->_number - 1;
01875             else
01876               for ( int i = 0; i < nbCellNodes; ++i )
01877                 *nodalConnOfCell++ = cell._nodes[i]->_number - 1;
01878           }
01879         prevNbCells += cells->size();
01880 
01881         // fill dimMesh
01882         TCellType cellType = dimCells.type();
01883         nodalConnOfCell = &connectivity[0];
01884         for ( size_t i = 0; i < cells->size(); ++i, nodalConnOfCell += nbCellNodes )
01885           dimMesh->insertNextCell( cellType, nbCellNodes, nodalConnOfCell );
01886       }
01887     dimMesh->finishInsertingCells();
01888     mesh->setMeshAtLevel( dim - meshDim, dimMesh );
01889     dimMesh->decrRef();
01890   }
01891 }
01892 
01893 //================================================================================
01898 //================================================================================
01899 
01900 void IntermediateMED::setGroups( ParaMEDMEM::MEDFileUMesh* mesh )
01901 {
01902   bool isMeshNameSet = false;
01903   const int meshDim = mesh->getMeshDimension();
01904   for ( int dim = 0; dim <= meshDim; ++dim )
01905     {
01906       const int meshDimRelToMaxExt = ( dim == 0 ? 1 : dim - meshDim );
01907 
01908       vector<const DataArrayInt *> medGroups;
01909       vector<MEDCouplingAutoRefCountObjectPtr<DataArrayInt> > refGroups;
01910       for ( size_t i = 0; i < _groups.size(); ++i )
01911         {
01912           Group& grp = _groups[i];
01913           if ( (int)getDim( &grp ) != dim )
01914             continue;
01915           // convert only named groups or field supports
01916           if ( grp.empty() || (grp._name.empty() && !grp._isProfile ))
01917             continue;
01918           //if ( grp._medGroup ) continue; // already converted
01919 
01920           // sort cells by ID and remember their initial order in the group
01921           TCellToOrderMap cell2order;
01922           unsigned orderInGroup = 0;
01923           vector< Group* > groupVec;
01924           if ( grp._groups.empty() ) groupVec.push_back( & grp );
01925           else                       groupVec = grp._groups;
01926           for ( size_t iG = 0; iG < groupVec.size(); ++iG )
01927             {
01928               Group* aG = groupVec[ iG ];
01929               for ( size_t iC = 0; iC < aG->_cells.size(); ++iC )
01930                 cell2order.insert( cell2order.end(), make_pair( aG->_cells[iC], orderInGroup++ ));
01931             }
01932           bool isSelfIntersect = ( orderInGroup != cell2order.size() );
01933           if ( isSelfIntersect ) // self intersecting group
01934             {
01935               ostringstream msg;
01936               msg << "Self intersecting sub-mesh: id = " << i+1
01937                   << ", name = |" << grp._name << "|" << endl
01938                   << " nb unique elements = " << cell2order.size() << endl
01939                   << " total nb elements  = " << orderInGroup;
01940               if ( grp._isProfile )
01941                 {
01942                   THROW_IK_EXCEPTION( msg.str() );
01943                 }
01944               else
01945                 {
01946                   cout << msg << endl;
01947                 }
01948             }
01949           // create a med group
01950           grp._medGroup = DataArrayInt::New();
01951           grp._medGroup->setName( grp._name.c_str() );
01952           grp._medGroup->alloc( orderInGroup, /*nbOfCompo=*/1 );
01953           int * idsPrt = grp._medGroup->getPointer();
01954           TCellToOrderMap::iterator cell2orderIt, cell2orderEnd = cell2order.end();
01955           for ( cell2orderIt = cell2order.begin(); cell2orderIt != cell2orderEnd; ++cell2orderIt )
01956             *idsPrt++ = (*cell2orderIt).first->_number - 1;
01957 
01958           // try to set the mesh name
01959           if ( !isMeshNameSet &&
01960                dim == meshDim &&
01961                !grp._name.empty() &&
01962                grp.size() == mesh->getSizeAtLevel( meshDimRelToMaxExt ))
01963             {
01964               mesh->setName( grp._name.c_str() );
01965               isMeshNameSet = true;
01966             }
01967           else if ( !grp._name.empty() )
01968             {
01969               medGroups.push_back( grp._medGroup );
01970             }
01971           // set relocation table
01972           setRelocationTable( &grp, cell2order );
01973 
01974           // Issue 0021311. Use case: a gibi group has references (recorded in pile 1)
01975           // and several names (pile 27) refer (pile 10) to this group.
01976           // We create a copy of this group per each named reference
01977           for ( unsigned iRef = 0 ; iRef < grp._refNames.size(); ++iRef )
01978             if ( !grp._refNames[ iRef ].empty() )
01979               {
01980                 refGroups.push_back( grp._medGroup->deepCpy() );
01981                 refGroups.back()->setName( grp._refNames[ iRef ].c_str() );
01982                 medGroups.push_back( refGroups.back() );
01983               }
01984         }
01985       mesh->setGroupsAtLevel( meshDimRelToMaxExt, medGroups );
01986     }
01987 }
01988 
01989 //================================================================================
01993 //================================================================================
01994 
01995 bool IntermediateMED::isOnAll( const Group* grp, int & dimRel ) const
01996 {
01997   int dim = getDim( grp );
01998 
01999   int nbElems = 0;
02000   CellsByDimIterator dimCells( *this, dim );
02001   while ( const set<Cell > * cells = dimCells.nextType() )
02002     nbElems += cells->size();
02003 
02004   const bool onAll = ( nbElems == grp->size() );
02005 
02006   if ( dim == 0 )
02007     dimRel = 0;
02008   else
02009     {
02010       int meshDim = 3;
02011       for ( ; meshDim > 0; --meshDim )
02012         {
02013           dimCells.init( meshDim );
02014           if ( dimCells.nextType() )
02015             break;
02016         }
02017       dimRel = dim - meshDim;
02018     }
02019   return onAll;
02020 }
02021 
02022 //================================================================================
02026 //================================================================================
02027 
02028 ParaMEDMEM::MEDFileFields * IntermediateMED::makeMEDFileFields(ParaMEDMEM::MEDFileUMesh* mesh)
02029 {
02030   if ( _nodeFields.empty() && _cellFields.empty() ) return 0;
02031 
02032   // set long names
02033   set< string > usedFieldNames;
02034   setFieldLongNames(usedFieldNames);
02035 
02036   MEDFileFields* fields = MEDFileFields::New();
02037 
02038   for ( size_t i = 0; i < _nodeFields.size(); ++i )
02039     setFields( _nodeFields[i], fields, mesh, i+1, usedFieldNames );
02040 
02041   for ( size_t i = 0; i < _cellFields.size(); ++i )
02042     setFields( _cellFields[i], fields, mesh, i+1, usedFieldNames );
02043 
02044   return fields;
02045 }
02046 
02047 //================================================================================
02051 //================================================================================
02052 
02053 void IntermediateMED::setFields( SauvUtilities::DoubleField* fld,
02054                                  ParaMEDMEM::MEDFileFields*  medFields,
02055                                  ParaMEDMEM::MEDFileUMesh*   mesh,
02056                                  const TID                   castemID,
02057                                  set< string >&              usedFieldNames)
02058 {
02059   if ( !fld || !fld->isMedCompatible() ) return;
02060 
02061   // if ( !fld->hasCommonSupport() ):
02062   //     each sub makes MEDFileFieldMultiTS
02063   // else:
02064   //     unite several subs into a MEDCouplingFieldDouble
02065 
02066   const bool uniteSubs = fld->hasCommonSupport();
02067   if ( !uniteSubs )
02068     cout << "Castem field #" << castemID << " " << fld->_name
02069          << " is incompatible with MED format, so we split it into several fields" << endl;
02070 
02071   for ( size_t iSub = 0; iSub < fld->_sub.size(); )
02072     {
02073       // set field name
02074       if ( !uniteSubs || fld->_name.empty() )
02075         makeFieldNewName( usedFieldNames, fld );
02076 
02077       // allocate values
02078       DataArrayDouble * values = DataArrayDouble::New();
02079       values->alloc( fld->getNbTuples(iSub), fld->_sub[iSub].nbComponents() );
02080 
02081       // set values
02082       double * valPtr = values->getPointer();
02083       if ( uniteSubs )
02084         {
02085           int nbElems = fld->_group->size();
02086           for ( int elemShift = 0; elemShift < nbElems; )
02087             elemShift += fld->setValues( valPtr, iSub++, elemShift );
02088           setTS( fld, values, medFields, mesh );
02089         }
02090       else
02091         {
02092           fld->setValues( valPtr, iSub++ );
02093           setTS( fld, values, medFields, mesh, iSub );
02094         }
02095     }
02096 }
02097 
02098 //================================================================================
02102 //================================================================================
02103 
02104 void IntermediateMED::setTS( SauvUtilities::DoubleField*  fld,
02105                              ParaMEDMEM::DataArrayDouble* values,
02106                              ParaMEDMEM::MEDFileFields*   medFields,
02107                              ParaMEDMEM::MEDFileUMesh*    mesh,
02108                              const int                    iSub)
02109 {
02110   // analyze a field support
02111   const Group* support = fld->getSupport();
02112   int dimRel;
02113   const bool onAll = isOnAll( support, dimRel );
02114   if ( !onAll && support->_name.empty() )
02115     {
02116       const_cast<Group*>(support)->_name += "PFL_" + fld->_name;
02117       support->_medGroup->setName( support->_name.c_str() );
02118     }
02119 
02120   // make a time-stamp
02121   MEDCouplingFieldDouble * timeStamp = MEDCouplingFieldDouble::New( fld->getMedType(),
02122                                                                     fld->getMedTimeDisc() );
02123   timeStamp->setName( fld->_name.c_str() );
02124   timeStamp->setDescription( fld->_description.c_str() );
02125   MEDCouplingAutoRefCountObjectPtr< MEDCouplingUMesh > dimMesh = mesh->getMeshAtLevel( dimRel );
02126   timeStamp->setMesh( dimMesh );
02127   for ( size_t i = 0; i < (size_t)fld->_sub[iSub].nbComponents(); ++i )
02128     values->setInfoOnComponent( i, fld->_sub[iSub]._comp_names[ i ].c_str() );
02129   timeStamp->setArray( values );
02130   values->decrRef();
02131 
02132   // get a field to add the time-stamp
02133   bool isNewMedField = false;
02134   if ( !fld->_curMedField || fld->_name != fld->_curMedField->getName() )
02135     {
02136       fld->_curMedField = MEDFileFieldMultiTS::New();
02137       isNewMedField = true;
02138     }
02139 
02140   // set an order
02141   timeStamp->setOrder( fld->_curMedField->getNumberOfTS() );
02142 
02143   // add the time-stamp
02144   if ( onAll )
02145     fld->_curMedField->appendFieldNoProfileSBT( timeStamp );
02146   else
02147     fld->_curMedField->appendFieldProfile( timeStamp, mesh, dimRel, support->_medGroup );
02148   timeStamp->decrRef();
02149 
02150   if ( isNewMedField ) // timeStamp must be added before this
02151     {
02152       medFields->pushField( fld->_curMedField );
02153       fld->_curMedField->decrRef();
02154     }
02155 }
02156 
02157 //================================================================================
02161 //================================================================================
02162 
02163 void IntermediateMED::makeFieldNewName(std::set< std::string >&    usedNames,
02164                                        SauvUtilities::DoubleField* fld )
02165 {
02166   string base = fld->_name;
02167   if ( base.empty() )
02168     {
02169       base = "F_";
02170     }
02171   else
02172     {
02173       string::size_type pos = base.rfind('_');
02174       if ( pos != string::npos )
02175         base = base.substr( 0, pos+1 );
02176       else
02177         base += '_';
02178     }
02179 
02180   int i = 1;
02181   do
02182     {
02183       fld->_name = base + SauvUtilities::toString( i++ );
02184     }
02185   while( !usedNames.insert( fld->_name ).second );
02186 }
02187 
02188 //================================================================================
02192 //================================================================================
02193 
02194 std::vector< double >& DoubleField::addComponent( int nb_values )
02195 {
02196   _comp_values.push_back( std::vector< double >() );
02197   std::vector< double >& res = _comp_values.back();
02198   res.resize( nb_values );
02199   return res;
02200 }
02201 //================================================================================
02205 //================================================================================
02206 
02207 const Group* DoubleField::getSupport( const int iSub ) const
02208 {
02209   return _group ? _group : _sub[iSub]._support;
02210 }
02211 
02212 //================================================================================
02216 //================================================================================
02217 
02218 bool DoubleField::isMultiTimeStamps() const
02219 {
02220   if ( _sub.size() < 2 )
02221     return false;
02222   bool sameSupports = true;
02223   Group* grp1 = _sub[0]._support;
02224   for ( size_t i = 1; i < _sub.size() && sameSupports; ++i )
02225     sameSupports = ( grp1 == _sub[i]._support );
02226 
02227   return sameSupports;
02228 }
02229 
02230 //================================================================================
02234 //================================================================================
02235 
02236 bool DoubleField::isMedCompatible() const
02237 {
02238   for ( size_t iSub = 0; iSub < _sub.size(); ++iSub )
02239     {
02240       if ( !getSupport(iSub) || !getSupport(iSub)->_medGroup )
02241         THROW_IK_EXCEPTION("SauvReader INTERNAL ERROR: NULL field support");
02242 
02243       if ( !_sub[iSub].isValidNbGauss() )
02244         {
02245           cout << "Skip field <" << _name << "> : different nb of gauss points in components" <<endl;
02246           return false;
02247         }
02248     }
02249   // check if there are no gauss or nbGauss() == nbCellNodes,
02250   // else we lack info on gauss point localization
02251   // TODO?
02252   return true;
02253 }
02254 
02255 //================================================================================
02259 //================================================================================
02260 
02261 bool DoubleField::hasSameComponentsBySupport() const
02262 {
02263   vector< _Sub_data >::const_iterator sub_data = _sub.begin();
02264   const _Sub_data& first_sub_data = *sub_data;
02265   for ( ++sub_data ; sub_data != _sub.end(); ++sub_data )
02266   {
02267     if ( first_sub_data._comp_names != sub_data->_comp_names )
02268       return false; // diff names of components
02269 
02270     if ( first_sub_data._nb_gauss != sub_data->_nb_gauss &&
02271          first_sub_data._support->_cellType == sub_data->_support->_cellType)
02272       return false; // diff nb of gauss points on same cell type
02273   }
02274   return true;
02275 }
02276 
02277 //================================================================================
02281 //================================================================================
02282 
02283 ParaMEDMEM::TypeOfField DoubleField::getMedType( const int iSub ) const
02284 {
02285   using namespace INTERP_KERNEL;
02286 
02287   const Group* grp = hasCommonSupport() ? _group : _sub[iSub]._support;
02288   if ( _sub[iSub].nbGauss() > 1 )
02289     {
02290       const CellModel& cm = CellModel::GetCellModel( _sub[iSub]._support->_cellType );
02291       return (int) cm.getNumberOfNodes() == _sub[iSub].nbGauss() ? ON_GAUSS_NE : ON_GAUSS_PT;
02292     }
02293   else
02294     {
02295       return getDim( grp ) == 0 ? ON_NODES : ON_CELLS;
02296     }
02297 }
02298 
02299 //================================================================================
02303 //================================================================================
02304 
02305 ParaMEDMEM::TypeOfTimeDiscretization DoubleField::getMedTimeDisc() const
02306 {
02307   return ONE_TIME;
02308   // NO_TIME = 4,
02309   // ONE_TIME = 5,
02310   // LINEAR_TIME = 6,
02311   // CONST_ON_TIME_INTERVAL = 7
02312 }
02313 
02314 //================================================================================
02318 //================================================================================
02319 
02320 int DoubleField::getNbTuples( const int iSub ) const
02321 {
02322   int nb = 0;
02323   if ( hasCommonSupport() && !_group->_groups.empty() )
02324     for ( size_t i = 0; i < _group->_groups.size(); ++i )
02325       nb += _sub[i].nbGauss() * _sub[i]._support->size();
02326   else
02327     nb = _sub[iSub].nbGauss() * getSupport(iSub)->size();
02328   return nb;
02329 }
02330 
02331 //================================================================================
02335 //================================================================================
02336 
02337 int DoubleField::setValues( double * valPtr, const int iSub, const int elemShift ) const
02338 {
02339   // find values for iSub
02340   int iComp = 0;
02341   for ( int iS = 0; iS < iSub; ++iS )
02342     iComp += _sub[iS].nbComponents();
02343   const vector< double > * compValues = &_comp_values[ iComp ];
02344 
02345   const vector< unsigned >& relocTable = getSupport( iSub )->_relocTable;
02346 
02347   // Set values
02348 
02349   const int nbElems      = _sub[iSub]._support->size();
02350   const int nbGauss      = _sub[iSub].nbGauss();
02351   const int nbComponents = _sub[iSub].nbComponents();
02352   const int nbValsByElem = nbComponents * nbGauss;
02353   // check nb values
02354   int nbVals = 0;
02355   for ( iComp = 0; iComp < nbComponents; ++iComp )
02356     nbVals += compValues[iComp].size();
02357   if ( nbVals != nbElems * nbValsByElem )
02358     THROW_IK_EXCEPTION("SauvMedConvertor.cxx: support size mismatches field size");
02359   // compute nb values in previous subs
02360   int valsShift = 0;
02361   for ( int iS = iSub-1, shift = elemShift; shift > 0; )
02362   {
02363     int nbE = _sub[iS]._support->size();
02364     shift -= nbE;
02365     valsShift += nbE * _sub[iS].nbComponents() * _sub[iS].nbGauss();
02366   }
02367   for ( int iE = 0; iE < nbElems; ++iE )
02368     {
02369       int iMed = valsShift + nbValsByElem * ( relocTable.empty() ? iE : relocTable[iE+elemShift]-elemShift );
02370       for ( iComp = 0; iComp < nbComponents; ++iComp )
02371         for ( int iG = 0; iG < nbGauss; ++iG )
02372           valPtr[ iMed + iG * nbComponents + iComp ] = compValues[iComp][ iE * nbGauss + iG ];
02373     }
02374   return nbElems;
02375 }
02376 
02377 //================================================================================
02381 //================================================================================
02382 
02383 IntermediateMED::~IntermediateMED()
02384 {
02385   for ( size_t i = 0; i < _nodeFields.size(); ++i )
02386     if ( _nodeFields[i] )
02387       delete _nodeFields[i];
02388   _nodeFields.clear();
02389 
02390   for ( size_t i = 0; i < _cellFields.size(); ++i )
02391     if ( _cellFields[i] )
02392       delete _cellFields[i];
02393   _cellFields.clear();
02394 
02395   for ( size_t i = 0; i < _groups.size(); ++i )
02396     if ( _groups[i]._medGroup )
02397       _groups[i]._medGroup->decrRef();
02398 }
02399 
02400 //================================================================================
02404 CellsByDimIterator::CellsByDimIterator( const IntermediateMED & medi, int dimm)
02405 {
02406   myImed = & medi;
02407   init( dimm );
02408 }
02412 void CellsByDimIterator::init(const int  dimm)
02413 {
02414   myCurType = -1;
02415   myTypeEnd = INTERP_KERNEL::NORM_HEXA20 + 1;
02416   myDim = dimm;
02417 }
02421 const std::set< Cell > * CellsByDimIterator::nextType()
02422 {
02423   while ( ++myCurType < myTypeEnd )
02424     if ( !myImed->_cellsByType[myCurType].empty() && ( myDim < 0 || dim(false) == myDim ))
02425       return & myImed->_cellsByType[myCurType];
02426   return 0;
02427 }
02431 int CellsByDimIterator::dim(const bool last) const
02432 {
02433   int typp = myCurType;
02434   if ( !last )
02435     while ( typp < myTypeEnd && myImed->_cellsByType[typp].empty() )
02436       ++typp;
02437   return typp < myTypeEnd ? getDimension( TCellType( typp )) : 4;
02438 }
02439 // END CellsByDimIterator ========================================================
02440